File: SVPolygon2dImpl.mesa
Last edited by Bier on June 1, 1984 4:51:36 pm PDT
Author: Eric Bier on February 18, 1987 6:21:46 pm PST
Contents: Routines for creating, manipulating, and playing with polygons. Part of the SolidViews package.
DIRECTORY
RealFns, SV2d, SVAngle, SVLines2d, SVPolygon2d, SVVector3d;
SVPolygon2dImpl: CEDAR PROGRAM
IMPORTS RealFns, SVAngle, SVLines2d, SVVector3d
EXPORTS SVPolygon2d =
BEGIN
Point2d: TYPE = SV2d.Point2d;
Path: TYPE = REF PathObj;
PathObj: TYPE = SV2d.PathObj;
Polygon: TYPE = REF PolygonObj;
PolygonObj: TYPE = SV2d.PolygonObj;
TrigLine: TYPE = REF TrigLineObj;
TrigLineObj: TYPE = SV2d.TrigLineObj;
TrigLineSeg: TYPE = SV2d.TrigLineSeg;
TrigPolygon: TYPE = REF TrigPolygonObj;
TrigPolygonObj: TYPE = SV2d.TrigPolygonObj;
Vector3d: TYPE = SVVector3d.Vector3d;
GLOBAL VARIABLES
unitSquare: Polygon;
root3, root3Over3, twoRoot3Over3: REAL;-- constants computed at load time in Init (see below).
globalVLine, globalCLine: TrigLine;
POLYGONS
CreatePoly:
PUBLIC
PROC [len:
NAT]
RETURNS [poly: Polygon] = {
poly ← NEW[PolygonObj[len]];
poly.len ← 0;
};
CopyPoly:
PUBLIC
PROC [poly: Polygon]
RETURNS [copy: Polygon] = {
copy ← CreatePoly[poly.maxVerts];
FOR i:
NAT
IN [0..poly.len)
DO
copy[i] ← poly[i];
ENDLOOP;
copy.len ← poly.len;
};
CircumHexagon:
PUBLIC
PROC [r:
REAL]
RETURNS [hex: Polygon] = {
Makes a hexagon with two of its sides parallel to the x axis.
hex ← CreatePoly[6];
hex ← AddPolyPoint[hex, [r*root3Over3, r]];-- front right
hex ← AddPolyPoint[hex, [r*twoRoot3Over3, 0]];-- right
hex ← AddPolyPoint[hex, [r*root3Over3, -r]];-- back right
hex ← AddPolyPoint[hex, [-r*root3Over3, -r]];-- back left
hex ← AddPolyPoint[hex, [-r*twoRoot3Over3, 0]];-- left
hex ← AddPolyPoint[hex, [-r*root3Over3, r]];-- front left
};
ClearPoly:
PUBLIC
PROC [poly: Polygon] = {
poly.len ← 0;
};
AddPolyPoint:
PUBLIC
PROC [poly: Polygon, point: Point2d]
RETURNS [polyPlusPoint: Polygon] = {
IF poly.len = poly.maxVerts
THEN {
polyPlusPoint ← CreatePoly[poly.maxVerts+1];
polyPlusPoint ← PartPolyGetsPartPoly[poly, 0, polyPlusPoint, 0, poly.maxVerts];
polyPlusPoint[poly.maxVerts] ← point;
polyPlusPoint.len ← poly.maxVerts+1;
}
ELSE {
poly[poly.len] ← point;
poly.len ← poly.len + 1;
polyPlusPoint ← poly;
};
}; -- end of AddPolyPoint
PutPolyPoint:
PUBLIC
PROC [poly: Polygon, index:
NAT, point: Point2d]
RETURNS [newPoly: Polygon] = {
IF index+1 > poly.maxVerts
THEN {
newPoly ← CreatePoly[index+1];
newPoly ← PartPolyGetsPartPoly[poly, 0, newPoly, 0, poly.maxVerts];
newPoly[index] ← point;
newPoly.len ← index+1;
}
ELSE {
IF index+1 > poly.len THEN poly.len ← index+1;
poly[index] ← point;
newPoly ← poly};
};
IsClockwisePoly:
PUBLIC
PROC [poly: Polygon]
RETURNS [
BOOL] = {
RETURN [SignedArea[poly]>0];
};
InvertPoly:
PUBLIC
PROC [poly: Polygon]
RETURNS [ylop: Polygon] = {
ylop ← NEW[PolygonObj[poly.len]];
FOR i:
NAT
IN[0..poly.len)
DO
ylop[poly.len - i -1] ← poly[i];
ENDLOOP;
};
InvertPolyInPlace:
PUBLIC
PROC [poly: Polygon] = {
tempPoint: Point2d;
FOR i:
NAT
IN[0..poly.len/2)
DO
tempPoint ← poly[poly.len - i - 1];
poly[poly.len - i - 1] ← poly[i];
poly[i] ← tempPoint;
ENDLOOP;
};
PartPolyGetsPartPath:
PUBLIC
PROC [fromPath: Path, fromStart:
NAT, toPoly: Polygon, toStart:
NAT, duration:
NAT]
RETURNS [newPoly: Polygon] = {
toPoly in range [toStart..toStart+duration-1] gets fromPath in range [fromStart..fromStart+duration-1];
j: NAT ← fromStart;
IF toStart+duration > toPoly.maxVerts
THEN {
newPoly ← CreatePoly[toStart+duration];
newPoly ← PartPolyGetsPartPoly[toPoly, 0, newPoly, 0, toPoly.len];
newPoly.len ← toStart+duration}
ELSE {
newPoly ← toPoly;
IF toStart+duration > newPoly.len THEN newPoly.len ← toStart+duration};
IF fromStart+duration > fromPath.len THEN ERROR; -- don't use bad data
FOR i:
NAT
IN[toStart..toStart+duration-1]
DO
newPoly[i] ← fromPath[j];
j ← j + 1;
ENDLOOP;
};
PartPolyGetsPartPoly:
PUBLIC
PROC [fromPoly: Polygon, fromStart:
NAT, toPoly: Polygon, toStart:
NAT, duration:
NAT]
RETURNS [newPoly: Polygon] = {
toPoly in range [toStart..toStart+duration-1] gets fromPoly in range [fromStart..fromStart+duration-1];
j: NAT ← fromStart;
IF toStart+duration > toPoly.maxVerts
THEN {
newPoly ← CreatePoly[toStart+duration];
newPoly ← PartPolyGetsPartPoly[toPoly, 0, newPoly, 0, toPoly.len];
newPoly.len ← toStart+duration}
ELSE {
newPoly ← toPoly;
IF toStart+duration > newPoly.len THEN newPoly.len ← toStart+duration};
IF fromStart+duration > fromPoly.len THEN ERROR; -- don't use bad data
FOR i:
NAT
IN[toStart..toStart+duration-1]
DO
newPoly[i] ← fromPoly[j];
j ← j + 1;
ENDLOOP;
};
PathToPolygon:
PUBLIC
PROC [path: Path]
RETURNS [poly: Polygon] = {
poly ← NEW[PolygonObj[path.len]];
FOR i:
NAT
IN[0..path.len)
DO
poly[i] ← path[i];
ENDLOOP;
poly.len ← path.len;
}; -- end of PathToPolygon
TRIG POLYGONS (a richer [and slower] structure)
PolygonToTrigPolygon:
PUBLIC
PROC [poly: Polygon]
RETURNS [trigPoly: TrigPolygon] = {
lastPoint, thisPoint: Point2d;
newSeg: TrigLineSeg;
trigPoly ← NEW[TrigPolygonObj[poly.len]];
lastPoint ← poly[0];
FOR i:
NAT
IN[1..poly.len)
DO
thisPoint ← poly[i];
newSeg ← SVLines2d.CreateTrigLineSeg[lastPoint, thisPoint];
trigPoly[i-1] ← newSeg;
lastPoint ← thisPoint;
ENDLOOP;
thisPoint ← poly[0];
newSeg ← SVLines2d.CreateTrigLineSeg[lastPoint, thisPoint];
trigPoly[poly.len-1] ← newSeg;
};
CreatePath:
PUBLIC
PROC [len:
NAT]
RETURNS [path: Path] = {
path ← NEW[PathObj[len]];
path.len ← 0;
};
CopyPath:
PUBLIC
PROC [path: Path]
RETURNS [copy: Path] = {
copy ← CreatePath[path.maxVerts];
FOR i:
NAT
IN [0..path.len)
DO
copy[i] ← path[i];
ENDLOOP;
copy.len ← path.len;
};
ClearPath:
PUBLIC
PROC [path: Path] = {
path.len ← 0;
};
AddPathPoint:
PUBLIC
PROC [path: Path, point: Point2d]
RETURNS [pathPlusPoint: Path] = {
IF path.len = path.maxVerts
THEN {
pathPlusPoint ← CreatePath[path.maxVerts+1];
pathPlusPoint ← PartPathGetsPartPath[path, 0, pathPlusPoint, 0, path.maxVerts];
pathPlusPoint[path.maxVerts] ← point;
pathPlusPoint.len ← path.maxVerts+1;
}
ELSE {
path[path.len] ← point;
path.len ← path.len + 1;
pathPlusPoint ← path;
};
}; -- end of AddPathPoint
InsertPathPoint:
PUBLIC
PROC [path: Path, point: Point2d]
RETURNS [newPath: Path] = {
IF path.len = path.maxVerts
THEN {
newPath ← CreatePath[path.maxVerts+1];
newPath[0] ← point;
newPath ← PartPathGetsPartPath[path, 0, newPath, 1, path.maxVerts];
newPath.len ← path.maxVerts+1;
}
ELSE {
ShiftUpPath[path, 0, 1];
path[0] ← point;
path.len ← path.len + 1;
newPath ← path;
};
}; -- end of InsertPathPoint
SplicePathPoint:
PUBLIC
PROC [path: Path, point: Point2d, index:
INTEGER]
RETURNS [newPath: Path]= {
IF index is -1, splice onto beginning, if index is path.len-1, splice onto end.
IF index > path.len-1 THEN ERROR AttemptToSplicePastEndOfPath;
IF index = -1 THEN {newPath ← InsertPathPoint[path, point]; RETURN};
IF index = path.len-1
THEN {
newPath ← AddPathPoint[path, point]; RETURN};
IF path.len = path.maxVerts
THEN {
newPath ← CreatePath[path.maxVerts+1];
newPath ← PartPathGetsPartPath[path, 0, newPath, 0, path.maxVerts];
ShiftUpPath[newPath, index+1, 1];
newPath[index+1] ← point;
}
ELSE {
ShiftUpPath[path, index+1, 1];
path[index+1] ← point;
newPath ← path;
};
}; -- end of SplicePathPoint
AttemptToSplicePastEndOfPath: PUBLIC ERROR = CODE;
DeletePathPoint:
PUBLIC
PROC [path: Path, index:
NAT]
RETURNS [newPath: Path] = {
IF index> path.len-1 THEN ERROR AttemptToDeleteNonExistentPoint;
IF path.len = 0 THEN ERROR PathEmpty;
IF index = path.len-1 THEN {path.len ← path.len -1; newPath ← path; RETURN};
ShiftDownPath[path, index+1, 1];
newPath ← path;
}; -- end of DeletePathPoint
AttemptToDeleteNonExistentPoint: PUBLIC ERROR = CODE;
PathEmpty: PUBLIC ERROR = CODE;
PutPathPoint:
PUBLIC
PROC [path: Path, index:
NAT, point: Point2d] = {
IF index+1 > path.maxVerts THEN ERROR;
IF index+1 > path.len THEN path.len ← index+1;
path[index] ← point;
};
ConcatPath:
PUBLIC
PROC [path1, path2: Path]
RETURNS [cat: Path] = {
cat ← CreatePath[path1.len+path2.len];
FOR i:
NAT
IN[0..path1.len)
DO
cat[i] ← path1[i];
ENDLOOP;
FOR i:
NAT
IN[path1.len..path1.len+path2.len)
DO
cat[i] ← path2[i-path1.len];
ENDLOOP;
};
SubPath:
PUBLIC
PROC [path: Path, lo, hi:
NAT]
RETURNS [subpath: Path] = {
IF hi < lo THEN ERROR;
subpath ← CreatePath[hi-lo+1];
FOR i:
NAT
IN[lo..hi]
DO
subpath[i - lo] ← path[i];
ENDLOOP;
subpath.len ← hi-lo+1;
};
SubPathOfPoly:
PUBLIC
PROC [poly: Polygon, lo, hi:
NAT]
RETURNS [subpath: Path] = {
IF hi < lo THEN ERROR;
subpath ← CreatePath[hi-lo+1];
FOR i:
NAT
IN[lo..hi]
DO
subpath[i - lo] ← poly[i];
ENDLOOP;
subpath.len ← hi-lo+1;
};
ShiftUpPath:
PUBLIC
PROC [path: Path, startAt:
NAT, by:
NAT] = {
startAt is the index of the lowest element which should be shifted up
IF path.len + by > path.maxVerts THEN ERROR;
FOR i:
NAT
DECREASING
IN [startAt..path.len)
DO
path[i+by] ← path[i];
ENDLOOP;
path.len ← path.len + by;
}; -- end of ShiftUpPath
ShiftDownPath:
PUBLIC
PROC [path: Path, startAt:
NAT, by:
NAT] = {
startAt the index of the highest element which should be shifted down. Must be greater than zero.
IF startAt - by < 0 THEN ERROR ShiftingDataOffLeftEnd;
IF startAt > path.len-1 THEN ERROR ShiftingDownNonExistentElements;
FOR i:
NAT
IN [startAt..path.len)
DO
path[i-by] ← path[i];
ENDLOOP;
path.len ← path.len - by;
}; -- end of ShiftDownPath
ShiftingDataOffLeftEnd: PUBLIC ERROR = CODE;
ShiftingDownNonExistentElements: PUBLIC ERROR = CODE;
PolygonToPath:
PUBLIC
PROC [poly: Polygon]
RETURNS [path: Path] = {
path ← NEW[PathObj[poly.len]];
FOR i:
NAT
IN[0..poly.len)
DO
path[i] ← poly[i];
ENDLOOP;
path.len ← poly.len;
}; -- end of PolygonToPath
PartPathGetsPartPath:
PUBLIC
PROC [fromPath: Path, fromStart:
NAT, toPath: Path, toStart:
NAT, duration:
NAT]
RETURNS [newPath: Path] = {
toPath in range [toStart..toStart+duration-1] gets fromPath in range [fromStart..fromStart+duration-1];
j: NAT ← fromStart;
IF toStart+duration > toPath.maxVerts
THEN {
newPath ← CreatePath[toStart+duration];
newPath ← PartPathGetsPartPath[toPath, 0, newPath, 0, toPath.len];
newPath.len ← toStart+duration}
ELSE {
newPath ← toPath;
IF toStart+duration > newPath.len THEN newPath.len ← toStart+duration};
IF fromStart+duration > fromPath.len THEN ERROR; -- don't use bad data
FOR i:
NAT
IN[toStart..toStart+duration-1]
DO
newPath[i] ← fromPath[j];
j ← j + 1;
ENDLOOP;
};
PointPolyClass: TYPE = SVPolygon2d.PointPolyClass;-- {in, on, out};
CirclePointInPoly:
PUBLIC
PROC [point: Point2d, poly: Polygon]
RETURNS [class: PointPolyClass] = {
Take point. Take a vertex V and its clockwise neighbor C from poly. Find the angle from V to C. Add up these angles for all such pairs V and C in poly. If the total is -360 degress, the point is in the poly. If the total is zero, the point is outside. Say +or- 3.6 degrees is close enough to zero and -360 +or- 3.6 is good enough for 360. Otherwise, error.
to find the angles, I could take the cross product of the vectors from point to V and to C, normalize, and the take ArcSin, but Cedar doesn't have ArcSin so for now, I will use the ArcTan machinery available from SV2d.CreateTrigLineSeg to find the angle of each vector with respect to the x-axis and then take the difference of these two angles.
v, c: Point2d;
thisTheta: REAL;
iPrime: NAT;
totalTheta: REAL ← 0;
v ← poly[0];
SVLines2d.FillTrigLineFromPoints[point, v, globalVLine];
FOR i:
NAT
IN[1..poly.len]
DO
iPrime ← IF i = poly.len THEN 0 ELSE i;
c ← poly[iPrime];
SVLines2d.FillTrigLineFromPoints[point, c, globalCLine];
thisTheta ← SVAngle.ShortestDifference[globalCLine.theta, globalVLine.theta];
IF thisTheta = 180 OR thisTheta = -180 THEN RETURN[on];
totalTheta ← totalTheta + thisTheta;
v ← c;
SVLines2d.CopyTrigLine[from: globalCLine, to: globalVLine];
ENDLOOP;
IF -3.6 < totalTheta AND totalTheta < 3.6 THEN RETURN [out];
IF -363.6 < totalTheta AND totalTheta < -356.4 THEN RETURN [in]
ELSE RETURN [in]; -- **** should be an error
};
SquarePointInPoly:
PUBLIC
PROC [point: Point2d, poly: Polygon]
RETURNS [class: PointPolyClass] = {
RETURN[in];
};
UpDown: TYPE = {up, on, down};
LeftRight: TYPE = {left, on, right};
YComparedToHorizLine:
PRIVATE
PROC [testY:
REAL, yLine:
REAL]
RETURNS [UpDown] = {
IF testY > yLine THEN RETURN[up];
IF testY = yLine THEN RETURN[on];
RETURN[down];
};
XComparedToVertLine:
PRIVATE
PROC [testX:
REAL, xLine:
REAL]
RETURNS [LeftRight] = {
IF testX > xLine THEN RETURN[right];
IF testX = xLine THEN RETURN[on];
RETURN[left];
};
LineSegHitsHorizLineLeftOfT:
PRIVATE
PROC [loPoint: Point2d, hiPoint: Point2d, t, lineY:
REAL]
RETURNS [LeftRight] = {
We wish to know if the (infinite) line passing through loPoint and hiPoint hits the horizontal line y = lineY to the left, right, or smack dab on the point (t, lineY).
Let (x0,y0) = loPoint and (x1,y1) = hiPoint then the line equation can be expressed as:
y - y0 = (y1-y0)/(x1-x0)(x - x0) unless x1 = x0.
IF x1 = x0 then we compare x0 to t and return result.
Otherwise, we solve for x where y = lineY:
lineY - y0 = [(y1-y0)/(x1-x0)](x - x0)
(lineY - y0)(x1-x0) = (y1-y0)(x - x0)
(lineY - y0)(x1-x0)/(y1-y0) = x - x0
x0 + (lineY - y0)(x1-x0)/(y1-y0) = x.
So, x > t => x0 + (lineY - y0)(x1-x0)/(y1-y0) > t.
taking 1 mult and 1 div which =approx. 3 mults.
x: REAL;
IF loPoint[1] = hiPoint[1] THEN RETURN[XComparedToVertLine[loPoint[1], t]];
x ← loPoint[1] + (lineY - loPoint[2])*(hiPoint[1]-loPoint[1])/(hiPoint[2]-loPoint[2]);
RETURN[XComparedToVertLine[x, t]];
};
IsEven:
PRIVATE
PROC [n:
NAT]
RETURNS [
BOOL] = {
RETURN [(n/2)*2 = n];
};
BoysePointInPoly:
PUBLIC
PROC [point: Point2d, poly: Polygon]
RETURNS [class: PointPolyClass] = {
Consider a horiz line and a vert line passing through p. How many edges of poly cross the horizontal line to the right of the vertical line? Call this number "hits". If hits is even then p is outside of poly. If hits is odd, the p is inside poly. If we find an edge of poly which passes through p, then p is on poly;
Call the vertical line vLine, (x = p[1]) and the horizontal line hLine (y = p[2]).
Method: There are nine starting cases. The first vert of poly may be up, on, or down of hLine. It may be left, on or right of vLine. Similarly, there are nine cases for the other endpoint. Most cases make it easy to determine whether this edge of the poly cross hLine to the right of vLine or not. The hard cases are:
1) up-left to down-right
2) down-left to up-right
For these cases, we use LineSegHitsHorizLineLeftOfT where t is p[1]
Hopefully, there will be some way of coding this without mentioning all 81 cases explicitly.
currentX, nextX: LeftRight;
currentY, nextY, resolveY: UpDown;
i: NAT;
crossings: NAT ← 0;
currentPoint: Point2d ← poly[0];
intersect: LeftRight;
currentX ← XComparedToVertLine[currentPoint[1], point[1]];
currentY ← YComparedToHorizLine[currentPoint[2], point[2]];
resolveY ← currentY;
FOR j:
NAT
IN [1..poly.len]
DO
i ← IF j = poly.len THEN 0 ELSE j;
nextX ← XComparedToVertLine[poly[i][1], point[1]];
nextY ← YComparedToHorizLine[poly[i][2], point[2]];
SELECT TRUE FROM
currentX = left
AND currentY = up =>
SELECT TRUE FROM
nextX = left AND nextY = on => {currentY ← on; resolveY ← up}; -- we stay left so no right crossings.
nextX = left => {resolveY ← currentY ← nextY}; -- we stay left so no right crossings
nextX = on AND nextY = up => {currentX ← on; resolveY ← up};
nextX = on AND nextY = on => {RETURN[on]};-- p is on a vertex of poly
nextX = on AND nextY = down => {currentX ← on; resolveY ← currentY ← down}; -- cross to the left
nextX = right AND nextY = up => {currentX ← right};-- cross vert line only
nextX = right AND nextY = on => {currentX ← right; currentY ← on; resolveY ← up}; -- we still consider ourselves up
nextX = right AND nextY = down => {-- one of the hard cases
currentX ← right; resolveY ← currentY ← down;
intersect ← LineSegHitsHorizLineLeftOfT[poly[i], currentPoint, point[1], point[2]];
IF intersect = right THEN crossings ← crossings + 1
ELSE IF intersect = on THEN RETURN[on];
};
ENDCASE => ERROR;
currentX = left
AND currentY = on =>
SELECT TRUE FROM
nextX = left AND nextY = on => {currentY ← on}; -- resolve y unchanged
nextX = left => resolveY ← currentY ← nextY;
nextX = on AND nextY = up => {currentX ← on; resolveY ← currentY ← up};
nextX = on AND nextY = on => {RETURN[on]};-- p is on a vertex of poly
nextX = on AND nextY = down => {currentX ← on; resolveY ← currentY ← down};
nextX = right AND nextY = up => {currentX ← right; resolveY ← currentY ← up}; -- cross vert line only
nextX = right AND nextY = on => {RETURN[on]}; -- we crossed thru p.
nextX = right AND nextY = down => {currentX ← right; resolveY ← currentY ← down};
ENDCASE => ERROR;
currentX = left
AND currentY = down =>
SELECT TRUE FROM
nextX = left AND nextY = on => {currentY ← on}; -- we stay left. resolveY unchanged
nextX = left => {resolveY ← currentY ← nextY}; -- we stay left. Take upDown of next
nextX = on AND nextY = up => {currentX ← on; resolveY ← currentY ← up};
nextX = on AND nextY = on => {RETURN[on]};-- p is on a vertex of poly
nextX = on AND nextY = down => {currentX ← on};
nextX = right
AND nextY = up => {-- one of the hard cases
currentX ← right; resolveY ← currentY ← up;
intersect ← LineSegHitsHorizLineLeftOfT[currentPoint, poly[i], point[1], point[2]];
IF intersect = right THEN crossings ← crossings + 1
ELSE IF intersect = on THEN RETURN[on];
};
nextX = right AND nextY = on => {currentX ← right; currentY ← on; resolveY ← down}; -- we still consider ourselves down
nextX = right AND nextY = down => {currentX ← right};
ENDCASE => ERROR;
currentX = on
AND currentY = up =>
SELECT TRUE FROM
nextX = left AND nextY = on => {currentX ← left; currentY ← on; resolveY ← up}; -- we go left.
nextX = left => {currentX ← left; resolveY ← currentY ← nextY};
nextX = on AND nextY = up => {};-- no change
nextX = on AND nextY = on => RETURN[on];-- p is on a vertex of poly
nextX = on AND nextY = down => RETURN[on]; -- we just crossed thru p
nextX = right AND nextY = up => {currentX ← right};
nextX = right AND nextY = on => {currentX ← right; currentY ← on; resolveY ← up}; -- we still consider ourselves up
nextX = right AND nextY = down => {-- an automatic crossing!
currentX ← right; resolveY ← currentY ← down;
crossings ← crossings + 1};
ENDCASE => ERROR;
currentX = on AND currentY = on => RETURN[on]; -- should never occur. Should be caught sooner.
currentX = on
AND currentY = down =>
SELECT TRUE FROM
nextX = left AND nextY = on => {currentX ← left; currentY ← on; resolveY ← down}; -- we go left.
nextX = left => {currentX ← left; resolveY ← currentY ← nextY};
nextX = on AND nextY = up => RETURN[on]; -- we just crossed thru p
nextX = on AND nextY = on => RETURN[on];-- p is on a vertex of poly
nextX = on AND nextY = down => {}; -- no change
nextX = right AND nextY = up => {-- automatic crossing!
currentX ← right; resolveY ← currentY ← up;
crossings ← crossings + 1};
nextX = right AND nextY = on => {currentX ← right; currentY ← on; resolveY ← down}; -- we still consider ourselves down
nextX = right AND nextY = down => {currentX ← right};
ENDCASE => ERROR;
currentX = right
AND currentY = up =>
SELECT TRUE FROM
nextX = left AND nextY = up => {currentX ← left};
nextX = left AND nextY = on => {currentX ← left; currentY ← on; resolveY ← up};
nextX = left
AND nextY = down => {-- one of the hard cases
currentX ← left; resolveY ← currentY ← down;
intersect ← LineSegHitsHorizLineLeftOfT[poly[i], currentPoint, point[1], point[2]];
IF intersect = right THEN crossings ← crossings + 1
ELSE IF intersect = on THEN RETURN[on];
};
nextX = on AND nextY = up => {currentX ← on};
nextX = on AND nextY = on => {RETURN[on]};-- p is on a vertex of poly
nextX = on
AND nextY = down => {-- automatic crossing!
currentX ← on; resolveY ← currentY ← down;
crossings ← crossings + 1;
};
nextX = right AND nextY = up => {}; -- no change
nextX = right AND nextY = on => {currentY ← on; resolveY ← up}; -- we still consider ourselves up
nextX = right AND nextY = down => { -- automatic crossing!
resolveY ← currentY ← down;
crossings ← crossings + 1;
};
ENDCASE => ERROR;
currentX = right
AND currentY = on =>
SELECT TRUE FROM
nextX = left
AND nextY = up => {-- this is a crossing if resolveY is down
IF resolveY = down THEN crossings ← crossings + 1;
currentX ← left; resolveY ← currentY ← up};
nextX = left AND nextY = on => {RETURN[on]};-- we pass thru p
nextX = left
AND nextY = down => {-- this is a crossing if resolveY is up
IF resolveY = up THEN crossings ← crossings + 1;
currentX ← left; resolveY ← currentY ← down};
nextX = on
AND nextY = up => {-- this is a crossing if resolveY is down
IF resolveY = down THEN crossings ← crossings + 1;
currentX ← on; resolveY ← currentY ← up};
nextX = on AND nextY = on => {RETURN[on]};-- p is on a vertex of poly
nextX = on
AND nextY = down => {-- this is a crossing if resolveY is up
IF resolveY = up THEN crossings ← crossings + 1;
currentX ← on; resolveY ← currentY ← down};
nextX = right
AND nextY = up => {-- this is a crossing if resolveY is down
IF resolveY = down THEN crossings ← crossings + 1;
resolveY ← currentY ← up};
nextX = right AND nextY = on => {}; -- no change
nextX = right
AND nextY = down => {-- this is a crossing if resolveY is up
IF resolveY = up THEN crossings ← crossings + 1;
resolveY ← currentY ← down};
ENDCASE => ERROR;
currentX = right
AND currentY = down =>
SELECT TRUE FROM
nextX = left
AND nextY = up => {-- one of the hard cases
currentX ← left; resolveY ← currentY ← up;
intersect ← LineSegHitsHorizLineLeftOfT[currentPoint, poly[i], point[1], point[2]];
IF intersect = right THEN crossings ← crossings + 1
ELSE IF intersect = on THEN RETURN[on];
};
nextX = left AND nextY = on => {currentX ← left; currentY ← on; resolveY ← down};
nextX = left AND nextY = down => {currentX ← left};
nextX = on
AND nextY = up => {-- automatic crossing!
currentX ← on; resolveY ← currentY ← up;
crossings ← crossings + 1};
nextX = on AND nextY = on => {RETURN[on]};-- p is on a vertex of poly
nextX = on AND nextY = down => {currentX ← on};
nextX = right
AND nextY = up => {-- automatic crossing!
resolveY ← currentY ← up;
crossings ← crossings + 1};
nextX = right AND nextY = on => {currentY ← on; resolveY ← down}; -- we still consider ourselves down
nextX = right AND nextY = down => {};-- no change
ENDCASE => ERROR;
ENDCASE => ERROR;
currentPoint ← poly[i];
ENDLOOP;
IF IsEven[crossings] THEN RETURN[out] ELSE RETURN[in];
};
PointToVector:
PRIVATE
PROC [point: Point2d]
RETURNS [vector: Vector3d] = {
vector[1] ← point[1];
vector[2] ← point[2];
vector[3] ← 0;
};
SignedArea:
PUBLIC
PROC [poly: Polygon]
RETURNS [area:
REAL] = {
Chose an arbitrary point P in the plane. Chose two consecutive vertices (v1, v2) on the polygon. Call the vector P to v1, D1. Call the distance P to v2, D2. The area of the triangle made by these three points is |D1|*|D2|*sin(angle between D1 and D2)/2.
This is just D2 x D1 where "x" is the cross product operator.
This will be positive if D2 is clockwise of D1.
If P is in the interior of poly for poly convex, and if we add up all of the areas we get by taking vertices pairwise in clockwise order, area will be positive. If we take them in counter-clockwise order, area will be negative. This turns out to be true whatever the position of P and even if poly is concave part of the time.
D1, D2: Vector3d;
iPlusOne: NAT;
areaVector: Vector3d;
partialArea: REAL;
area ← 0;
let P be the origin so the vector is the same as its endpoint
FOR i: NAT IN [0..poly.len-1] DO
iPlusOne ← IF i = poly.len-1 THEN 0 ELSE i + 1;
D1 ← PointToVector[poly[i]];
D2 ← PointToVector[poly[iPlusOne]];
areaVector ← SVVector3d.CrossProduct[D2, D1];-- will only have a z component
partialArea ← areaVector[3];
area ← area + partialArea;
ENDLOOP;
}; -- end of signed area
ClockwisePerimeterAroundUnitSquare:
PUBLIC
PROC [from, to: Point2d]
RETURNS [perim:
REAL] = {
Find out which sides from and to are on. Number the sides as follows: 0 = left, 1 = top, 2 = right, 3 = bottom.
unitsquare is a global polygon. unitsquare[0] = lower left; [1] = upper left; [2] = upper right; [3] = lower right;
Side: TYPE = NAT;
side: Side;
fromSide, toSide: Side;
perimPart: REAL;
SELECT TRUE FROM
from[1] = -0.5 => fromSide ← 0;
from[2] = 0.5 => fromSide ← 1;
from[1] = 0.5 => fromSide ← 2;
from[2] = -0.5 => fromSide ← 3;
ENDCASE => ERROR;
SELECT TRUE FROM
to[1] = -0.5 => toSide ← 0;
to[2] = 0.5 => toSide ← 1;
to[1] = 0.5 => toSide ← 2;
to[2] = -0.5 => toSide ← 3;
ENDCASE => ERROR;
side ← fromSide;
IF toSide = fromSide
THEN {
IF IsClockwisePtoQAlongSquareSide[from, to, toSide] THEN
perim ← DistanceClockwiseAlongSquareSide[from, to, toSide]
ELSE perim ← -DistanceClockwiseAlongSquareSide[to, from, toSide];
RETURN};
toSide # fromSide so add on the clockwise distance to the next vertex and keep going.
perim ← DistanceClockwiseAlongSquareSide[from, unitSquare[fromSide+1], fromSide];
FOR i:
NAT
IN[1..3]
DO
side ← fromSide + i;
IF side > 3 THEN side ← side -4;
IF side = toSide THEN GOTO ToFound;
perim ← perim + 1;
REPEAT
ToFound => {-- find distance from last vertex to this point
perimPart ← DistanceClockwiseAlongSquareSide[unitSquare[side], to, side];
perim ← perim + perimPart;
};
ENDLOOP;
IF perim > 2 THEN perim ← perim - 4;
}; -- end of ClockwisePerimeterAroundUnitSquare
IsClockwisePtoQAlongSquareSide:
PROC [p, q: Point2d, side:
NAT]
RETURNS [
BOOL] = {
SELECT side FROM
0 => -- left side. p to q is clockwise if q above p;
RETURN[q[2] > p[2]];
1 => -- top side. p to q is clockwise if q right of p;
RETURN[q[1] > p[1]];
2 => -- right side. p to q is clockwise if q below p;
RETURN[q[2] < p[2]];
3 => -- bottom side. p to q is clockwise if q left of p;
RETURN[q[1] < p[1]];
ENDCASE => ERROR;
}; -- end of IsClockwisePtoQAlongSquareSide
DistanceClockwiseAlongSquareSide:
PROC [p, q: Point2d, side:
NAT]
RETURNS [
REAL] = {
Assumes that q is clockwise from p;
SELECT side FROM
0 => -- left side. Distance p to q is difference in y's;
RETURN[q[2] - p[2]];
1 => -- top side. Distance p to q is difference in x's;
RETURN[q[1] - p[1]];
2 => -- right side. Distance p to q is difference in y's;
RETURN[p[2] - q[2]];
3 => -- bottom side. Distance p to q is difference in x's;
RETURN[p[1] - q[1]];
ENDCASE => ERROR;
}; -- end of DistanceClockwiseAlongSquareSide
Init:
PROC = {
unitSquare ← CreatePoly[4];
unitSquare ← AddPolyPoint[unitSquare, [-.5, -.5]];
unitSquare ← AddPolyPoint[unitSquare, [-.5, .5]];
unitSquare ← AddPolyPoint[unitSquare, [.5, .5]];
unitSquare ← AddPolyPoint[unitSquare, [.5, -.5]];
root3 ← RealFns.SqRt[3];
root3Over3 ← root3/3.0;
twoRoot3Over3 ← 2*root3Over3;
globalVLine ← SVLines2d.CreateEmptyTrigLine[];
globalCLine ← SVLines2d.CreateEmptyTrigLine[];
}; -- end of Init
Init[];
END.