<<>> <> <> <> <> 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 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]; }; <> 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.