<<>> <> <> <> DIRECTORY G2dBasic, G2dVector, G3dBasic, G3dPatch, G3dPlane, G3dVector, Real, RealFns, G3dTriangle; G3dTriangleImpl: CEDAR MONITOR IMPORTS G2dBasic, G2dVector, G3dPatch, G3dPlane, G3dVector, Real, RealFns EXPORTS G3dTriangle ~ BEGIN <> Triple: TYPE ~ G3dBasic.Triple; Patch: TYPE ~ G3dPatch.Patch; Plane: TYPE ~ G3dPlane.Plane; IntTriple: TYPE ~ G3dTriangle.IntTriple; IntTriples: TYPE ~ G3dTriangle.IntTriples; Segment: TYPE ~ G3dTriangle.Segment; Segments: TYPE ~ G3dTriangle.Segments; Triangle: TYPE ~ G3dTriangle.Triangle; Triangles: TYPE ~ G3dTriangle.Triangles; TriangleProc: TYPE ~ G3dTriangle.TriangleProc; TriIntersection: TYPE ~ G3dTriangle.TriIntersection; <> TrianglesFromPatch: PUBLIC PROC [patch: Patch, sRes, tRes: INT] RETURNS [triangles: REF Triangles] ~ { AddTriangle: PROC [p1, p2, p3: Triple] ~ { triangles[triangles.length] ¬ MakeTriangle[p1, p2, p3]; triangles.length ¬ triangles.length+1; }; triangles ¬ NEW[Triangles[2*sRes*tRes]]; FOR s: INT IN [0..sRes) DO s0: REAL ¬ REAL[s]/REAL[sRes]; s1: REAL ¬ REAL[s+1]/REAL[sRes]; FOR t: INT IN [0..tRes) DO t0: REAL ¬ REAL[t]/REAL[tRes]; t1: REAL ¬ REAL[t+1]/REAL[tRes]; p00: Triple ¬ G3dPatch.Position[patch, s0, t0]; p01: Triple ¬ G3dPatch.Position[patch, s0, t1]; p10: Triple ¬ G3dPatch.Position[patch, s1, t0]; p11: Triple ¬ G3dPatch.Position[patch, s1, t1]; AddTriangle[p00, p01, p10]; AddTriangle[p10, p01, p11]; ENDLOOP; ENDLOOP; }; BoundaryFromPatch: PUBLIC PROC [patch: Patch, sRes, tRes: INT, cullColinearPoints: BOOL ¬FALSE] RETURNS [segments: REF Segments] ~ { <> AddSeg: PROC [s, t: REAL] ~ { pp: Triple ¬ G3dPatch.Position[patch, s, t]; segments ¬ AddToSegments[segments, [p, pp, G3dVector.Distance[p, pp]]]; p ¬ pp; }; p: Triple ¬ G3dPatch.Position[patch, 0.0, 0.0]; segments ¬ NEW[Segments[2*(sRes+tRes)]]; FOR t: INT IN [1..tRes] DO -- s = 0, t from 0 to 1 AddSeg[0.0, REAL[t]/REAL[tRes]]; ENDLOOP; FOR s: INT IN [1..sRes] DO -- s from 0 to 1, t = 1 AddSeg[REAL[s]/REAL[sRes], 1.0]; ENDLOOP; FOR t: INT DECREASING IN [0..tRes-1] DO -- t s = 1, from 1 to 0 AddSeg[1.0, REAL[t]/REAL[tRes]]; ENDLOOP; FOR s: INT DECREASING IN [0..sRes-1] DO -- s from 1 to 0, t = 0 AddSeg[REAL[s]/REAL[sRes], 0.0]; ENDLOOP; IF cullColinearPoints THEN CullColinearPoints[segments]; }; CullColinearPoints: PROC [segments: REF Segments] ~ { FOR i1: INT ¬ 0, i1+1 WHILE i1 < INT[segments.length] DO i2: INT ¬ (i1+1) MOD segments.length; s1: Segment ¬ segments[i1]; s2: Segment ¬ segments[i2]; IF G3dVector.Collinear[s1.p1, s1.p2, s2.p2] THEN { segments[i1].p2 ¬ s2.p2; segments.length ¬ segments.length-1; FOR i: INT IN [i2..segments.length) DO segments[i] ¬ segments[i+1]; ENDLOOP; }; ENDLOOP; }; AddToSegments: PROC [segments: REF Segments, s: Segment] RETURNS [REF Segments] ~ { LengthenSegments: PROC [segments: REF Segments] RETURNS [new: REF Segments] ~ { newLength: NAT ¬ MAX[Real.Ceiling[1.3*segments.maxLength], 3]; new ¬ NEW[Segments[newLength]]; FOR i: NAT IN [0..segments.length) DO new[i] ¬ segments[i]; ENDLOOP; new.length ¬ segments.length; }; IF segments = NIL THEN segments ¬ NEW[Segments[1]]; IF segments.length = segments.maxLength THEN segments ¬ LengthenSegments[segments]; segments[segments.length] ¬ s; segments.length ¬ segments.length+1; RETURN[segments]; }; <<>> DistanceToTriangles: PUBLIC PROC [triangles: REF Triangles, p: Triple] RETURNS [minD: REAL] ~ { minAbsD, minSqD: REAL ¬ 100000.0; FOR i: INT IN [0..triangles.length) DO t: REF Triangle ¬ triangles[i]; d: REAL ¬ G3dPlane.DistanceToPoint[p, t.plane, FALSE]; absD: REAL ¬ ABS[d]; IF absD < minAbsD THEN { q: Triple ~ G3dPlane.ProjectPointToPlane[p, t.plane]; IF InsideTriangle[q, t] THEN {minAbsD ¬ absD; minD ¬ d; minSqD ¬ d*d} ELSE { Test: PROC [p1, p2: Triple] ~ { near: Triple ¬ G3dVector.NearestToSegment[p1, p2, p].point; sqD: REAL ¬ G3dVector.SquareDistance[near, p]; IF sqD < minSqD THEN { minSqD ¬ sqD; minAbsD ¬ RealFns.SqRt[sqD]; minD ¬ IF d < 0.0 THEN -minAbsD ELSE minAbsD; }; }; Test[t.p1, t.p2]; Test[t.p2, t.p3]; Test[t.p3, t.p1]; }; }; ENDLOOP; }; TriangleCenter: PUBLIC PROC [t: REF Triangle] RETURNS [c: Triple] ~ { c ¬ G3dVector.Midpoint[G3dVector.Midpoint[t.p1, t.p2], G3dVector.Midpoint[t.p2, t.p3]]; }; AverageTriangleNormals: PUBLIC PROC [triangles: REF Triangles] RETURNS [n: Triple ¬ []] ~ { FOR i: INT IN [0..triangles.length) DO n ¬ G3dVector.Add[n, G3dPlane.NormalFromPlane[triangles[i].plane]]; ENDLOOP; n ¬ G3dVector.Unit[G3dVector.Div[n, triangles.length]]; }; MaxDeviationFromNormal: PUBLIC PROC [normal: Triple, triangles: REF Triangles] RETURNS [deg: REAL ¬ 1.0] ~ { FOR i: INT IN [0..triangles.length) DO deg ¬ MIN[deg, G3dVector.Dot[normal, G3dPlane.NormalFromPlane[triangles[i].plane]]]; ENDLOOP; deg ¬ G2dBasic.ArcCosDeg[deg]; }; <<>> ApplyToTrianglesSides: PUBLIC PROC [triangles: REF Triangles, action: TriangleProc] ~ { FOR i: INT IN [0..triangles.length) DO t: REF Triangle ¬ triangles[i]; IF NOT action[t.p1, t.p2, t] OR NOT action[t.p2, t.p3, t] OR NOT action[t.p3, t.p1, t] THEN EXIT; ENDLOOP; }; MakeTriangle: PUBLIC PROC [p1, p2, p3: Triple, scratch: REF Triangle ¬ NIL] RETURNS [t: REF Triangle] ~ { t ¬ IF scratch = NIL THEN NEW[Triangle] ELSE scratch; t.p1 ¬ p1; t.p2 ¬ p2; t.p3 ¬ p3; SetTriangle[t]; }; SetTriangle: PUBLIC PROC [t: REF Triangle] ~ { q1, q2, q3: G3dBasic.Pair; t.plane ¬ G3dPlane.FromThreePoints[t.p1, t.p2, t.p3, TRUE]; t.majorPlane ¬ G3dPlane.GetMajorPlane[t.plane]; q1 ¬ G3dPlane.ProjectPointToMajorPlane[t.p1, t.majorPlane]; q2 ¬ G3dPlane.ProjectPointToMajorPlane[t.p2, t.majorPlane]; q3 ¬ G3dPlane.ProjectPointToMajorPlane[t.p3, t.majorPlane]; IF (SELECT t.majorPlane FROM xy => t.plane.z, xz => t.plane.y, ENDCASE => t.plane.x) > 0.0 THEN { t.l1 ¬ G2dVector.Line[q1, q2]; t.l2 ¬ G2dVector.Line[q2, q3]; t.l3 ¬ G2dVector.Line[q3, q1]; } ELSE { -- invert from normal sense t.l1 ¬ G2dVector.Line[q2, q1]; t.l2 ¬ G2dVector.Line[q3, q2]; t.l3 ¬ G2dVector.Line[q1, q3]; }; }; SetTriangles: PUBLIC PROC [triangles: REF Triangles] ~ { FOR i: INT IN [0..triangles.length) DO SetTriangle[triangles[i]]; ENDLOOP; }; InsideTriangle: PUBLIC PROC [p: Triple, t: REF Triangle] RETURNS [b: BOOL] ~ { pair: G3dBasic.Pair ¬ G3dPlane.ProjectPointToMajorPlane[p, t.majorPlane]; <> <> <> <> <> b ¬ G2dVector.DistanceToLine[pair, t.l1] > -0.0001 AND G2dVector.DistanceToLine[pair, t.l2] > -0.0001 AND G2dVector.DistanceToLine[pair, t.l3] > -0.0001; }; IntersectTriangle: PUBLIC PROC [p1, p2: Triple, t: REF Triangle] RETURNS [i: TriIntersection] ~ { d1: REAL ¬ p1.x*t.plane.x+p1.y*t.plane.y+p1.z*t.plane.z+t.plane.w; d2: REAL ¬ p2.x*t.plane.x+p2.y*t.plane.y+p2.z*t.plane.z+t.plane.w; i.intersect ¬ IF (d1 > 0) = (d2 > 0) THEN FALSE ELSE InsideTriangle[i.point ¬ G3dVector.Interp[d1/(d1-d2), p1, p2], t]; }; Normal: PUBLIC PROC [p1, p2, p3: Triple] RETURNS [n: Triple] ~ { n ¬ G3dVector.Cross[G3dVector.Sub[p2, p1], G3dVector.Sub[p3, p1]]; }; scratchTriangles: ARRAY [0..10) OF REF Triangles ¬ ALL[NIL]; ObtainTriangles: PUBLIC ENTRY PROC RETURNS [REF Triangles] ~ { FOR i: INT IN [0..10) DO scratch: REF Triangles ¬ scratchTriangles[i]; IF scratch # NIL THEN {scratchTriangles[i] ¬ NIL; scratch.length ¬ 0; RETURN[scratch]}; ENDLOOP; RETURN[NEW[Triangles[12]]]; }; ReleaseTriangles: PUBLIC ENTRY PROC [scratch: REF Triangles] ~ { FOR i: INT IN [0..10) DO IF scratchTriangles[i] = NIL THEN {scratchTriangles[i] ¬ scratch; EXIT}; ENDLOOP; }; AddTriangle: PUBLIC PROC [t: REF Triangle, triangles: REF Triangles] RETURNS [r: REF Triangles] ~ { LengthenTriangles: PROC [triangles: REF Triangles] RETURNS [new: REF Triangles] ~ { newLength: NAT ¬ MAX[Real.Ceiling[1.3*triangles.maxLength], 3]; new ¬ NEW[Triangles[newLength]]; FOR i: NAT IN [0..triangles.length) DO new[i] ¬ triangles[i]; ENDLOOP; new.length ¬ triangles.length; }; IF (r ¬ triangles) = NIL THEN r ¬ NEW[Triangles[1]]; IF r.length = r.maxLength THEN r ¬ LengthenTriangles[r]; r[r.length] ¬ t; r.length ¬ r.length+1; }; CopyTriangles: PUBLIC PROC [triangles: REF Triangles] RETURNS [out: REF Triangles] ~ { out ¬ NEW[Triangles[triangles.length]]; FOR i: INT IN [0..triangles.length) DO out[i] ¬ triangles[i]; ENDLOOP; out.length ¬ triangles.length; }; AddIntTriple: PUBLIC PROC [i: IntTriple, intTriples: REF IntTriples] RETURNS [r: REF IntTriples] ~ { LengthenIntTriples: PROC [ints: REF IntTriples] RETURNS [new: REF IntTriples] ~ { newLength: NAT ¬ MAX[Real.Ceiling[1.3*ints.maxLength], 3]; new ¬ NEW[IntTriples[newLength]]; FOR i: NAT IN [0..ints.length) DO new[i] ¬ ints[i]; ENDLOOP; new.length ¬ ints.length; }; IF (r ¬ intTriples) = NIL THEN r ¬ NEW[IntTriples[1]]; IF r.length = r.maxLength THEN r ¬ LengthenIntTriples[r]; r[r.length] ¬ i; r.length ¬ r.length+1; }; END.