DIRECTORY Polygons3d, Real, Vector2, Vector3d; Polygons3dImpl: CEDAR PROGRAM IMPORTS Real, Vector2, Vector3d EXPORTS Polygons3d ~ BEGIN OPEN Polygons3d; TriangleNormal: PUBLIC PROC [p0, p1, p2: Triple] RETURNS [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]; cross: Triple _ [v1.y*v2.z-v1.z*v2.y, v1.z*v2.x-v1.x*v2.z, v1.x*v2.y-v1.y*v2.x]; m: REAL _ cross.x*cross.x+cross.y*cross.y+cross.z*cross.z; RETURN[IF m # 0.0 THEN [cross.x/m, cross.y/m, cross.z/m] ELSE cross]; }; PolygonNormal: PUBLIC PROC [vertices: TripleSequence] RETURNS [Triple] ~ { normal: Triple _ [0.0, 0.0, 0.0]; FOR i: NAT IN[0..vertices.length) DO -- Newells' method p0: Triple _ vertices[i]; p1: Triple _ vertices[(i+1) MOD vertices.length]; normal _ [ normal.x+(p0.y-p1.y)*(p0.z+p1.z), normal.y+(p0.z-p1.z)*(p0.x+p1.x), normal.z+(p0.x-p1.x)*(p0.y+p1.y)]; ENDLOOP; RETURN[normal]; }; PolygonNormals: PUBLIC PROC [ polygons: REF NatTable, vertices: TripleSequence, normals: TripleSequence _ NIL] RETURNS [TripleSequence] ~ { IF polygons # NIL AND vertices # NIL THEN { IF normals = NIL OR normals.length < polygons.length THEN normals _ NEW[TripleSequenceRep[polygons.length]]; normals.length _ polygons.length; FOR n: NAT IN [0..polygons.length) DO poly: REF NatSequence _ polygons[n]; nSides: INTEGER _ poly.length; normals[n] _ [0.0, 0.0, 0.0]; FOR i: NAT IN[0..nSides) DO -- Newells' method p0: Triple _ vertices[poly[i]]; p1: Triple _ vertices[poly[(i+1) MOD nSides]]; normals[n] _ [ normals[n].x+(p0.y-p1.y)*(p0.z+p1.z), normals[n].y+(p0.z-p1.z)*(p0.x+p1.x), normals[n].z+(p0.x-p1.x)*(p0.y+p1.y)]; ENDLOOP; ENDLOOP; }; RETURN[normals]; }; ApplyToFrontFacingPolygons: PUBLIC PROC [ proc: PROC[nPoly: NAT], polygons: REF NatTable, pairs: PairSequence, view: Matrix, normals: TripleSequence] ~ { IF polygons = NIL OR pairs = NIL OR view = NIL OR normals = NIL THEN RETURN; FOR n: NAT IN[0..polygons.length) DO IF normals[n].x*view[0][2]+normals[n].y*view[1][2]+normals[n].z*view[2][2] >= 0.0 THEN proc[n]; ENDLOOP; }; CopyPolygon: PUBLIC PROC [poly: REF NatSequence] RETURNS [copy: REF NatSequence] ~ { copy _ NEW[NatSequence[poly.length]]; copy.length _ poly.length; FOR n: NAT IN [0..poly.length) DO copy[n] _ poly[n]; ENDLOOP; }; ReversePolygon: PUBLIC PROC [poly: REF NatSequence] RETURNS [REF NatSequence] ~ { FOR n: NAT IN [0..poly.length/2) DO nn: NAT _ poly.length-n-1; temp: NAT _ poly[n]; poly[n] _ poly[nn]; poly[nn] _ temp; ENDLOOP; RETURN[poly]; }; CenterOfPolygon: PUBLIC PROC [poly: REF NatSequence, points: TripleSequence] RETURNS [v: Triple _ origin] ~ { FOR n: NAT IN [0..poly.length) DO v _ Vector3d.Add[v, points[poly[n]]]; ENDLOOP; IF poly.length # 0 THEN v _ Vector3d.Div[v, poly.length]; }; CentersOfPolygons: PUBLIC PROC [ polygons: REF NatTable, points: TripleSequence, 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] _ CenterOfPolygon[polygons[n], points]; ENDLOOP; RETURN[centers]; }; ProjectToXYPlane: PUBLIC PROC [ poly: REF NatSequence, points: TripleSequence, o, x, y: Triple] RETURNS [pairs: PairSequence] ~ { pairs _ NEW[PairSequenceRep[poly.length]]; pairs.length _ poly.length; FOR n: NAT IN [0..pairs.length) DO q: Triple ~ Vector3d.Sub[points[poly[n]], o]; pairs[n] _ [Vector3d.Dot[q, x], Vector3d.Dot[q, y]]; ENDLOOP; }; GetMinMaxPairs: PUBLIC PROC [pairs: PairSequence] RETURNS [mm: MinMaxPairs] ~ { mm _ [[10000.0, 10000.0], [-10000.0, -10000.0]]; FOR n: NAT IN [0..pairs.length) DO p: Pair ~ pairs[n]; mm.min _ [MIN[mm.min.x, p.x], MIN[mm.min.y, p.y]]; mm.max _ [MAX[mm.max.x, p.x], MAX[mm.max.y, p.y]]; ENDLOOP; }; GetMinMaxTriples: PUBLIC PROC [triples: TripleSequence] RETURNS [mm: MinMaxTriples] ~ { mm _ [[10000.0, 10000.0, 10000.0], [-10000.0, -10000.0, -10000.0]]; FOR n: NAT IN [0..triples.length) DO p: Triple ~ triples[n]; mm.min _ [MIN[mm.min.x, p.x], MIN[mm.min.y, p.y], MIN[mm.min.z, p.z]]; mm.max _ [MAX[mm.max.x, p.x], MAX[mm.max.y, p.y], MAX[mm.max.z, p.z]]; ENDLOOP; }; SquarePairDistance: PUBLIC PROC [p0, p1: Pair] RETURNS [REAL] ~ { RETURN[Vector2.Square[Vector2.Sub[p0, p1]]]; }; SquareTripleDistance: PUBLIC PROC [p0, p1: Triple] RETURNS [REAL] ~ { RETURN[Vector3d.SquareLength[Vector3d.Sub[p0, p1]]]; }; LengthenTriangleSequence: PUBLIC PROC [t: TriangleSequence] RETURNS [TriangleSequence] ~ { old: TriangleSequence ~ t; t _ NEW[TriangleSequenceRep[Real.RoundI[1.3*old.length]]]; FOR i: NAT IN [0..old.length) DO t[i] _ old[i]; ENDLOOP; t.length _ old.length; RETURN[t]; }; Triangulate2Polygons: PUBLIC PROC [ poly0, poly1: REF NatSequence, points: TripleSequence] RETURNS [triangles: TriangleSequence] ~ { NewTriangle: PROC [a, b, c: NAT] ~ { IF triangles.length >= triangles.maxLength THEN triangles _ LengthenTriangleSequence[triangles]; 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: MinMaxPairs _ GetMinMaxPairs[apairs]; bmm: MinMaxPairs _ GetMinMaxPairs[bpairs]; mm: MinMaxPairs _ [ [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 _ SquarePairDistance[apairs[n], bpairs[nn]]) < max THEN {astart _ a0 _ n; bstart _ b0 _ nn; max _ s}; ENDLOOP; ENDLOOP; }; a: REF NatSequence ~ poly0; b: REF NatSequence ~ ReversePolygon[CopyPolygon[poly1]]; acenter: Triple _ CenterOfPolygon[a, points]; bcenter: Triple _ CenterOfPolygon[b, points]; center: Triple _ Vector3d.Mul[Vector3d.Add[acenter, bcenter], 0.5]; z: Triple _ Vector3d.Normalize[Vector3d.Sub[bcenter, acenter]]; x: Triple _ Vector3d.Ortho[z, yAxis]; y: Triple _ Vector3d.Normalize[Vector3d.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 _ poly0.length+poly1.length+2; -- +2 is for good luck ScaleAndCenter[]; SetStartingPoints[]; a1 _ (a0+1) MOD a.length; b1 _ (b0+1) MOD b.length; triangles _ NEW[TriangleSequenceRep[points.length]]; DO IF SquarePairDistance[apairs[a0], bpairs[b1]] < SquarePairDistance[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; }; END. τPolygons3dImpl.mesa Copyright c 1985 by Xerox Corporation. All rights reserved. Bloomenthal, February 26, 1987 7:27:20 pm PST Normal Procedures RETURN[[Vector3d.Cross[Vector3d.Sub[p1, p0], Vector3d.Sub[p2, p1]]]; RETURN[[ (p0.y*p1.z + p1.y*p2.z + p2.y*p0.z - p0.y*p2.z - p1.y*p0.z - p2.y*p1.z), (p0.z*p1.x + p1.z*p2.x + p2.z*p0.x - p0.z*p2.x - p1.z*p0.x - p2.z*p1.x), (p0.x*p1.y + p1.x*p2.y + p2.x*p0.y - p0.x*p2.y - p1.x*p0.y - p2.x*p1.y)]]; Miscellaneous Procedures Triangle Procedures Κ …˜šœ™Jšœ Οmœ1™™DJšœ/˜/J˜/JšœP˜PJšœžœ3˜:Jšžœžœ žœ#žœ˜Ešžœ™J™IJ™IJ™KJ™—J˜J˜—š‘ œž œžœ ˜JJ˜!J˜š žœžœžœžœΟc˜>J˜Jšœžœ˜1˜ J˜!J˜!J˜"—Jšžœ˜J˜—Jšžœ ˜J˜J˜—š‘œžœžœ˜Jšœ žœ ˜Jšœ˜Jšœžœ˜Jšžœ˜Jšœ˜š žœ žœžœ žœžœ˜+šžœ žœžœ!˜4Jšžœ žœ%˜7—Jšœ!˜!J˜šžœžœžœž˜%Jšœžœ˜$Jšœžœ˜J˜J˜J˜š žœžœžœ žœ’˜5J˜Jšœ!žœ ˜.˜J˜%J˜%J˜&—Jšžœ˜—Jšžœ˜—J˜—Jšžœ ˜J˜J˜—š‘œž œ˜)Jšœžœžœ˜Jšœ žœ ˜Jšœ˜Jšœ ˜ Jšœ˜Jšœ˜Jšžœ žœžœ žœžœžœžœ žœžœžœ˜Lšžœžœžœž˜$šžœO˜QJšžœ ˜ —Jšžœ˜—J˜——š ™š ‘ œžœžœžœžœžœ˜TJšœžœ˜%J˜šžœžœžœž˜!J˜Jšžœ˜—J˜J˜—š ‘œž œžœžœžœ˜Qšžœžœžœž˜#Jšœžœ˜Jšœžœ ˜J˜J˜Jšžœ˜—Jšžœ˜ J˜J˜—š‘œž œžœ%˜LJšžœ˜ šžœžœžœž˜!Jšœ%˜%Jšžœ˜—Jšžœžœ"˜9J˜J˜—š‘œžœžœ˜ Jšœ žœ ˜Jšœ˜Jšœžœ˜Jšžœ˜Jšœ˜Jšžœ žœžœ žœžœžœžœ˜3šžœ žœžœ$˜7Jšžœ žœ%˜7—Jšœ!˜!šžœžœžœž˜%Jšœ2˜2Jšžœ˜—Jšžœ ˜J˜J˜—š‘œžœžœ˜Jšœžœ6˜?Jšžœ˜Jšœ˜Jšœžœ˜*J˜šžœžœžœž˜"J˜-Icodešœ4˜4Lšžœ˜—J˜J˜—š‘œž œžœ˜OL˜0šžœžœžœž˜"L˜Lšœ žœžœ˜2Lšœ žœžœ˜2Lšžœ˜—L˜L˜—š‘œž œžœ˜WL˜Cšžœžœžœž˜$Lšœ˜Lšœ žœžœžœ˜FLšœ žœžœžœ˜FLšžœ˜—L˜L˜—š‘œž œžœžœ˜AJšžœ&˜,Jšœ˜J˜—š‘œž œžœžœ˜EJšžœ.˜4Jšœ˜——š ™š‘œž œ˜;Jšžœ˜J˜Jšœžœ3˜:Jš žœžœžœžœžœ˜8J˜Jšžœ˜ J˜J˜—š‘œžœžœ˜#Jšœžœ%˜6Jšžœ˜%Jšœ˜š‘ œžœ žœ˜$šžœ(˜*Jšžœ1˜5—Jšœ(˜(J˜&J˜—š‘œžœ˜š‘œžœ˜/šžœžœžœž˜"Jšœ4˜4Jšžœ˜—J˜—L˜*L˜*šœ˜Lšœžœžœ˜7Lšœžœžœ˜8—Lšœžœ˜#Lšœžœ˜#Lš œžœžœžœžœ˜HLš œžœžœžœžœ˜HJšœ˜Jšœ˜J˜—š‘œžœ˜Jšœžœ ˜šžœžœžœž˜šžœžœžœž˜šžœ6˜8Jšžœ.˜2—Jšžœ˜—Jšžœ˜—J˜—Jšœžœ˜Jšœžœ2˜8J˜-J˜-J•StartOfExpansion![v: Vector3d.Triple, s: REAL]˜CJ˜?J˜%J˜5JšœA˜AJšœA˜AJšœ žœ˜(Jšœ žœ˜Jšœžœ ’˜HJ˜J˜J˜Jšœ žœ ˜Jšœ žœ ˜Jšœ žœ%˜4šž˜šžœX˜Zšžœ˜Jšœ!˜!J˜Jšœ žœ ˜J˜—šžœ˜Jšœ!˜!J˜Jšœ žœ ˜J˜—Jšžœ žœ žœžœ˜)Jšžœ.žœžœ˜:—Jšžœ˜—J˜—L˜—Jšžœ˜J˜J˜—…—*'£