-- File: SVPolygon2dImpl.mesa
-- Last edited by Bier on December 18, 1982 1:11 am
-- Author: Eric Bier on August 20, 1982 12:52 pm
-- Contents: Routines for creating, manipulating, and playing with polygons. Part of the SolidViews package.
DIRECTORY
RealFns,
SV2d,
SVAngle,
SVLines2d,
SVPolygon2d,
SVVector3d;
SVPolygon2dImpl: 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;
Vector: TYPE = SVVector3d.Vector;
-- 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: Vector] = {
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.
-- These is just D1 x D2 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: Vector;
iPlusOne: NAT;
areaVector: Vector;
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.