Hi,
As maybe some of you remember we had a long thread about stretched arcs
a while back.
I made erroneously an assumption with respect to the outside/inside of an arc
as Karl pointed out when I thought it was possible with an algebraic solution.
I gave the problem a good workout and came up with an algorithm that works
most of the time. The solution I went for is along the line trying to find the
normal from the arc that intersects the point and then check distance.
This is no easy task and my simple approach has problems with arcs that are
heavy stretched and extends barley over the narrow side.
Attached is a pdf doc (InsidePointProblem.pdf) explaining the problem which
occurs for points located on the inside of the arc and on the other side of the
axis from where the normal to the arc that intersects the point is.
Since I am not a geda developer I made the testing in Octave (and also Matlab
for better debugging capabilities)
I am attaching the octave files where the main one is:
isPointOnArc.m
which is the actual function and this file contains all help functions.
This file has been written using simple octave syntax that almost directly
translates to C-code, e.g. no fancy matrices operations and no usage of built
in functions. The other two are part of the testsuite with Test_isPointOnArc.m
the main test function, these does not translate to C-code as easily
(but almost.)
I am leaving this at it is since this was an itch I needed to scratch and I am
satisfied for the moment. The problem that still exists is something I might
have a look at later.
I dare say that this algorithm is better with respect to determining
intersection than todays pcb code but slightly more computational expensive.
It is a lot of code but not very heavy.
If for some reason some geda developer were to use this and need to ask
questions I will definitively try to answer. I have made a lot of comments
but it is most likely not enough.
B.R
/Igor
Attachment:
InsidePointProblem.pdf
Description: Adobe PDF document
function res = isPointOnArc(arc, point, R)
% Checks if the point is intersecting the given arc
% within the distance given by radius, e.g it is the
% same check if a circle with radius R is intersecting
% the arc.
%
% The function returns:
% 1 (true) if point is interesecting
% 0 (false) if point is not intersecting
% a negative number indicating something is wrong.
% -1 did not find normal to point on the outside of the arc
% even if the point is within arc extension.
% arc is a struct with members
% X
% Y
% Width
% Height
% Thickness
% StartAng
% Delta
%
% point is a struct with members
% X
% Y
%
% R is the max distance from the point for which an
% intersect can happen, e.g if the point is located
% at a distance larger then R from the arc it is not
% intersecting
res = 0;
a = arc.Width;
a2 = a*a;
b = arc.Height;
b2 = b*b;
c = arc.Thickness/2;
r = c+R;
r2 = r*r;
% StartAng and Delta are the so called parameterised
% angles, e.g does not correlate to the true physical
startAng = pi*arc.StartAng/180;
endAng = startAng + pi*arc.Delta/180;
% translate point to arc coordinates
x = point.X - arc.X;
x2 = x*x;
y = point.Y - arc.Y;
y2 = y*y;
% check if point is located outside the corresponding
% ellipse major bounding box.
if (x < -a-r) || (x > a+r) || (y < -b-r) || (y > b+r)
return;
end
% OK point is inside major box, start investigate
% Check if point is inside or outside arc center line
toCenterVal = 1 - x2/a2 - y2/b2;
if toCenterVal > 0
pointInside = 1; % point is located inside the centerline
elseif toCenterVal < 0
pointInside = 0; % point is located outside the centerline
else
pointInside = -1; % point is located on the centerline
end
% Check arc endpoints
startX = a*cos(startAng);
startY = b*sin(startAng);
px = x-startX;
py = y-startY;
if px*px+py*py <= r2
res = 1; % point is inside arc startpoint circle
return;
end
endX = a*cos(endAng);
endY = b*sin(endAng);
px = x-endX;
py = y-endY;
if px*px+py*py <= r2
res = 1; % point is inside arc endpoint circle
return;
end
% continue with the rest of the arc.
% we know that the function checkQuadrant returns a number
% from a list depending on placement of the point where the list
% is as:
% 0 : centerpoint
% 1 : positive x-axis
% 2 : first quadrant
% 3 : positive y-axis
% 4 : second quadrant
% 5 : negative x-axis
% 6 : third quadrant
% 7 : negative y-axis
% 8 : fourth quadrant
pointQuadrant = checkQuadrant(x,y);
pointOnX = 0;
pointOnY = 0;
if 0 == pointQuadrant
% point is at center position, enough to check if radius
% is larger than smallest of Heigth, Width
if (R>a) || (R>b)
res = 1;
else
res = 0;
end
return;
elseif (1== pointQuadrant) || (5 == pointQuadrant)
% point is on xaxis
pointOnX = 1;
elseif (3 == pointQuadrant) || (7 == pointQuadrant)
% point is on yaxis
pointOnY = 1;
end
% point is not at center, get physical angles
arcStartAlpha = getPointAngle(startX, startY);
arcEndAlpha = getPointAngle(endX, endY);
pointAlpha = getPointAngle(x,y);
% get arc quadrants
arcStartQuadrant = checkQuadrant(startX,startY);
arcEndQuadrant = checkQuadrant(endX,endY);
% check outside or centerline location
if pointInside <= 0
% point is outside or on centerline
direction = 1;
pointCovered = checkArcCoverage4Point( ...
a, b, x, y, arcStartQuadrant, arcEndQuadrant, pointQuadrant, ...
startX, startY, startAng, endX, endY, endAng, direction);
% if point is not coverd by arc it cant intersect
if 0 == pointCovered
return;
end
% If the point is on either axis and within the arc extension
% it is bound to intersect since it is inside major box which
% touches at the axis crossings but on or outside the centerline
if pointOnX || pointOnY
res = 1;
return;
end
% OK now we need to find out where the normal to the point
% is located on the arc. We do this by solving the equation
% (Bx-Ax)*(y-Ay)-(By-Ay)*(x-Ax) which is 0 on the normal
% and positive ccw of normal and negative cw of normal
% This is done by making a start guess of the parameterised
% ellipse angle where the normal is located and use this guess
% and evaluating the value and the corresponding derivative.
% Solution is said to be found if the resulting value is within
% given tolerance which is in % and based on the ellipse
% thickness and search radius. This means that for special
% case of finding a point with zero search radius on arc
% with zero thickness the point must fullfill the ellipse
% equation and be inside the arc extension.
if -1==pointInside
res = 1; % point is on arc center
return;
else
% point is outside ellipse center start searching.
[found, currVal] = search4Normal(a, b, x, y, pointAlpha, ...
pointQuadrant, r, direction);
if found
% we have found the normal to the point, calculate distance
ax = a*cos(currVal);
ay = b*sin(currVal);
d2 = (x-ax)^2+(y-ay)^2;
if d2 <= r2
res = 1;
else
res = 0;
end
return;
else
% We did not find normal, this is an error that should be catched
res = -1; % not found normal for point outside arc.
return;
end
end % else for: if -1==pointInside
else % if pointInside <= 0
% point is inside arc centerline
% There is a particular problem using the same type of algorithm for
% the point inside the arc since at the long end depending on the ratio
% between the two axises and the thickness+search radius there will be
% a swithch of the end of the inside normal with length half thicknes +
% search radius where it will cross the axis and be located in another
% quadrant than where the start of the normal is located. This
% phenomena will break the current search algorithm wich is based on
% point location, e.g if the point is located on the switched side and
% would have been covered by the sweep of the normal from the other
% quadrant but the arc does not extend far enough into the quadrant
% where the point is located it will not be seen as intersecting even
% if it is. This must be handled as a special case for inside points
% but I am for the moment lost on how to implement that. If it was
% possible to solve for where the end of the inside normal actually
% crosses the axis one can make check if the arc start/end is within
% this region and have the search try to find a normal touching the
% point as final step since this is a particular case which will not
% happen so often that only do this as last step if search has not
% given a true/false on intersection earlier.
% point is inside of arc
direction = 0;
pointCovered = checkArcCoverage4Point( ...
a, b, x, y, arcStartQuadrant, arcEndQuadrant, pointQuadrant, ...
startX, startY, startAng, endX, endY, endAng, direction);
if pointCovered
% If the point is on either axis and within the arc extension
% it is easy to ceck for intersect since the min allowed distance is
% along the axis, e.g normal is known.
if pointOnX
if ((a-x)<=r) || ((a+x)<=r)
res = 1;
return;
else
return;
end
end
if pointOnY
if ((b-y)<=r) || ((b+y)<=r)
res = 1;
return;
else
return;
end
end
% point is inside ellipse and not on axis start searching.
[found, currVal] = search4Normal(a, b, x, y, pointAlpha, ...
pointQuadrant, r, direction);
if found
% we have found the normal to the point, calculate distance
ax = a*cos(currVal);
ay = b*sin(currVal);
d2 = (x-ax)^2+(y-ay)^2;
if d2 <= r2
res = 1;
else
res = 0;
end
return;
else
% We did not find normal, this is an error that should be catched
res = -1; % not found normal for point outside arc.
return;
end
end % if pointCovered
end % else for: if pointInside <= 0
%end % matlab syntax
endfunction % octave syntax
% subfunctions below
function [found, currVal] = search4Normal(a, b, x, y, pointAlpha, ...
pointQuadrant, r, direction)
MAX_COUNT = 10; % Max number of iterations for solver
TOL = 0.01; % Tolerance in % for solver
maxtol = r*TOL/100; %used as tolerance for estimate if point is on normal
% we start with first estimate being the transformed
% point angle to ellipse parameterised angle
% If point is in first or third quadrant use this as max Limit for solution.
% If point is in second or fourth quadrant the angle will be min Limit
startVal = convertPhysicalToParam(pointAlpha, a, b);
[minLim, maxLim] = getSearchLimits(startVal, pointQuadrant, direction);
minVal = minLim;
maxVal = maxLim;
% make initial evaluation of equation to solve for
[solVal, soldt] = solve4normal(x, y, a, b, startVal, direction);
found = 0;
count = 0;
currVal = startVal;
while ~found && (count < MAX_COUNT)
% increment loop counter
count = count+1;
% check latest solution
if abs(solVal) <= maxtol
found = 1;
break;
end
% remember on which side the point is.
if 1 == direction
% normal is directed outward
if solVal > 0
pointIsCCW = 1; % left
else
pointIsCCW = 0; % right
end
elseif 0 == direction
% normal is directed inwards
if solVal < 0
pointIsCCW = 1; % right
else
pointIsCCW = 0; % left
end
end
% make a step in given direction.
% Steplength is calculated as
% -solVal/soldt but since the expression in question can have local
% maxima/minima between current estimate and the correct solution we
% need to make sure we take this into account.
if (solVal < 0) && (1 == direction)
% point is CW/Right of current normal, move cw, e.g negative
if soldt < 0
step = -solVal/soldt;
else
step = -(currVal-minVal)/2;
end
elseif (solVal > 0) && (0 == direction)
% point is CW/Left of current normal, move cw, e.g negative
if soldt > 0
step = -solVal/soldt;
else
step = - (currVal-minVal)/2;
end
elseif (solVal > 0) && (1 == direction)
% point is CCW/Left of current normal, move ccw, e.g positive
if soldt < 0
step = -solVal/soldt;
else
step = (maxVal-currVal)/2;
end
elseif (solVal < 0) && (0 == direction)
% point is CCW/Right of current normal, move ccw, e.g positive
if soldt > 0
step = -solVal/soldt;
else
step = (maxVal-currVal)/2;
end
elseif 0 == soldt
% no derivative, just take short step out from this point in point
% direction
if (solval < 0) && (1 == direction)
% move in negative direction
step = -(currVal-minVal)/10;
elseif (solval < 0) && (0 == direction)
% move in negative direction
step = -(currVal-minVal)/10;
elseif (solval > 0) && (1 == direction)
% move in positve direction
step = (maxVal-currVal)/10;
elseif (solval < 0) && (0 == direction)
% move in positve direction
step = (maxVal-currVal)/10;
end
end
% take the step
nextVal = currVal + step;
% check against boundarys
if nextVal > maxLim
nextVal = maxLim;
elseif nextVal < minLim
nextVal = minLim;
end
% make the evaluation with new estimatiation
[solVal, soldt] = solve4normal(x, y, a, b, nextVal, direction);
% set up the new boundarys
if (solVal < 0) && (1 == pointIsCCW) && (1 == direction)
% we have crossed the 0 by moving too much ccw, point was
% ccw before last estimation. We can reset the min/max
% values
minVal = currVal; % currVal is the previous estimate used which was too low
maxVal = nextVal; % nextVal is the last estimate which is too high
elseif (solVal > 0) && (1 == pointIsCCW) && (0 == direction)
% we have crossed the 0 by moving too much ccw, point was
% ccw before last estimation. We can reset the min/max
% values
minVal = currVal; % currVal is the previous estimate used which was too low
maxVal = nextVal; % nextVal is the last estimate which is too high
elseif (solVal > 0) && (0 == pointIsCCW) && (1 == direction)
% we have crossed the 0 by moving too much cw, point was
% cw before last estimation. We can reset the min/max
% values
minVal = nextVal; % nextVal is the last value used
maxVal = currVal; % currVal is the previous estimate
elseif (solVal < 0) && (0 == pointIsCCW) && (0 == direction)
% we have crossed the 0 by moving too much cw, point was
% cw before last estimation. We can reset the min/max
% values
minVal = nextVal; % nextVal is the last value used
maxVal = currVal; % currVal is the previous estimate
elseif (solVal < 0) && (0 == pointIsCCW) && (1 == direction)
% step in ccw direction too short set minVal
minVal = nextVal;
elseif (solVal > 0) && (0 == pointIsCCW) && (0 == direction)
% step in ccw direction too short set minVal
minVal = nextVal;
elseif (solVal > 0) && (1 == pointIsCCW) && (1 == direction)
% step in cw direction too short set maxVal
maxVal = nextVal;
elseif (solVal < 0) && (1 == pointIsCCW) && (0 == direction)
% step in cw direction too short set maxVal
maxVal = nextVal;
end
currVal = nextVal;
end % end while !found && count < MAX_COUNT
%end % matlab syntax
endfunction % octave syntax
function [minLim, maxLim] = getSearchLimits(startVal, pointQuadrant, direction)
switch pointQuadrant
case 2
if 1 == direction
minLim = 0;
maxLim = startVal;
elseif 0 == direction
minLim = startVal;
maxLim = pi/2;
end
case 4
if 1 == direction
minLim = startVal;
maxLim = pi;
elseif 0 == direction
minLim = pi/2;
maxLim = startVal;
end
case 6
if 1 == direction
minLim = pi;
maxLim = startVal;
elseif 0 == direction
minLim = startVal;
maxLim = 3*pi/2;
end
case 8
if 1 == direction
minLim = startVal;
maxLim = 2*pi;
elseif 0 == direction
minLim = 3*pi/2;
maxLim = startVal;
end
%end % matlab syntax
endswitch % octave syntax
%end % matlab syntax
endfunction % octave syntax
function pointCovered = checkArcCoverage4Point( ...
a, b, x, y, arcStartQuadrant, arcEndQuadrant, pointQuadrant, ...
startX, startY, startAng, endX, endY, endAng, direction)
% This function is almost duplicate of checkArcCoverage4PointOutside and
% could probably with some work be combined with an added argument relating
% to which side to check for.
pointCovered = 0;
% we will check for coverage evaluating the result of the expression
% (Bx-Ax)*(Cy-Ay)-(By-Ay)*(Cx-Ax)
% where the vevtor B-A gives the line direction and C is the point
% The expression will be positive if point is to the left of the line and
% negative if point is to the right of the line with 0 denominating that
% the point is on the line.
if arcStartQuadrant == arcEndQuadrant
% point must be in same quadrant and to the right of start normal and
% to the left of end normal.
if pointQuadrant ~= arcStartQuadrant
% not convered
return;
else
pointCovered = evaluateCoverage(a, b, x, y, startX, startY, ...
startAng, endX, endY, endAng, 0);
return;
end
elseif arcStartQuadrant < arcEndQuadrant
% arc does not wrap
if (arcStartQuadrant > pointQuadrant) || ...
(arcEndQuadrant < pointQuadrant)
% does not cover
return;
elseif (arcStartQuadrant < pointQuadrant) && ...
(arcEndQuadrant > pointQuadrant)
% does cover
pointCovered = 1;
return;
elseif arcStartQuadrant == pointQuadrant
val = evaluatePointDirection(a, b, x, y, startX, startY, startAng, direction);
if 1 == direction
% point is outside use outward normal
if val >= 0
% point is CCW(left) of arcStart means covered
pointCovered = 1;
return;
else
% does not cover
return;
end
elseif 0 == direction
% point is inside use inward normal
if val <= 0
% point is CCW(right) of arcStart means covered
pointCovered = 1;
return;
else
% does not cover
return;
end
end
elseif arcEndQuadrant == pointQuadrant
val = evaluatePointDirection(a, b, x, y, endX, endY, endAng, direction);
if 1 == direction
% point is outside use outward normal
if val <= 0
% point is CW(right) of arcEnd means covered
pointCovered = 1;
return;
else
% does not cover
return;
end
elseif 0 == direction
if val >= 0
% point is CW(left) of arcEnd means covered
pointCovered = 1;
return;
else
% does not cover
return;
end
end
end
else
% arc wraps around positive x axis
if (arcStartQuadrant > pointQuadrant) && ...
(arcEndQuadrant < pointQuadrant)
% does not cover
return;
elseif (arcEndQuadrant > pointQuadrant) || ...
(arcStartQuadrant < pointQuadrant)
% does cover
pointCovered = 1;
return;
elseif arcStartQuadrant == pointQuadrant
val = evaluatePointDirection(a, b, x, y, startX, startY, startAng, direction);
if 1 == direction
% point is outside use outward normal
if val >= 0
% point is CCW(left) of arcStart means covered
pointCovered = 1;
return;
else
% does not cover
return;
end
elseif 0 == direction
% point is inside use inward normal
if val <= 0
% point is CCW(right) of arcStart means covered
pointCovered = 1;
return;
else
% does not cover
return;
end
end
elseif arcEndQuadrant == pointQuadrant
val = evaluatePointDirection(a, b, x, y, endX, endY, endAng, direction);
if 1 == direction
% point is outside use outward normal
if val <= 0
% point is CW(right) of arcEnd means covered
pointCovered = 1;
return;
else
% does not cover
return;
end
elseif 0 == direction
% point is inside use inward normal
if val >= 0
% point is CW(left) of arcEnd means covered
pointCovered = 1;
return;
else
% does not cover
return;
end
end
end
end
%end % matlab syntax
endfunction % octave syntax
function pointCovered = evaluateCoverage(a, b, x, y, startX, startY, ...
startAng, endX, endY, endAng, direction)
pointCovered = 0;
% check point placement versus start/end normals.
startVal = evaluatePointDirection(a, b, x, y, startX, startY, startAng, direction);
endVal = evaluatePointDirection(a, b, x, y, endX, endY, endAng, direction);
if (startVal>=0) && (endVal<=0)
pointCovered = 1;
return;
end
%end % matlab syntax
endfunction % octave syntax
function Val = evaluatePointDirection(a, b, x, y, ax, ay, Ang, direction)
% check point placement versus normal in given direction.
if 1 == direction
% we use outside normal
bx = ax + b*cos(Ang);
by = ay + a*sin(Ang);
elseif 0 == direction
% we use inside normal
bx = ax - b*cos(Ang);
by = ay - a*sin(Ang);
end
Val = (bx-ax)*(y-ay)-(by-ay)*(x-ax);
%end % matlab syntax
endfunction % octave syntax
function res = convertPhysicalToParam(alpha, a, b)
if alpha <= pi/2
% first quadrant
res = atan(a*tan(alpha)/b);
elseif alpha <= pi
% second quadrant
res = pi+atan(a*tan(alpha)/b);
elseif alpha <= 3*pi/2
% third quadrant
res = pi+atan(a*tan(alpha)/b);
else
% fourth quadrant
res = 2*pi+atan(a*tan(alpha)/b);
end
%end % matlab syntax
endfunction % octave syntax
function alpha = getPointAngle(x,y)
% returns the physical angle calculated from positive x axis in ccw direction
q = checkQuadrant(x,y);
switch q
case 0
alpha = 0; % actually not defined but....
case 1
alpha = 0;
case 2
alpha = atan(y/x);
case 3
alpha = pi/2;
case 4
alpha = pi+atan(y/x);
case 5
alpha = pi;
case 6
alpha = pi+atan(y/x);
case 7
alpha = 3*pi/2;
case 8
alpha = 2*pi+atan(y/x);
end
%endswitch % octave syntax
end
%endfunction % octave syntax
function quadrant = checkQuadrant(x,y)
% 1 Check which quadrant point is in or if it is located
% on one of the axes
if (0==x) && (0==y)
quadrant = 0; % point has same center coordinates as arc.
elseif (x>0) && (0==y)
quadrant = 1; % point is located on positive x-axis
elseif (x>0) && (y>0)
quadrant = 2; % point is located in first quadrant
elseif (0==x) && (y>0)
quadrant = 3; % point is located on positive y-axis
elseif (x<0) && (y>0)
quadrant = 4; % point is located in second quadrant
elseif (x<0) && (0==y)
quadrant = 5; % point is located on negative x-axis
elseif (x<0) && (y<0)
quadrant = 6; % point is located in third quadrant
elseif (0==x) && (y<0)
quadrant = 7; % point is located on negative y-axis
elseif (x>0) && (y<0)
quadrant = 8; % point is located in fourth quadrant
end
%end % matlab syntax
endfunction % octave syntax
function [res, ddt] = solve4normal(x, y, a, b, t, direction)
% this function returns the value of point x, y location
% with respect to the ellips (a,b) normal at parameterised angle t
% The normal is directed:
% outwards : direction = 1
% inwards : direction = 0;
ax=a*cos(t);
ay=b*sin(t);
if 1 == direction
bx=ax+b*cos(t);
by=ay+a*sin(t);
elseif 0 == direction
bx=ax-b*cos(t);
by=ay-a*sin(t);
end
res = (bx-ax)*(y-ay)-(by-ay)*(x-ax);
p1 = b*sin(t);
p2 = a*cos(t);
if 1 == direction
ddt = -p1*(y-p1)-p2*(x-p2)-a^2*sin(t)^2-b^2*cos(t)^2;
elseif 0 == direction
ddt = p1*(y-p1)+p2*(x-p2)+a^2*sin(t)^2+b^2*cos(t)^2;
end
%end % matlab syntax
endfunction % octave syntax
% USE_CPUTIME = 0;
check_PointOutsideNotIntersect = 1;
check_PointOutsideIntersect = 1;
check_PointInsideNotIntersect = 1;
check_PointInsideIntersect = 1;
Width = 1500;
Height = 200;
Thickness = 100;
StartAng = 45;
Delta = 180;
R = 20;
[arc, pOutNis, pOutIs, pInNis, pInIs] = PrepareArcData(Width, Height, Thickness, StartAng, Delta, R);
% check the points outside the arc
% profile on
% profile off
if check_PointOutsideNotIntersect
for a = 1:length(pOutNis)
% if 1 == USE_CPUTIME
% t=cputime;
% else
% tic
% end
% profile resume
val = isPointOnArc(arc, pOutNis(a), R);
% profile off
% if 1 == USE_CPUTIME
% calctime = cputime-t;
% else
% calctime = toc;
% end
hitOutNis(a) = val; %#ok<AGROW>
% hitOutNisTime(a) = calctime; %#ok<AGROW>
if 0 == hitOutNis(a)
% Correct this point should not intersect arc
plot(pOutNis(a).X, pOutNis(a).Y, 'g.') % add green dot to cyan line
elseif 1 == hitOutNis(a)
% Not correct this point should not intersect arc
plot(pOutNis(a).X, pOutNis(a).Y, 'r+') % add red + to cyan line
else
% A negative number indicates some type of error in isPointOnArc
plot(pOutNis(a).X, pOutNis(a).Y, 'ko') % add black o to cyan line
end
end
end
if check_PointOutsideIntersect
for a = 1:length(pOutIs)
% if 1 == USE_CPUTIME
% t=cputime;
% else
% tic
% end
val = isPointOnArc(arc, pOutIs(a), R);
% if 1 == USE_CPUTIME
% calctime = cputime-t;
% else
% calctime = toc;
% end
hitOutIs(a) = val; %#ok<AGROW>
% hitOutIsTime(a) = calctime; %#ok<AGROW>
if 0 == hitOutIs(a)
% Not Correct this point should intersect arc
plot(pOutIs(a).X, pOutIs(a).Y, 'r+') % add red + to yellow line
elseif 1 == hitOutIs(a)
% Correct this point should intersect arc
plot(pOutIs(a).X, pOutIs(a).Y, 'g.') % add green dot to yellow line
else
% A negative number indicates some type of error in isPointOnArc
plot(pOutIs(a).X, pOutIs(a).Y, 'ko') % add black o to yellow line
end
end
end
if check_PointInsideNotIntersect
for a = 1:length(pInNis)
% if 1 == USE_CPUTIME
% t=cputime;
% else
% tic
% end
val = isPointOnArc(arc, pInNis(a), R);
% if 1 == USE_CPUTIME
% calctime = cputime-t;
% else
% calctime = toc;
% end
hitInNis(a) = val; %#ok<AGROW>
% hitInNisTime(a) = calctime; %#ok<AGROW>
if 0 == hitInNis(a)
% Correct this point should not intersect arc
plot(pInNis(a).X, pInNis(a).Y, 'g.') % add green dot to cyan line
elseif 1 == hitInNis(a)
% Not correct this point should intersect arc
plot(pInNis(a).X, pInNis(a).Y, 'r+') % add red + to cyan line
else
% A negative number indicates some type of error in isPointOnArc
plot(pInNis(a).X, pInNis(a).Y, 'ko') % add black o to cyan line
end
end
end
if check_PointInsideIntersect
for a = 1:length(pInIs)
% if 1 == USE_CPUTIME
% t=cputime;
% else
% tic
% end
val = isPointOnArc(arc, pInIs(a), R);
% if 1 == USE_CPUTIME
% calctime = cputime-t;
% else
% calctime = toc;
% end
hitInIs(a) = val; %#ok<AGROW>
% hitInIsTime(a) = calctime; %#ok<AGROW>
if 0 == hitInIs(a)
% Not correct this point should intersect arc
plot(pInIs(a).X, pInIs(a).Y, 'r+') % add red + to yellow line
elseif 1 == hitInIs(a)
% Correct this point should intersect arc
plot(pInIs(a).X, pInIs(a).Y, 'g.') % add green dot to yellow line
else
% A negative number indicates some type of error in isPointOnArc
plot(pInIs(a).X, pInIs(a).Y, 'ko') % add black o to yellow line
end
end
end
% profile viewer
function [arc, pOutNis, pOutIs, pInNis, pInIs] = PrepareArcData(a, b, T, startAng, Delta, R)
% [arc, pOutNis, pOutIs, pInNis, pInIs] = PrepareArcData(a, b, T, startAng, Delta, R)
%
% Prepares data for a given arc in terms of:
% a : half width along x axis
% b : half height along y axis
% T : arc thickness
% startAng : start angle in degrees for the arc with 0 along positive x-axis
% and ccw as positive direction. Note the angle represents
% the parameterised angle for ellipse equation:
% x = a*cos(t), y=b*sin(t)
% Delta : arc extension in degrees
% R : value used to create the output vectors for points that is intersecting
% and not intersecting.
%
% Outputs:
% arc : Prepared struct suitable for input to isPointOnArc
% pOutNis : array of point structs where point is outside of arc not intersecting
% pOutIs : array of point structs where point is outside of arc and intersecting
% pInNis : array of point structs where point is inside of arc not intersecting
% pInIs : array of point structs where point is inside of arc and intersecting
% create arc struct
arc.X = 0;
arc.Y = 0;
arc.Width = a;
arc.Height = b;
arc.Thickness = T;
arc.StartAng = startAng;
arc.Delta = Delta;
% help vars
rIn = T/2+R*0.9; % make the inside point 10% from R
rOut = T/2+R*1.1; % make the outside point 10% from R
% converto to radians
startAlpha = pi*startAng/180;
endAlpha = startAlpha + pi*Delta/180;
% polar vector
alpha = (startAlpha:pi/1000:endAlpha)';
% centerline
arcCenter.x = a*cos(alpha);
arcCenter.y = b*sin(alpha);
% normal to arc facing outwards
arcNormal.x = b*cos(alpha);
arcNormal.y = a*sin(alpha);
% create the arcdata and the corresponding points
% preallocate the output arrays
% TODO
% Check if arc has ends, e.g. is not a full circle/ellipse
dOut = 0;
dIn = 0;
if Delta < 360
% arc is not full circle. Edge is special case since it
% is a halfcircle.
alfa = getAlfa(arcNormal.x(1), arcNormal.y(1));
alphaOut = (alfa-pi/2:pi/1000:alfa)';
dOut = length(alphaOut);
alphaIn = (alfa-pi/2:-pi/1000:alfa-pi)';
dIn = length(alphaIn);
startHalfOut.x = arcCenter.x(1)+T/2*cos(alphaOut);
startHalfOut.y = arcCenter.y(1)+T/2*cos(alphaOut);
arcOuter.x = arcCenter.x(1)+T/2*cos(alphaOut);
arcOuter.y = arcCenter.y(1)+T/2*sin(alphaOut);
pointOuterIn.x = arcCenter.x(1)+rIn*cos(alphaOut);
pointOuterIn.y = arcCenter.y(1)+rIn*sin(alphaOut);
pointOuterOut.x = arcCenter.x(1)+rOut*cos(alphaOut);
pointOuterOut.y = arcCenter.y(1)+rOut*sin(alphaOut);
for a = 1:dOut
pOutIs(a).X = pointOuterIn.x(a); %#ok<AGROW>
pOutIs(a).Y = pointOuterIn.y(a); %#ok<AGROW>
pOutNis(a).X = pointOuterOut.x(a); %#ok<AGROW>
pOutNis(a).Y = pointOuterOut.y(a); %#ok<AGROW>
end
startHalfIn.x = arcCenter.x(1)+T/2*cos(alphaIn);
startHalfIn.y = arcCenter.y(1)+T/2*sin(alphaIn);
arcInner.x = arcCenter.x(1)+T/2*cos(alphaIn);
arcInner.y = arcCenter.y(1)+T/2*sin(alphaIn);
pointInnerIn.x = arcCenter.x(1)+rIn*cos(alphaIn);
pointInnerIn.y = arcCenter.y(1)+rIn*sin(alphaIn);
pointInnerOut.x = arcCenter.x(1)+rOut*cos(alphaIn);
pointInnerOut.y = arcCenter.y(1)+rOut*sin(alphaIn);
for a = 1:dIn
pInIs(a).X = pointInnerIn.x(a); %#ok<AGROW>
pInIs(a).Y = pointInnerIn.y(a); %#ok<AGROW>
pInNis(a).X = pointInnerOut.x(a); %#ok<AGROW>
pInNis(a).Y = pointInnerOut.y(a); %#ok<AGROW>
end
end
for a = 1:length(alpha)
alfa = getAlfa(arcNormal.x(a), arcNormal.y(a));
% outside arc points
arcOuter.x(a+dOut) = arcCenter.x(a)+T/2*cos(alfa);
arcOuter.y(a+dOut) = arcCenter.y(a)+T/2*sin(alfa);
% outside points intersecting with arc , use rIn as radius
pointOuterIn.x(a+dOut) = arcCenter.x(a)+rIn*cos(alfa);
pOutIs(a+dOut).X = pointOuterIn.x(a+dOut); %#ok<AGROW>
pointOuterIn.y(a+dOut) = arcCenter.y(a)+rIn*sin(alfa);
pOutIs(a+dOut).Y = pointOuterIn.y(a+dOut); %#ok<AGROW>
% outside points Not intersecting with arc , use rOut as radius
pointOuterOut.x(a+dOut) = arcCenter.x(a)+rOut*cos(alfa);
pOutNis(a+dOut).X = pointOuterOut.x(a+dOut); %#ok<AGROW>
pointOuterOut.y(a+dOut) = arcCenter.y(a)+rOut*sin(alfa);
pOutNis(a+dOut).Y = pointOuterOut.y(a+dOut); %#ok<AGROW>
% inside arc points
arcInner.x(a+dIn) = arcCenter.x(a)-T/2*cos(alfa);
arcInner.y(a+dIn) = arcCenter.y(a)-T/2*sin(alfa);
% inside points intersecting with arc , use rIn as radius
pointInnerIn.x(a+dIn) = arcCenter.x(a)-rIn*cos(alfa);
pInIs(a+dIn).X = pointInnerIn.x(a+dIn); %#ok<AGROW>
pointInnerIn.y(a+dIn) = arcCenter.y(a)-rIn*sin(alfa);
pInIs(a+dIn).Y = pointInnerIn.y(a+dIn); %#ok<AGROW>
% inside points Not intersecting with arc , use rOut as radius
pointInnerOut.x(a+dIn) = arcCenter.x(a)-rOut*cos(alfa);
pInNis(a+dIn).X = pointInnerOut.x(a+dIn); %#ok<AGROW>
pointInnerOut.y(a+dIn) = arcCenter.y(a)-rOut*sin(alfa);
pInNis(a+dIn).Y = pointInnerOut.y(a+dIn); %#ok<AGROW>
end % for a = 1:length(alpha)
dOut = length(arcOuter.x);
dIn = length(arcInner.x);
if Delta < 360
% arc is not full circle. Edge is special case since it
% is a halfcircle.
alfa = getAlfa(arcNormal.x(end), arcNormal.y(end));
alphaOut = (alfa:pi/1000:alfa+pi/2)';
alphaIn = (alfa+pi:-pi/1000:alfa+pi/2)';
endHalfOut.x = arcCenter.x(end)+T/2*cos(alphaOut);
endHalfOut.y = arcCenter.y(end)+T/2*sin(alphaOut);
endHalfIn.x = arcCenter.x(end)+T/2*cos(alphaIn);
endHalfIn.y = arcCenter.y(end)+T/2*sin(alphaIn);
for a = 1: length(alphaOut)
arcOuter.x(a+dOut) = arcCenter.x(end)+T/2*cos(alphaOut(a));
arcOuter.y(a+dOut) = arcCenter.y(end)+T/2*sin(alphaOut(a));
pointOuterIn.x(a+dOut) = arcCenter.x(end)+rIn*cos(alphaOut(a));
pointOuterIn.y(a+dOut) = arcCenter.y(end)+rIn*sin(alphaOut(a));
pointOuterOut.x(a+dOut) = arcCenter.x(end)+rOut*cos(alphaOut(a));
pointOuterOut.y(a+dOut) = arcCenter.y(end)+rOut*sin(alphaOut(a));
pOutIs(a+dOut).X = pointOuterIn.x(a+dOut); %#ok<AGROW>
pOutIs(a+dOut).Y = pointOuterIn.y(a+dOut); %#ok<AGROW>
pOutNis(a+dOut).X = pointOuterOut.x(a+dOut); %#ok<AGROW>
pOutNis(a+dOut).Y = pointOuterOut.y(a+dOut); %#ok<AGROW>
end
for a = 1: length(alphaIn)
arcInner.x(a+dIn) = arcCenter.x(end)+T/2*cos(alphaIn(a));
arcInner.y(a+dIn) = arcCenter.y(end)+T/2*sin(alphaIn(a));
pointInnerIn.x(a+dIn) = arcCenter.x(end)+rIn*cos(alphaIn(a));
pointInnerIn.y(a+dIn) = arcCenter.y(end)+rIn*sin(alphaIn(a));
pointInnerOut.x(a+dIn) = arcCenter.x(end)+rOut*cos(alphaIn(a));
pointInnerOut.y(a+dIn) = arcCenter.y(end)+rOut*sin(alphaIn(a));
pInIs(a+dIn).X = pointInnerIn.x(a+dIn); %#ok<AGROW>
pInIs(a+dIn).Y = pointInnerIn.y(a+dIn); %#ok<AGROW>
pInNis(a+dIn).X = pointInnerOut.x(a+dIn); %#ok<AGROW>
pInNis(a+dIn).Y = pointInnerOut.y(a+dIn); %#ok<AGROW>
end
end
figure(1)
clf
plot(arcCenter.x, arcCenter.y, 'b')
hold on
plot(arcOuter.x, arcOuter.y, 'm')
plot(arcInner.x, arcInner.y, 'm')
plot(pointOuterIn.x, pointOuterIn.y, 'y')
plot(pointInnerIn.x, pointInnerIn.y, 'y')
plot(pointOuterOut.x, pointOuterOut.y, 'c')
plot(pointInnerOut.x, pointInnerOut.y, 'c')
end % matlab syntax
%endfunction % octave syntax
% subfunction
function alfa = getAlfa(x, y)
if (x >= 0) && (y >= 0) %#ok<ALIGN>
% first quadrant
alfa = atan(y/x);
elseif (x < 0) && (y >= 0)
% second quadrant
alfa = pi+atan(y/x);
elseif (x < 0) && (y < 0)
% third quadrant
alfa = -pi+atan(y/x);
elseif (x >= 0) && (y < 0)
% fourth quadrant
alfa = pi+atan(y/x);
end
%end % matlab syntax
endfunction % octave syntax
_______________________________________________ geda-user mailing list geda-user@xxxxxxxxxxxxxx http://www.seul.org/cgi-bin/mailman/listinfo/geda-user