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
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
Types
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;
Attributes
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;
};
Store all polygon edges into hash table, keyed by the two vertex indices of that edge
FOR n: NAT IN [0..polygons.length) DO Store[polygons[n]]; ENDLOOP;
fetch the "values" of each edge from hash table, building neighbor pointer lists
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];
};
Miscellaneous Procedures
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];
};
Inside Procedures
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]
~ {
Based on code from Pat Hanrahan:
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 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)
(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];
};
Area Procedures
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];
};
Angle Procedures
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]];
};
};
Center Procedures
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];
};
Normal Procedures
TriangleNormal: PUBLIC PROC [p0, p1, p2: Triple, unitize: BOOL ¬ FALSE]
RETURNS [normal: Triple]
~ {
G3dVector.Cross[G3dVector.Sub[p1, p0], G3dVector.Sub[p2, p1];
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];
};
Nearness Procedures
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];
};
Callback Procedures
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;
Process polys with no normals (Dorado memory limit):
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];
};
};
Triangle Procedures
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;
};
Sequence Operations
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.