
DIRECTORY G2dBasic, G2dVector, G3dBasic, G3dMatrix, G3dPlane, G3dPolygon, G3dVector, Real, RealFns, RefTab, Vector2;

G3dPolygonImpl: CEDAR PROGRAM
IMPORTS G2dBasic, G2dVector, G3dMatrix, G3dPlane, G3dVector, Real, RealFns, RefTab, Vector2
EXPORTS G3dPolygon

~ BEGIN
Border:							TYPE ~ G2dBasic.Border;
NatSequence:						TYPE ~ G2dBasic.NatSequence;
Pair:								TYPE ~ G2dBasic.Pair;
PairSequence:					TYPE ~ G2dBasic.PairSequence;
Triple:							TYPE ~ G2dBasic.Triple;
TripleSequence:					TYPE ~ G2dBasic.TripleSequence;
TripleSequenceRep:				TYPE ~ G2dBasic.TripleSequenceRep;
NatPair:							TYPE ~ G3dBasic.NatPair;
SurfaceSequence:				TYPE ~ G3dBasic.SurfaceSequence;
Quad:								TYPE ~ G3dBasic.Quad;
QuadSequenceRep:				TYPE ~ G3dBasic.QuadSequenceRep;
RealSequence:					TYPE ~ G3dBasic.RealSequence;
RealSequenceRep:				TYPE ~ G3dBasic.RealSequenceRep;
Segment:							TYPE ~ G3dBasic.Segment;
Matrix:							TYPE ~ G3dMatrix.Matrix;
MajorPlane:						TYPE ~ G3dPlane.MajorPlane;
Plane:								TYPE ~ G3dPlane.Plane;
NearPolygon:						TYPE ~ G3dPolygon.NearPolygon;
NearTriangle:					TYPE ~ G3dPolygon.NearTriangle;
PolygonProc:						TYPE ~ G3dPolygon.PolygonProc;
Polygon:							TYPE ~ G3dPolygon.Polygon;
PolygonSequence:				TYPE ~ G3dPolygon.PolygonSequence;
PolygonSequenceRep:			TYPE ~ G3dPolygon.PolygonSequenceRep;
Triangle:							TYPE ~ G3dPolygon.Triangle;
TriangleSequence:				TYPE ~ G3dPolygon.TriangleSequence;
TriangleSequenceRep:			TYPE ~ G3dPolygon.TriangleSequenceRep;
Nat3Sequence:					TYPE ~ G3dPolygon.Nat3Sequence;
Nat3SequenceRep:				TYPE ~ G3dPolygon.Nat3SequenceRep;
SetTriangle: PUBLIC PROC [triangle: Triangle, accs, lines: BOOL Ź FALSE] ~ {
triangle.plane Ź
G3dPlane.Unit[G3dPlane.FromThreePoints[triangle.p0, triangle.p1, triangle.p2]];
triangle.majorPlane Ź G3dPlane.GetMajorPlane[triangle.plane];
IF accs THEN {
triangle.a0 Ź G3dVector.NearnessAccelerator[triangle.p0, triangle.p1];
triangle.a1 Ź G3dVector.NearnessAccelerator[triangle.p1, triangle.p2];
triangle.a2 Ź G3dVector.NearnessAccelerator[triangle.p2, triangle.p0];
};
IF lines THEN {
q0: Pair Ź G3dPlane.ProjectPointToMajorPlane[triangle.p0, triangle.majorPlane];
q1: Pair Ź G3dPlane.ProjectPointToMajorPlane[triangle.p1, triangle.majorPlane];
q2: Pair Ź G3dPlane.ProjectPointToMajorPlane[triangle.p2, triangle.majorPlane];
IF (SELECT triangle.majorPlane FROM
xy => triangle.plane.z, xz => triangle.plane.y, ENDCASE => triangle.plane.x) > 0.0
THEN {
triangle.l0 Ź G2dVector.Line[q0, q1];
triangle.l1 Ź G2dVector.Line[q1, q2];
triangle.l2 Ź G2dVector.Line[q2, q0];
}
ELSE { -- invert from normal sense
triangle.l0 Ź G2dVector.Line[q1, q0];
triangle.l1 Ź G2dVector.Line[q2, q1];
triangle.l2 Ź G2dVector.Line[q0, q2];
};
};
};

SetPolygon: PUBLIC PROC [
polygon: Polygon,
vertices: TripleSequence Ź NIL,
normal, center, area, accs, lines: BOOL Ź FALSE]
~ {
IF polygon.points = NIL AND vertices = NIL THEN RETURN;
IF polygon.points = NIL THEN {
IF polygon.points = NIL OR polygon.points.maxLength < polygon.indices.length
THEN polygon.points Ź NEW[TripleSequenceRep[polygon.indices.length]];
polygon.points.length Ź polygon.indices.length;
FOR n: NAT IN [0..polygon.points.length) DO
polygon.points[n] Ź vertices[polygon.indices[n]];
ENDLOOP;
};
IF polygon.plane = []
THEN polygon.plane Ź G3dPlane.Unit[G3dPlane.FromPolygon[polygon.points]];
IF polygon.majorPlane = unknown
THEN polygon.majorPlane Ź G3dPlane.GetMajorPlane[polygon.plane];
IF normal THEN polygon.normal Ź [polygon.plane.x, polygon.plane.y, polygon.plane.z];
IF accs THEN {
polygon.accs Ź NEW[QuadSequenceRep[polygon.points.length]];
polygon.accs.length Ź polygon.points.length;
};
IF lines THEN {
polygon.lines Ź NEW[TripleSequenceRep[polygon.points.length]];
polygon.lines.length Ź polygon.points.length;
};
IF accs OR lines THEN FOR n: NAT IN [0..polygon.points.length) DO
nn: NAT Ź (n+1) MOD polygon.points.length;
p1: Triple Ź polygon.points[n];
p2: Triple Ź polygon.points[nn];
IF lines THEN polygon.lines[n] Ź SELECT polygon.majorPlane FROM
xy			=> G2dVector.Line[[p1.x, p1.y],	[p2.x, p2.y]],
xz			=> G2dVector.Line[[p1.x, p1.z],	[p2.x, p2.z]],
ENDCASE	=> G2dVector.Line[[p1.y, p1.z],	[p2.y, p2.z]];
IF accs THEN polygon.accs[n] Ź G3dVector.NearnessAccelerator[p1, p2];
ENDLOOP;
};

SetPolygonNeighbors: PUBLIC PROC [polygons: PolygonSequence] ~ {
refNatPair: REF NatPair Ź NEW[NatPair];
neighbors: RefTab.Ref Ź RefTab.Create[equal: Equal, hash: Hash];
Store: PROC [poly: Polygon] ~ {
i: NatSequence Ź poly.indices;
FOR n: NAT IN [0..i.length) DO
[] Ź RefTab.Store[neighbors, NEW[NatPair Ź [i[n], i[(n+1) MOD i.length]]], poly];
ENDLOOP;
};
SetNeighbor: PROC [poly: Polygon] ~ {
i: NatSequence Ź poly.indices;
poly.neighbors Ź NEW[PolygonSequenceRep[i.length]];
poly.neighbors.length Ź i.length;
FOR n: NAT IN [0..i.length) DO
refNatPair­ Ź [i[(n+1) MOD i.length], i[n]];
poly.neighbors[n] Ź NARROW[RefTab.Fetch[neighbors, refNatPair].val];
ENDLOOP;
};
FOR n: NAT IN [0..polygons.length) DO Store[polygons[n]]; ENDLOOP;
FOR n: NAT IN [0..polygons.length) DO SetNeighbor[polygons[n]]; ENDLOOP;
};
Equal: PROC [key1, key2: REF ANY] RETURNS [BOOL] ~ {
RETURN[NARROW[key1, REF NatPair]­ = NARROW[key2, REF NatPair]­];
};

Hash: PROC [key: REF ANY] RETURNS [CARDINAL] ~ {
address: REF NatPair ~ NARROW[key];
RETURN[address.x+address.y];
};
IsConvex: PUBLIC PROC [polygon: Polygon] RETURNS [BOOL] ~ {
GetPair: PROC [n: NAT] RETURNS [Pair] ~ {
RETURN[G3dPlane.ProjectPointToMajorPlane[polygon.points[n], polygon.majorPlane]];
};
IF polygon.indices = NIL OR polygon.indices.length < 3
THEN RETURN[TRUE]
ELSE {
first: BOOL;
p0: Pair Ź GetPair[polygon.indices.length-2];
p1: Pair Ź GetPair[polygon.indices.length-1];
v0: Pair Ź [p1.x-p0.x, p1.y-p0.y];
FOR n: NAT IN [0..polygon.indices.length) DO
p2: Pair Ź GetPair[n];
v1: Pair Ź [p2.x-p1.x, p2.y-p1.y];
cross: REAL Ź Vector2.Cross[v0, v1];
IF n = 0 THEN first Ź cross > 0 ELSE IF first # (cross > 0) THEN RETURN[FALSE];
p0 Ź p1;
p1 Ź p2;
v0 Ź v1;
ENDLOOP;
RETURN[TRUE];
};
};

NPoints: PUBLIC PROC [points: TripleSequence, polygon: NatSequence Ź NIL] RETURNS [NAT]
~ {
RETURN[SELECT TRUE FROM
points = NIL => 0,
polygon # NIL => polygon.length,
ENDCASE => points.length];
};

GetPointN: PROC [n: NAT, points: TripleSequence, polygon: NatSequence Ź NIL]
RETURNS [Triple] ~ INLINE {
RETURN[IF polygon = NIL THEN points[n] ELSE points[polygon[n]]];
};

PolygonFlatness: PUBLIC PROC [
points: TripleSequence,
polygon: NatSequence Ź NIL,
polygonPlane: Plane Ź []]
RETURNS [max: REAL Ź 0.0]
~ {
nPoints: NAT Ź NPoints[points, polygon];
IF nPoints > 3 THEN {
plane: Plane Ź IF polygonPlane = []
THEN G3dPlane.FromPolygon[points, polygon]
ELSE G3dPlane.Unit[polygonPlane];
p1: Triple Ź GetPointN[nPoints-1, points, polygon];
FOR n: NAT IN [0..nPoints) DO
p0: Triple Ź p1;
v: Triple Ź G3dVector.Unit[
G3dVector.Sub[p1Ź GetPointN[n, points, polygon], p0]];
max Ź MAX[max, ABS[G3dVector.Dot[[plane.x, plane.y, plane.z], v]]];
ENDLOOP;
};
max Ź G2dBasic.ArcCosDeg[1.0-max];
};

PolygonReverse: PUBLIC PROC [polygon: NatSequence] RETURNS [NatSequence] ~ {
IF polygon # NIL THEN FOR n: NAT IN [0..polygon.length/2) DO
nn: NAT Ź polygon.length-n-1;
temp: NAT Ź polygon[n];
polygon[n] Ź polygon[nn];
polygon[nn] Ź temp;
ENDLOOP;
RETURN[polygon];
};
InsidePolygon: PUBLIC PROC [q: Triple, polygon: Polygon] RETURNS [Border] ~ {
RETURN[SELECT polygon.majorPlane FROM
xy => InsidePolygonXY[q, polygon.points, polygon.indices],
xz => InsidePolygonXZ[q, polygon.points, polygon.indices],
yz => InsidePolygonYZ[q, polygon.points, polygon.indices],
ENDCASE => ERROR];
};

InsidePolygonXY: PUBLIC PROC [q: Triple, points: TripleSequence, indices: NatSequence Ź NIL]
RETURNS [Border]
~ {
zcross: REAL;
odd: BOOL Ź FALSE;
nPoints: NAT Ź IF indices # NIL THEN indices.length ELSE points.length;
d2: Pair Ź [points[nPoints-1].x-q.x, points[nPoints-1].y-q.y];
FOR n: NAT IN [0..nPoints) DO
d1: Pair Ź d2;
i2: NAT Ź IF indices # NIL THEN indices[n] ELSE n;
d2 Ź [points[i2].x-q.x, points[i2].y-q.y];
IF (d1.y > 0 AND d2.y > 0) OR (d1.y < 0 AND d2.y < 0) OR (d1.x < 0 AND d2.x < 0)
THEN LOOP;	-- no chance to cross 
zcross Ź d2.y*d1.x-d1.y*d2.x;  
IF zcross = 0.0 AND										   -- point on line containing edge
(d1.x <= 0 OR d2.x <= 0)
THEN RETURN[edge];
IF (d1.y > 0 OR d2.y > 0) AND (zcross < 0) # (d1.y-d2.y < 0)	-- remaining special case
THEN odd Ź NOT odd;
ENDLOOP;
RETURN[IF odd THEN inside ELSE outside];
};

InsidePolygonXZ: PUBLIC PROC [q: Triple, points: TripleSequence, indices: NatSequence Ź NIL]
RETURNS [Border]
~ {
zcross: REAL;
odd: BOOL Ź FALSE;
nPoints: NAT Ź IF indices # NIL THEN indices.length ELSE points.length;
d2: Pair Ź [points[nPoints-1].x-q.x, points[nPoints-1].z-q.z];
FOR n: NAT IN [0..nPoints) DO
d1: Pair Ź d2;
i2: NAT Ź IF indices # NIL THEN indices[n] ELSE n;
d2 Ź [points[i2].x-q.x, points[i2].z-q.z];
IF (d1.y >0 AND d2.y >0) OR (d1.y <0 AND d2.y <0) OR (d1.x <0 AND d2.x <0) THEN LOOP; 
zcross Ź d2.y*d1.x-d1.y*d2.x;  
IF zcross = 0.0 AND (d1.x <= 0 OR d2.x <= 0) THEN RETURN[edge];
IF (d1.y > 0 OR d2.y > 0) AND (zcross < 0) # (d1.y-d2.y < 0) THEN odd Ź NOT odd;
ENDLOOP;
RETURN[IF odd THEN inside ELSE outside];
};

InsidePolygonYZ: PUBLIC PROC [q: Triple, points: TripleSequence, indices: NatSequence Ź NIL]
RETURNS [Border]
~ {
zcross: REAL;
odd: BOOL Ź FALSE;
nPoints: NAT Ź IF indices # NIL THEN indices.length ELSE points.length;
d2: Pair Ź [points[nPoints-1].y-q.y, points[nPoints-1].z-q.z];
FOR n: NAT IN [0..nPoints) DO
d1: Pair Ź d2;
i2: NAT Ź IF indices # NIL THEN indices[n] ELSE n;
d2 Ź [points[n].y-q.y, points[n].z-q.z];
IF (d1.y >0 AND d2.y >0) OR (d1.y <0 AND d2.y <0) OR (d1.x <0 AND d2.x <0) THEN LOOP; 
zcross Ź d2.y*d1.x-d1.y*d2.x;  
IF zcross = 0.0 AND (d1.x <= 0 OR d2.x <= 0) THEN RETURN[edge];
IF (d1.y > 0 OR d2.y > 0) AND (zcross < 0) # (d1.y-d2.y < 0) THEN odd Ź NOT odd;
ENDLOOP;
RETURN[IF odd THEN inside ELSE outside];
};
TriangleArea: PUBLIC PROC [p0, p1, p2: Triple] RETURNS [REAL] ~ {
RETURN[0.5*G3dVector.Length[G3dVector.Cross[
G3dVector.Sub[p1, p0], G3dVector.Sub[p2, p0]]]];
};

PolygonArea: PUBLIC PROC [
points: TripleSequence,
polygon: NatSequence Ź NIL,
polygonPlane: Plane Ź []]
RETURNS [area: REAL Ź 0.0]
~ {
nPoints: NAT Ź NPoints[points, polygon];
SELECT nPoints FROM
0, 1, 2	=> NULL;
3 => IF polygon = NIL
THEN RETURN[TriangleArea[points[0], points[1], points[2]]]
ELSE {
p0: Triple Ź points[polygon[0]];
p1: Triple Ź points[polygon[1]];
p2: Triple Ź points[polygon[2]];
RETURN[TriangleArea[p0, p1, p2]];
};
ENDCASE	=> {
plane: Plane Ź IF polygonPlane # []
THEN polygonPlane
ELSE G3dPlane.FromPolygon[points, polygon];
majorPlane: MajorPlane Ź G3dPlane.GetMajorPlane[plane];
pair1: Pair Ź G3dPlane.ProjectPointToMajorPlane[
GetPointN[nPoints-1, points, polygon], majorPlane];
correctionFactor: REAL Ź SELECT majorPlane FROM
xy			=> 0.5/plane.z,
xz			=> 0.5/plane.y,
ENDCASE	=> 0.5/plane.x;
FOR n: NAT IN [0..points.length) DO
pair0: Pair Ź pair1;
pair1 Ź G3dPlane.ProjectPointToMajorPlane[
GetPointN[n, points, polygon], majorPlane];
area Ź area+(pair1.x-pair0.x)*(pair0.y+pair1.y);
ENDLOOP;
area Ź correctionFactor*area;
};
};

PolygonAreas: PUBLIC PROC [
points: TripleSequence,
polygons: SurfaceSequence,
areas: RealSequence Ź NIL]
RETURNS [RealSequence]
~ {
IF polygons = NIL OR points = NIL THEN RETURN[NIL];
IF areas = NIL OR areas.length < polygons.length
THEN areas Ź NEW[RealSequenceRep[polygons.length]];
areas.length Ź polygons.length;
FOR n: NAT IN [0..polygons.length) DO
areas[n] Ź PolygonArea[points, polygons[n].vertices];
ENDLOOP;
RETURN[areas];
};
PolygonAngles: PUBLIC PROC [
points: TripleSequence,
polygon: NatSequence Ź NIL,
angles: RealSequence Ź NIL]
RETURNS [RealSequence]
~ {
nPoints: NAT Ź NPoints[points, polygon];
IF nPoints >= 3 THEN {
p0: Triple Ź GetPointN[nPoints-1, points, polygon];
p1: Triple Ź GetPointN[0, points, polygon];
v0: Triple Ź G3dVector.Unit[G3dVector.Sub[p1, p0]];
IF angles = NIL OR angles.length< nPoints THEN angles Ź NEW[RealSequenceRep[nPoints]];
angles.length Ź nPoints;
FOR n: NAT IN [0..nPoints) DO
p2: Triple Ź GetPointN[(n+1) MOD nPoints, points, polygon];
v1: Triple Ź G3dVector.Unit[G3dVector.Sub[p2, p1]];
angles[n] Ź G2dBasic.ArcCos[G3dVector.Dot[v0, v1]];
v0 Ź v1;
p0 Ź p1;
p1 Ź p2;
ENDLOOP;
};
RETURN[angles];
};

TriangleAngles: PUBLIC PROC [p0, p1, p2: Triple] RETURNS [Triple] ~ {
aSqd: REAL ~ G3dVector.SquareDistance[p1, p2];
bSqd: REAL ~ G3dVector.SquareDistance[p0, p2];
cSqd: REAL ~ G3dVector.SquareDistance[p0, p1];
IF aSqd = 0.0 OR bSqd = 0.0 OR cSqd = 0.0
THEN RETURN[[0.0, 0.0, 0.0]]
ELSE {
a: REAL ~ RealFns.SqRt[aSqd];
b: REAL ~ RealFns.SqRt[bSqd];
c: REAL ~ RealFns.SqRt[cSqd];
angleA: REAL ~ G2dBasic.ArcCos[(bSqd+cSqd-aSqd)/(2.0*b*c)];
angleB: REAL ~ G2dBasic.ArcCos[(cSqd+aSqd-bSqd)/(2.0*c*a)];
angleC: REAL ~ 180.0-angleA-angleB;
RETURN[[angleA, angleB, angleC]];
};
};
TriangleCenter: PUBLIC PROC [p0, p1, p2: Triple] RETURNS [Triple] ~ {
RETURN[[(1.0/3.0)*(p0.x+p1.x+p2.x), (1.0/3.0)*(p0.y+p1.y+p2.y), (1.0/3.0)*(p0.z+p1.z+p2.z)]];
};

SolidPolygonCenter: PROC [points: TripleSequence, nPoints: NAT, polygon: NatSequence]
RETURNS [c: Triple] ~ {
ref: ATOM Ź NIL;
TwiceTriangleArea: PROC [p0, p1, p2: Triple] RETURNS [a: REAL] ~ {
t: Triple Ź
G3dVector.Cross[G3dVector.Sub[p1, p0], G3dVector.Sub[p2, p0]];
IF ref = NIL THEN ref Ź SELECT MAX[t.x,t.y,t.z] FROM t.x => $X, t.y => $Y, ENDCASE => $Z;
a Ź G3dVector.Length[t];
IF (SELECT ref FROM $X => t.x, $Y => t.y, ENDCASE => t.z) < 0.0 THEN a Ź -a;
};
ThriceTriangleCenter: PROC [p0, p1, p2: Triple] RETURNS [Triple] ~ {
RETURN[[p0.x+p1.x+p2.x, p0.y+p1.y+p2.y, p0.z+p1.z+p2.z]];
};
p0: Triple Ź GetPointN[0, points, polygon];
p2: Triple Ź GetPointN[1, points, polygon];
twiceArea, twiceAreaSum: REAL Ź 0.0;
thriceCenter, sixTimesWeightedSum: Triple Ź [];
FOR n: NAT IN [2..nPoints) DO
p1: Triple Ź p2;
p2 Ź GetPointN[n, points, polygon];
twiceArea Ź ABS[TwiceTriangleArea[p0, p1, p2]];
thriceCenter Ź ThriceTriangleCenter[p0, p1, p2];
sixTimesWeightedSum Ź [
sixTimesWeightedSum.x+twiceArea*thriceCenter.x,
sixTimesWeightedSum.y+twiceArea*thriceCenter.y,
sixTimesWeightedSum.z+twiceArea*thriceCenter.z];
twiceAreaSum Ź twiceAreaSum+twiceArea;
ENDLOOP;
IF twiceAreaSum = 0.0
THEN RETURN[WirePolygonCenter[points, nPoints, polygon]]
ELSE {
f: REAL Ź 1.0/(3.0*twiceAreaSum);
c Ź [f*sixTimesWeightedSum.x, f*sixTimesWeightedSum.y, f*sixTimesWeightedSum.z];
};
};

WirePolygonCenter: PROC [points: TripleSequence, nPoints: NAT, polygon: NatSequence]
RETURNS [v: Triple] ~ {
totalLength: REAL Ź 0.0;
point1: Triple Ź GetPointN[nPoints-1, points, polygon];
FOR n: NAT IN [0..nPoints) DO
point0: Triple Ź point1;
length: REAL Ź
G3dVector.Distance[point0, point1 Ź GetPointN[n, points, polygon]];
v.x Ź v.x+length*(point0.x+point1.x);
v.y Ź v.y+length*(point0.y+point1.y);
v.z Ź v.z+length*(point0.z+point1.z);
totalLength Ź totalLength+length;
ENDLOOP;
v Ź G3dVector.Div[v, totalLength+totalLength];
};

PolygonCenter: PUBLIC PROC [
points: TripleSequence,
polygon: NatSequence Ź NIL,
solid: BOOL Ź TRUE]
RETURNS [Triple]
~ {
nPoints: NAT Ź NPoints[points, polygon];
SELECT nPoints FROM
< 3 => RETURN[[]];
3 => IF polygon = NIL
THEN RETURN[TriangleCenter[points[0], points[1], points[2]]]
ELSE {
p0: Triple Ź points[polygon[0]];
p1: Triple Ź points[polygon[1]];
p2: Triple Ź points[polygon[2]];
RETURN[TriangleCenter[p0, p1, p2]];
};
ENDCASE => IF solid
THEN RETURN[SolidPolygonCenter[points, nPoints, polygon]]
ELSE RETURN[WirePolygonCenter[points, nPoints, polygon]];
};

PolygonCenters: PUBLIC PROC [
points: TripleSequence,
polygons: SurfaceSequence,
solid: BOOL Ź TRUE,
centers: TripleSequence Ź NIL]
RETURNS [TripleSequence]
~ {
IF polygons = NIL OR points = NIL THEN RETURN[NIL];
IF centers = NIL OR centers.maxLength < polygons.length
THEN centers Ź NEW[TripleSequenceRep[polygons.length]];
centers.length Ź polygons.length;
FOR n: NAT IN [0..polygons.length) DO
centers[n] Ź PolygonCenter[points, polygons[n].vertices, solid];
ENDLOOP;
RETURN[centers];
};
TriangleNormal: PUBLIC PROC [p0, p1, p2: Triple, unitize: BOOL Ź FALSE]
RETURNS [normal: Triple]
~ {
v1: Triple Ź [p1.x-p0.x, p1.y-p0.y, p1.z-p0.z];
v2: Triple Ź [p2.x-p1.x, p2.y-p1.y, p2.z-p1.z];
normal Ź [v1.z*v2.y-v1.y*v2.z, v1.x*v2.z-v1.z*v2.x, v1.y*v2.x-v1.x*v2.y]; -- right-handed
IF unitize THEN normal Ź G3dVector.Unit[normal];
};

PolygonNormal: PUBLIC PROC [
points: TripleSequence,
polygon: NatSequence Ź NIL,
normalize: BOOL Ź FALSE]
RETURNS [normal: Triple Ź []]
~ {
nPoints: NAT Ź NPoints[points, polygon];
SELECT nPoints FROM
0, 1, 2 => RETURN;
3 => IF polygon = NIL
THEN normal Ź TriangleNormal[points[0], points[1], points[2]]
ELSE {
p0: Triple Ź points[polygon[0]];
p1: Triple Ź points[polygon[1]];
p2: Triple Ź points[polygon[2]];
normal Ź TriangleNormal[p0, p1, p2];
};
ENDCASE	=> {
p1: Triple Ź GetPointN[nPoints-1, points, polygon];
FOR n: NAT IN[0..nPoints) DO								-- Newells' method
p0: Triple Ź p1;
p1 Ź GetPointN[n, points, polygon];
normal Ź [		-- negated on June 22, 1988, for new RH coord sys.
normal.x+(p1.y-p0.y)*(p0.z+p1.z),
normal.y+(p1.z-p0.z)*(p0.x+p1.x),
normal.z+(p1.x-p0.x)*(p0.y+p1.y)];
ENDLOOP;
};
IF normalize THEN normal Ź G3dVector.Unit[normal];
};

PolygonNormals: PUBLIC PROC [
points: TripleSequence,
polygons: SurfaceSequence,
normals: TripleSequence Ź NIL,
unitize: BOOL Ź FALSE]
RETURNS [TripleSequence]
~ {
IF polygons # NIL AND points # NIL THEN {
IF normals = NIL OR normals.maxLength < polygons.length
THEN normals Ź NEW[TripleSequenceRep[polygons.length]];
normals.length Ź polygons.length;
FOR n: NAT IN [0..polygons.length) DO
normals[n] Ź PolygonNormal[points, polygons[n].vertices, unitize];
ENDLOOP;
};
RETURN[normals];
};
NearestToTriangle: PUBLIC PROC [triangle: Triangle, q: Triple] RETURNS [near: NearTriangle]
~ {
sense: BOOL;
e0, e1, e2: REAL;
qq: Triple ~ G3dPlane.ProjectPointToPlane[q, triangle.plane];
pair: Pair Ź G3dPlane.ProjectPointToMajorPlane[qq, triangle.majorPlane];
IF triangle.l0 = [] OR triangle.l1 = [] OR triangle.l2 = [] THEN SetTriangle[triangle,, TRUE];
sense Ź (e0 Ź G2dVector.DistanceToLine[pair, triangle.l0]) < 0.0;
near.inside Ź FALSE;
IF ((e1 Ź G2dVector.DistanceToLine[pair, triangle.l1]) < 0.0) # sense THEN RETURN;
IF ((e2 Ź G2dVector.DistanceToLine[pair, triangle.l2]) < 0.0) # sense THEN RETURN;
near Ź [qq, TRUE, e1, e2, e0];
};
NearestToPolygon: PUBLIC PROC [polygon: Polygon, q: Triple] RETURNS [near: NearPolygon] ~ {
qq: Triple ~ G3dPlane.ProjectPointToPlane[q, polygon.plane];
IF (near.inside Ź InsidePolygon[qq, polygon] = inside)
THEN
near.distance Ź G3dVector.Distance[near.point Ź qq, q]
ELSE {
p1: Triple Ź polygon.points[0];
sqDistanceMin: REAL Ź 1000000.0;
FOR n: NAT IN [0..polygon.points.length) DO
p0: Triple Ź p1;
point: Triple Ź G3dVector.NearestToSegment[p0, p1 Ź polygon.points[(n+1) MOD polygon.points.length], q, polygon.accs[n]].point;
sqDistance: REAL Ź G3dVector.SquareDistance[point, q];
IF sqDistance < sqDistanceMin THEN {
sqDistanceMin Ź sqDistance;
near.point Ź point;
};
ENDLOOP;
near.distance Ź RealFns.SqRt[sqDistanceMin];
};
};

NearestTriangleEdge: PUBLIC PROC [triangle: Triangle, q: Triple] RETURNS [Segment] ~ {
d0, d1, d2: REAL;
IF triangle.a0 = [] OR triangle.a1 = [] OR triangle.a2 = [] THEN SetTriangle[triangle, TRUE];
d0 Ź SqDistFromEdge[triangle.p0, triangle.p1, q, triangle.a0];
d1 Ź SqDistFromEdge[triangle.p1, triangle.p2, q, triangle.a1];
d2 Ź SqDistFromEdge[triangle.p2, triangle.p0, q, triangle.a2];
RETURN[SELECT MIN[d0, d1, d2] FROM
d0			=> [triangle.p0, triangle.p1],
d1			=> [triangle.p1, triangle.p2],
ENDCASE	=> [triangle.p2, triangle.p0]];
};

NearestPolygonEdge: PUBLIC PROC [polygon: Polygon, q: Triple] RETURNS [s: Segment] ~ {
min: REAL Ź 100000.0;
p1: Triple Ź polygon.points[0];
IF polygon.accs = NIL THEN SetPolygon[polygon: polygon, accs: TRUE];
FOR n: NAT IN [0..polygon.points.length) DO
nn: NAT Ź (n+1) MOD polygon.points.length;
p0: Triple Ź p1;
sqDist: REAL ~ SqDistFromEdge[p0, p1 Ź polygon.points[nn], q, polygon.accs[n]];
IF sqDist < min THEN {
min Ź sqDist;
s Ź [p0, p1];
};
ENDLOOP;
};

SqDistFromEdge: PROC [q0, q1, q: Triple, acc: Quad] RETURNS [REAL] ~ {
n: G3dVector.NearSegment Ź G3dVector.NearestToSegment[q0, q1, q, acc];
d: Triple ~ [n.point.x-q.x, n.point.y-q.y, n.point.z-q.z];
RETURN[d.x*d.x+d.y*d.y+d.z*d.z];
};
ApplyToPolygons: PUBLIC PROC [polygonProc: PolygonProc, polygons: SurfaceSequence] ~ {
IF polygons # NIL THEN
FOR n: NAT IN[0..polygons.length) DO
IF NOT polygonProc[n, polygons[n].vertices] THEN RETURN;
ENDLOOP;
};

ApplyToFrontFacingPolygons: PUBLIC PROC [
polygonProc: PolygonProc,
polygons: SurfaceSequence,
view: Matrix,
faceCenters: TripleSequence,
normals: TripleSequence]
~ {
ENABLE G3dMatrix.singular => NULL;
IF polygons # NIL AND view # NIL AND normals # NIL THEN {
inverse: Matrix;
persp: BOOL Ź G3dMatrix.HasPerspective[view];
IF persp AND faceCenters = NIL THEN persp Ź FALSE;		-- a kluge
IF persp THEN inverse Ź G3dMatrix.Invert[view, G3dMatrix.ObtainMatrix[]];
FOR n: NAT IN [0..normals.length) DO
visible: BOOL Ź IF persp
THEN G3dVector.FrontFacingWithPerspective[normals[n], faceCenters[n], inverse]
ELSE G3dVector.FrontFacingNoPerspective[normals[n], view];
IF visible THEN {IF NOT polygonProc[n, polygons[n].vertices] THEN RETURN};
ENDLOOP;
FOR n: NAT IN [normals.length..polygons.length) DO
IF NOT polygonProc[n, polygons[n].vertices] THEN RETURN;
ENDLOOP;
IF persp THEN G3dMatrix.ReleaseMatrix[inverse];
};
};
Triangulate2Polygons: PUBLIC PROC [
polygon0, polygon1: NatSequence,
points: TripleSequence]
RETURNS [triangles: Nat3Sequence]
~ {
ProjectToXYPlane: PROC [poly: NatSequence, points: TripleSequence, o, x, y: Triple]
RETURNS [pairs: PairSequence]
~ {
pairs Ź NEW[G3dBasic.PairSequenceRep[poly.length]];
pairs.length Ź poly.length;
FOR n: NAT IN [0..pairs.length) DO
q: Triple ~ G3dVector.Sub[points[poly[n]], o];
pairs[n] Ź [G3dVector.Dot[q, x], G3dVector.Dot[q, y]];
ENDLOOP;
};
NewTriangle: PROC [a, b, c: NAT] ~ {
IF triangles.length >= triangles.maxLength THEN {
old: Nat3Sequence Ź triangles;
triangles Ź NEW[Nat3SequenceRep[Real.Round[1.3*triangles.length]]];
FOR i: NAT IN [0..old.length) DO triangles[i] Ź old[i]; ENDLOOP;
triangles.length Ź old.length;
};
triangles[triangles.length] Ź [a, b, c];
triangles.length Ź triangles.length+1;
};
ScaleAndCenter: PROC ~ {
PutInUnitSquare: PROC [pairs: PairSequence] ~ {
FOR n: NAT IN [0..pairs.length) DO
pairs[n] Ź [sx*(pairs[n].x-cx), sy*(pairs[n].y-cy)];
ENDLOOP;
};
amm: G2dBasic.Box Ź G2dVector.MinMaxSequence[apairs];
bmm: G2dBasic.Box Ź G2dVector.MinMaxSequence[bpairs];
mm: G2dBasic.Box Ź [
[MIN[amm.min.x, bmm.min.x], MIN[amm.min.y, bmm.min.y]],
[MAX[amm.max.x, bmm.max.x], MAX[amm.max.y, bmm.max.y]]];
cx: REAL Ź 0.5*(mm.min.x+mm.max.x);
cy: REAL Ź 0.5*(mm.min.y+mm.max.y);
sx: REAL Ź IF mm.min.x # mm.max.x THEN 2.0/(mm.max.x-mm.min.x) ELSE 0.0;
sy: REAL Ź IF mm.min.y # mm.max.y THEN 2.0/(mm.max.y-mm.min.y) ELSE 0.0;
PutInUnitSquare[apairs];
PutInUnitSquare[bpairs];
};
SetStartingPoints: PROC ~ {
s, max: REAL Ź 10000.0;
FOR n: NAT IN [0..a.length) DO
FOR nn: NAT IN [0..b.length) DO
IF (s Ź G2dVector.SquareDistance[apairs[n], bpairs[nn]]) < max
THEN {astart Ź a0 Ź n; bstart Ź b0 Ź nn; max Ź s};
ENDLOOP;
ENDLOOP;
};
a: NatSequence ~ polygon0;
b: NatSequence ~ PolygonReverse[G2dBasic.CopyNatSequence[polygon1]];
acenter: Triple Ź PolygonCenter[points, a];
bcenter: Triple Ź PolygonCenter[points, b];
center: Triple Ź G3dVector.Mul[G3dVector.Add[acenter, bcenter], 0.5];
z: Triple Ź G3dVector.Unit[G3dVector.Sub[bcenter, acenter]];
x: Triple Ź G3dVector.Ortho[z, G3dBasic.yAxis];
y: Triple Ź G3dVector.Unit[G3dVector.Cross[x, z]];
apairs: PairSequence Ź ProjectToXYPlane[a, points, center, x, y];
bpairs: PairSequence Ź ProjectToXYPlane[b, points, center, x, y];
a0, a1, astart, b0, b1, bstart: NAT Ź 0;
nTriangles: NAT Ź 0;
maxNTriangles: NAT Ź polygon0.length+polygon1.length+2;	-- +2 is for good luck
ScaleAndCenter[];
SetStartingPoints[];
a1 Ź (a0+1) MOD a.length;
b1 Ź (b0+1) MOD b.length;
triangles Ź NEW[Nat3SequenceRep[points.length]];
DO
IF G2dVector.SquareDistance[apairs[a0], bpairs[b1]] <
G2dVector.SquareDistance[apairs[a1], bpairs[b0]]
THEN {
NewTriangle[a[a0], b[b0], b[b1]];
b0 Ź b1;
b1 Ź (b1+1) MOD b.length;
}
ELSE {
NewTriangle[a[a0], b[b0], a[a1]];
a0 Ź a1;
a1 Ź (a1+1) MOD a.length;
};
IF a0 = astart AND b0 = bstart THEN EXIT;
IF (nTriangles Ź nTriangles+1) >= maxNTriangles THEN EXIT;
ENDLOOP;
};
CopyPolygonSequence: PUBLIC PROC [polygons: PolygonSequence]
RETURNS [PolygonSequence] ~ {
copy: PolygonSequence Ź NIL;
IF polygons # NIL THEN {
copy Ź NEW[PolygonSequenceRep[polygons.length]];
copy.length Ź polygons.length;
FOR n: NAT IN [0..polygons.length) DO copy[n] Ź polygons[n]; ENDLOOP;
};
RETURN[copy];
};

AddToPolygonSequence: PUBLIC PROC [polygons: PolygonSequence, polygon: Polygon]
RETURNS [PolygonSequence]
~ {
IF polygons = NIL THEN polygons Ź NEW[PolygonSequenceRep[1]];
IF polygons.length = polygons.maxLength
THEN polygons Ź LengthenPolygonSequence[polygons];
polygons[polygons.length] Ź polygon;
polygons.length Ź polygons.length+1;
RETURN[polygons];
};

LengthenPolygonSequence: PUBLIC PROC [polygons: PolygonSequence, amount: REAL Ź 1.3]
RETURNS [new: PolygonSequence] ~ {
newLength: NAT Ź MAX[Real.Round[amount*polygons.length], 3];
new Ź NEW[PolygonSequenceRep[newLength]];
FOR i: NAT IN [0..polygons.length) DO new[i] Ź polygons[i]; ENDLOOP;
new.length Ź polygons.length;
};

CopyTriangleSequence: PUBLIC PROC [triangles: TriangleSequence]
RETURNS [TriangleSequence] ~ {
copy: TriangleSequence Ź NIL;
IF triangles # NIL THEN {
copy Ź NEW[TriangleSequenceRep[triangles.length]];
copy.length Ź triangles.length;
FOR n: NAT IN [0..triangles.length) DO copy[n] Ź triangles[n]; ENDLOOP;
};
RETURN[copy];
};

AddToTriangleSequence: PUBLIC PROC [triangles: TriangleSequence, triangle: Triangle]
RETURNS [TriangleSequence]
~ {
IF triangles = NIL THEN triangles Ź NEW[TriangleSequenceRep[1]];
IF triangles.length = triangles.maxLength
THEN triangles Ź LengthenTriangleSequence[triangles];
triangles[triangles.length] Ź triangle;
triangles.length Ź triangles.length+1;
RETURN[triangles];
};
LengthenTriangleSequence: PUBLIC PROC [triangles: TriangleSequence, amount: REAL Ź 1.3]
RETURNS [new: TriangleSequence] ~ {
newLength: NAT Ź MAX[Real.Round[amount*triangles.length], 3];
new Ź NEW[TriangleSequenceRep[newLength]];
FOR i: NAT IN [0..triangles.length) DO new[i] Ź triangles[i]; ENDLOOP;
new.length Ź triangles.length;
};

END.

   \  
G3dPolygonImpl.mesa
Copyright Ó 1985, 1992 by Xerox Corporation.  All rights reserved.
Bloomenthal, October 8, 1992 5:37 pm PDT
Glassner, July 5, 1989 7:05:24 pm PDT
Types
Attributes
Store all polygon edges into hash table, keyed by the two vertex indices of that edge
fetch the "values" of each edge from hash table, building neighbor pointer lists

Miscellaneous Procedures
Inside Procedures
Based on code from Pat Hanrahan:
((d1.x>=0 AND d2.x>=0) OR (d1.x<=0 AND d2.x<=0))	-- test needed in case edge is horiz.
((d1.x >= 0) # (d2.x > 0))								-- optimized (fails if point = vertex)
Area Procedures
Angle Procedures
Center Procedures
Normal Procedures
G3dVector.Cross[G3dVector.Sub[p1, p0], G3dVector.Sub[p2, p1];
Nearness Procedures

Callback Procedures
Process polys with no normals (Dorado memory limit):
Triangle Procedures
Sequence Operations

ĘĚ  "cedarcode" styleNewlineDelimiter
 J	Ďe6BJ(J%J JĎk	ktJ ĐblJT[JJ JheadĎlJ%J.J	"J/J%J3J8J'J4J	"J4J/J4J
	'J&J,J
$J	0J1J0J)J6J!;J+J8J"=J1J6 


Ďn#LJOOJ=JFJFJFJJOJOJOĎs˘˘#
˘˘˘RJ%J%J%JĎc"J%J%J%JJJJ Ą
JJJ#0JJ72LJ,EJ/+J1JJJEIJ<@JFTJ);J,JJ+>J-JAJ	*JJ ?J3J3J/6J9EJJJ Ą @J
'J@ĄJJQJJĄ%JJ3J!J,J*DJJJUJ
BJPJ
HJJ Ą4J
@JJ Ą0J	#JJ 
Ą;Ą)JKQJ6JJJ-J-J",JJ"J$JOJJJJJJJJ 
Ą1WJJ	J
 JJJ Ą	1LJJ
@JJ ĄJJJJJJ	(#J&*J!J3JJ6J1CJJJ"JJ ĄL<JJJJJJ
J ĄM%J::J::J::JJJ Ą<\J	JJ JJJ	GJ>JJ2J*

PJŁ!JŁ+@J



Ł&VJ"Ł&HJJ ŁVJJJ
(JJ Ą<\J	JJJJ	GJ>JJ2J*J
	


	ŁVJJ
?J
 PJJ
(JJ Ą<\J	JJJJ	GJ>JJ2J(J
	


	ŁVJJ
?J
 PJJ
(J ĄA&,J0JJ ĄJJJJJJ	(	JJ/:J J J J!J#JJ'+J70J3/JJJ#J*J+J0JJJJJ ĄJJJJJJ
3	0J	#3J%J5JJJ ĄJJJJJJ	(J3J+J3J


VJJ;J3J3JJJJJJ	JJ ĄEJ$.J$.J$.)JJJJJ/;J/;J#J!JJ ĄEJW]JJ Ą#UJJĄBJ>JYJJ
LJĄDJ39JJ+J+J$J/JJ#J /J0J/J/J0J&JJ-8J!JPJJJ Ą#TJJJ7JJCJ%J%J%J!JJ.JJ ĄJJJJ	JJ	(	JJ1<J J J J#JJ.9J.9JJ ĄJJJJJJJ
3$7J%7J!%J@JJ
J 
ĄGJJJ=J/J/JJŁYJ	!0JJ ĄJJJJJJ	(	JJ9=J J J J$JJ3	Ł6JJ#Ł2>J!J!J"JJJ!2JJ ĄJJJJ	JJ

)$7J%7J!%JBJJJ
J Ą [JJJJ==JHJ
^JAJJDRJDRJJJ Ą˘Đks˘	˘˘˘¤˘˘˘˘[J<46J6JJ +JJI3J&6$JJJJJ,JJJ Ą!VJJ
]J>J>J>"J#J#J 'JJ ĄVJJJ$D+J	*JJCOJJJJJJ Ą FJFJ:J J Ą:V$J&8JJJ Ą)JJJJJJJ"9JJ"-J	Ł
>J<I$	JJNJ6:J	&JJJ4#2J&8JJ"/JJ Ą#J  JJ!JĄ=SJJJ(3J"J.J6JJĄ$)1JJ4CJ
@JJJ(J&JĄĄ/"J4JJJ5J5J7J8J#J#J	HJ	HJJJĄJ<>J.2JJJJJDJ+J+JStartOfExpansion![v: Vector3d.Triple, s: REAL]EJ<J/J2JAJAJ (JJ&ŁNJJJ
J
J!035J0J!JJ
JJ!JJ
JJ)J.:JJ Ą<JJJ&0JJ
EJJJJ Ą.OJJJ=%'J.2J$J$JJJ Ą%TJ"J(<J )J
DJJJ Ą?JJJ(2JJ
GJJJJ Ą2TJJJ@')J15J'J&JJJ Ą'WJ#J)=J!*J
FJJJ JJ      g4  \  