<<>> <> <> <> DIRECTORY G3dBasic, G3dPlane, G3dTriangle, G3dVector, Real, RealFns, Rope, ImplicitCubeTet; ImplicitCubeTetImpl: CEDAR MONITOR IMPORTS G3dPlane, G3dTriangle, G3dVector, Real, RealFns EXPORTS ImplicitCubeTet ~ BEGIN <> Triple: TYPE ~ G3dBasic.Triple; Plane: TYPE ~ G3dPlane.Plane; IntTriple: TYPE ~ G3dTriangle.IntTriple; Segment: TYPE ~ G3dTriangle.Segment; Segments: TYPE ~ G3dTriangle.Segments; Triangle: TYPE ~ G3dTriangle.Triangle; Triangles: TYPE ~ G3dTriangle.Triangles; ROPE: TYPE ~ Rope.ROPE; Cube: TYPE ~ ImplicitCubeTet.Cube; Corner: TYPE ~ ImplicitCubeTet.Corner; EventType: TYPE ~ ImplicitCubeTet.EventType; Event: TYPE ~ ImplicitCubeTet.Event; Events: TYPE ~ ImplicitCubeTet.Events; TetEdge: TYPE ~ ImplicitCubeTet.TetEdge; TetFace: TYPE ~ ImplicitCubeTet.TetFace; TetIntersection: TYPE ~ ImplicitCubeTet.TetIntersection; TetIntersections: TYPE ~ ImplicitCubeTet.TetIntersections; TetIntersectProc: TYPE ~ ImplicitCubeTet.TetIntersectProc; TetProc: TYPE ~ ImplicitCubeTet.TetProc; Tetrahedron: TYPE ~ ImplicitCubeTet.Tetrahedron; ThreeEdges: TYPE ~ ImplicitCubeTet.ThreeEdges; ThreeCorners: TYPE ~ ImplicitCubeTet.ThreeCorners; TwoCorners: TYPE ~ ImplicitCubeTet.TwoCorners; sqRt3: REAL Ź RealFns.SqRt[3.0]; <> MidCorner: PUBLIC PROC [c1, c2: REF Corner] RETURNS [p: Triple] ~ { p.x Ź IF c1.id.i # c2.id.i THEN 0.5*(c1.pos.x+c2.pos.x) ELSE c1.pos.x; p.y Ź IF c1.id.j # c2.id.j THEN 0.5*(c1.pos.y+c2.pos.y) ELSE c1.pos.y; p.z Ź IF c1.id.k # c2.id.k THEN 0.5*(c1.pos.z+c2.pos.z) ELSE c1.pos.z; }; InterpCorners: PUBLIC PROC [t: REAL, c1, c2: REF Corner] RETURNS [p: Triple] ~ { p.x Ź IF c1.id.i # c2.id.i THEN c1.pos.x+t*(c2.pos.x-c1.pos.x) ELSE c1.pos.x; p.y Ź IF c1.id.j # c2.id.j THEN c1.pos.y+t*(c2.pos.y-c1.pos.y) ELSE c1.pos.y; p.z Ź IF c1.id.k # c2.id.k THEN c1.pos.z+t*(c2.pos.z-c1.pos.z) ELSE c1.pos.z; }; <> IntersectWithSurface: PUBLIC PROC [tet: REF Tetrahedron, action: TetIntersectProc] ~ { Pos: PROC [c: REF Corner, i: INT] RETURNS [j: INT] ~ {j Ź IF c.state = positive THEN i ELSE 0}; << [Artwork node; type 'Artwork on' to command tool] >> SELECT (Pos[tet.a, 8]+Pos[tet.b, 4]+Pos[tet.c, 2]+Pos[tet.d, 1]) FROM 1 => action[3, bd, cd, ad, ab]; 2 => action[3, ac, cd, bc, ab]; 3 => action[4, ad, bd, bc, ac]; 4 => action[3, ab, bc, bd, ab]; 5 => action[4, ad, ab, bc, cd]; 6 => action[4, ab, ac, cd, bd]; 7 => action[3, ab, ac, ad, ab]; 8 => action[3, ab, ad, ac, ab]; 9 => action[4, ab, bd, cd, ac]; 10 => action[4, ab, ad, cd, bc]; 11 => action[3, ab, bd, bc, ab]; 12 => action[4, ad, ac, bc, bd]; 13 => action[3, cd, ac, bc, ab]; 14 => action[3, bd, ad, cd, ab]; ENDCASE; -- 0, 15 }; FaceName: PUBLIC PROC [f: TetFace] RETURNS [r: ROPE] ~ { r Ź SELECT f FROM abc=>"abc", abd=>"abd", acd=>"acd", ENDCASE=>"bcd"; }; EdgeName: PUBLIC PROC [e: TetEdge] RETURNS [r: ROPE] ~ { r Ź SELECT e FROM ab=>"ab",ac=>"ac",ad=>"ad",bc=>"bc",bd=>"bd",ENDCASE=>"cd"; }; CornerName: PUBLIC PROC [c: REF Corner, tet: REF Tetrahedron] RETURNS [r: ROPE] ~ { r Ź SELECT c FROM tet.a => "a", tet.b => "b", tet.c => "c", tet.d => "d", ENDCASE => "?"; }; MakeTetrahedron: PUBLIC PROC [ a, b, c, d: REF Corner, which: INT, scratch: REF Tetrahedron Ź NIL] RETURNS [tet: REF Tetrahedron] ~ { IF (tet Ź scratch) = NIL THEN tet Ź NEW[Tetrahedron]; tet.which Ź which; tet.a Ź a; tet.b Ź b; tet.c Ź c; tet.d Ź d; tet.triangles[abc] Ź G3dTriangle.MakeTriangle[a.pos, c.pos, b.pos, tet.triangles[abc]]; tet.triangles[abd] Ź G3dTriangle.MakeTriangle[a.pos, b.pos, d.pos, tet.triangles[abd]]; tet.triangles[acd] Ź G3dTriangle.MakeTriangle[a.pos, d.pos, c.pos, tet.triangles[acd]]; tet.triangles[bcd] Ź G3dTriangle.MakeTriangle[b.pos, c.pos, d.pos, tet.triangles[bcd]]; }; <<>> TrianglesInTetrahedron: PUBLIC PROC [ candidates: REF Triangles, tet: REF Tetrahedron, scratch: REF Triangles] RETURNS [ret: REF Triangles] ~ { IF (ret Ź scratch) # NIL THEN ret.length Ź 0; FOR i: INT IN [0..candidates.length) DO t: REF Triangle Ź candidates[i]; IF TriangleInTetrahedron[t, tet] THEN ret Ź G3dTriangle.AddTriangle[t, ret]; ENDLOOP; }; TriangleInTetrahedron: PUBLIC PROC [t: REF Triangle, tet: REF Tetrahedron] RETURNS [BOOL] ~ { <> PlaneInTetrahedron: PROC [p: Plane] RETURNS [b: BOOL] ~ { Pos: PROC [c: REF Corner, type: CornerType] RETURNS [b: BOOL] ~ { d: REAL Ź c.pos.x*p.x+c.pos.y*p.y+c.pos.z*p.z+p.w; cInfo[type] Ź [TRUE, d, b Ź d > 0.0, c]; }; s: BOOL Ź Pos[tet.a, a]; b Ź Pos[tet.b, b] # s OR Pos[tet.c, c] # s OR Pos[tet.d, d] # s; }; <> CornerInfo: TYPE ~ RECORD [set: BOOL, distance: REAL, pos: BOOL, c: REF Corner]; CornerType: TYPE ~ {a, b, c, d}; cInfo: ARRAY CornerType OF CornerInfo Ź ALL[[FALSE, 0.0, FALSE, NIL]]; IF PlaneInTetrahedron[t.plane] THEN { Code: TYPE ~ ARRAY TetFace OF BOOL; SetCodes: PROC RETURNS [ARRAY [0..3) OF Code] ~ { Set: PROC [t: Triple] RETURNS [r: Code] ~ { FOR face: TetFace IN TetFace DO plane: Plane Ź tet.triangles[face].plane; -- planes point inwards r[face] Ź t.x*plane.x+t.y*plane.y+t.z*plane.z+plane.w < 0.0; ENDLOOP; }; RETURN[[Set[t.p1], Set[t.p2], Set[t.p3]]]; }; NoneSet: PROC [c: Code] RETURNS [b: BOOL Ź TRUE] ~ { FOR face: TetFace IN TetFace DO IF c[face] THEN RETURN[FALSE]; ENDLOOP; }; Cross: PROC [p1, p2: Triple] RETURNS [b: BOOL Ź FALSE] ~ { FOR face: TetFace IN TetFace DO IF G3dTriangle.IntersectTriangle[p1, p2, tet.triangles[face]].intersect THEN RETURN[TRUE]; ENDLOOP; }; SetInfo: PROC [c: REF Corner] RETURNS [i: CornerInfo] ~ { d: REAL Ź c.pos.x*t.plane.x+c.pos.y*t.plane.y+c.pos.z*t.plane.z+t.plane.w; RETURN[[TRUE, d, d > 0.0, c]]; }; Hit: PROC [t1, t2: CornerType] RETURNS [b: BOOL] ~ { c1: CornerInfo Ź cInfo[t1]; c2: CornerInfo Ź cInfo[t2]; b Ź G3dTriangle.InsideTriangle[InterpCorners[ }; <> out: ARRAY [0..3) OF Code Ź SetCodes[]; FOR face: TetFace IN TetFace DO IF out[0][face] AND out[1][face] AND out[2][face] THEN RETURN[FALSE]; ENDLOOP; <> FOR i: INT IN [0..3) DO IF NoneSet[out[i]] THEN RETURN[TRUE]; ENDLOOP; <> IF Cross[t.p1, t.p2] OR Cross[t.p2, t.p3] OR Cross[t.p3, t.p1] THEN RETURN[TRUE]; <> -- ab -- IF cInfo[a].pos # cInfo[b].pos AND Hit[a, b] THEN RETURN[TRUE]; -- ac -- IF NOT cInfo[c].set THEN cInfo[c] Ź SetInfo[tet.c]; IF cInfo[a].pos # cInfo[c].pos AND Hit[a, c] THEN RETURN[TRUE]; -- ad -- IF NOT cInfo[d].set THEN cInfo[d] Ź SetInfo[tet.d]; IF cInfo[a].pos # cInfo[d].pos AND Hit[a, d] THEN RETURN[TRUE]; -- bc -- IF cInfo[b].pos # cInfo[c].pos AND Hit[b, c] THEN RETURN[TRUE]; -- bd -- IF cInfo[b].pos # cInfo[d].pos AND Hit[b, d] THEN RETURN[TRUE]; -- cd -- IF cInfo[c].pos # cInfo[d].pos AND Hit[c, d] THEN RETURN[TRUE]; }; RETURN[FALSE]; }; InsideTetrahedron: PUBLIC PROC [t: REF Tetrahedron, p: Triple] RETURNS [BOOL] ~ { Test: PROC [face: TetFace] RETURNS [BOOL] ~ { RETURN[G3dPlane.Side[p, t.triangles[face].plane] = positive]; }; RETURN[Test[abc] AND Test[abd] AND Test[acd] AND Test[bcd]]; }; EdgesFromFace: PUBLIC PROC [f: TetFace] RETURNS [r: ThreeEdges] ~ { r Ź SELECT f FROM abc => [ab, ac, bc], abd => [ab, ad, bd], acd => [ac, ad, cd], ENDCASE => [bc, bd, cd]; }; CornersFromFace: PUBLIC PROC [f: TetFace, t: REF Tetrahedron] RETURNS [r: ThreeCorners] ~ { r Ź SELECT f FROM abc => [t.a, t.b, t.c], abd => [t.a, t.b, t.d], acd => [t.a, t.c, t.d], ENDCASE => [t.b, t.c, t.d]; }; CornersFromEdge: PUBLIC PROC [e: TetEdge, t: REF Tetrahedron] RETURNS [r: TwoCorners] ~ { r Ź SELECT e FROM ab => [t.a, t.b], ac => [t.a, t.c], ad => [t.a, t.d], bc => [t.b, t.c], bd => [t.b, t.d], ENDCASE => [t.c, t.d]; }; MidEdge: PUBLIC PROC [tet: REF Tetrahedron, e: TetEdge] RETURNS [Triple] ~ { tc: TwoCorners Ź CornersFromEdge[e, tet]; RETURN[MidCorner[tc.a, tc.b]]; }; LRFaces: PUBLIC PROC [e: TetEdge] RETURNS [l, r: TetFace] ~ {RETURN[LFace[e], RFace[e]]}; <<>> LFace: PUBLIC PROC [e: TetEdge] RETURNS [f: TetFace] ~ { f Ź SELECT e FROM ab=>abd, ac=>abc, ad=>acd, bc=>bcd, bd=>abd, ENDCASE=>bcd; }; RFace: PUBLIC PROC [e: TetEdge] RETURNS [f: TetFace] ~ { f Ź SELECT e FROM ab=>abc, ac=>acd, ad=>abd, bc=>abc, bd=>bcd, ENDCASE=>acd; }; <> Decompose: PUBLIC PROC [c: REF Cube, action: TetProc, scratch: REF Tetrahedron Ź NIL] ~ { odd: BOOL Ź (c.id.i+c.id.j+c.id.k) MOD 2 # 0; FOR i: INT IN [0..5) DO IF NOT action[InnerGetTetN[c, odd, i, scratch]] THEN EXIT; ENDLOOP; }; <> <> <> << [Artwork node; type 'Artwork on' to command tool] >> GetTetN: PUBLIC PROC [c: REF Cube, which: INT, scratch: REF Tetrahedron Ź NIL] RETURNS [t: REF Tetrahedron] ~ { t Ź InnerGetTetN[c, (c.id.i+c.id.j+c.id.k) MOD 2 # 0, which, scratch]; }; InnerGetTetN: PROC [cube: REF Cube, odd: BOOL, which: INT, scratch: REF Tetrahedron Ź NIL] RETURNS [tet: REF Tetrahedron] ~ { Inner: PROC [a, b, c, d: REF Corner] ~ { <> tet Ź MakeTetrahedron[a, b, c, d, which, scratch]; tet.cube Ź cube; }; IF odd THEN SELECT which FROM << [Artwork node; type 'Artwork on' to command tool] >> 0 => Inner[cube.lbn, cube.lbf, cube.ltn, cube.rbn]; 1 => Inner[cube.ltf, cube.ltn, cube.lbf, cube.rtf]; 2 => Inner[cube.rbf, cube.rbn, cube.rtf, cube.lbf]; 3 => Inner[cube.rtn, cube.rtf, cube.rbn, cube.ltn]; 4 => Inner[cube.lbf, cube.ltn, cube.rbn, cube.rtf]; ENDCASE => ERROR ELSE SELECT which FROM << [Artwork node; type 'Artwork on' to command tool] >> 0 => Inner[cube.lbf, cube.rbf, cube.ltf, cube.lbn]; 1 => Inner[cube.ltn, cube.ltf, cube.rtn, cube.lbn]; 2 => Inner[cube.rbn, cube.rbf, cube.lbn, cube.rtn]; 3 => Inner[cube.rtf, cube.rtn, cube.ltf, cube.rbf]; 4 => Inner[cube.lbn, cube.ltf, cube.rtn, cube.rbf]; ENDCASE => ERROR; }; IDFromPoint: PUBLIC PROC [p: Triple, size: REAL, offset: Triple Ź []] RETURNS [id: IntTriple] ~ { <> <> id Ź [Real.Round[(p.x-offset.x)/size+0.5], Real.Round[(p.y-offset.y)/size+0.5], Real.Round[(p.z-offset.z)/size+0.5]]; }; PointFromID: PUBLIC PROC [id: IntTriple, size: REAL] RETURNS [p: Triple] ~ { p Ź [id.i*size, id.j*size, id.k*size]; }; NewCube: PROC [id: IntTriple, size: REAL, offset: Triple Ź [], scratch: REF Cube Ź NIL] RETURNS [c: REF Cube] ~ { c Ź IF scratch # NIL THEN scratch ELSE NEW[Cube]; c.id Ź id; c.diag Ź size*sqRt3; c.l Ź offset.x+(REAL[id.i]-0.5)*size; c.r Ź c.l+size; c.b Ź offset.y+(REAL[id.j]-0.5)*size; c.t Ź c.b+size; c.n Ź offset.z+(REAL[id.k]-0.5)*size; c.f Ź c.n+size; }; SetCube: PUBLIC PROC [ id: IntTriple, lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf: REF Corner, offset: Triple Ź [], scratch: REF Cube Ź NIL] RETURNS [c: REF Cube] ~ { c Ź NewCube[id, rbn.pos.x-lbn.pos.x, offset, scratch]; c.lbn Ź lbn; c.lbf Ź lbf; c.ltn Ź ltn; c.ltf Ź ltf; c.rbn Ź rbn; c.rbf Ź rbf; c.rtn Ź rtn; c.rtf Ź rtf; }; MakeCube: PUBLIC PROC [id: IntTriple, size: REAL, offset: Triple Ź [], scratch: REF Cube Ź NIL] RETURNS [c: REF Cube] ~ { Set: PROC [x, y, z: REAL, i, j, k: INT] RETURNS [c: REF Corner] ~ { c Ź NEW[Corner Ź [[i, j, k], [x, y, z], 0, unset, FALSE]]; }; c Ź NewCube[id, size, offset, scratch]; c.lbn Ź Set[c.l, c.b, c.n, id.i, id.j, id.k]; c.lbf Ź Set[c.l, c.b, c.f, id.i, id.j, id.k+1]; c.ltn Ź Set[c.l, c.t, c.n, id.i, id.j+1, id.k]; c.ltf Ź Set[c.l, c.t, c.f, id.i, id.j+1, id.k+1]; c.rbn Ź Set[c.r, c.b, c.n, id.i+1, id.j, id.k]; c.rbf Ź Set[c.r, c.b, c.f, id.i+1, id.j, id.k+1]; c.rtn Ź Set[c.r, c.t, c.n, id.i+1, id.j+1, id.k]; c.rtf Ź Set[c.r, c.t, c.f, id.i+1, id.j+1, id.k+1]; }; TriangleInCube: PUBLIC PROC [t: REF Triangle, c: REF Cube] RETURNS [BOOL] ~ { <> PlaneInCube: PROC [p: Plane] RETURNS [in: BOOL] ~ { Pos: PROC [c: REF Corner, type: CornerType] RETURNS [b: BOOL] ~ { d: REAL Ź c.pos.x*p.x+c.pos.y*p.y+c.pos.z*p.z+p.w; cInfo[type] Ź [TRUE, d, b Ź d > 0.0, c]; }; s: BOOL Ź Pos[c.lbn, lbn]; in Ź IF ABS[cInfo[lbn].distance] > c.diag THEN FALSE ELSE Pos[c.lbf, lbf] # s OR Pos[c.ltn, ltn] # s OR Pos[c.ltf, ltf] # s OR Pos[c.rbn, rbn] # s OR Pos[c.rbf, rbf] # s OR Pos[c.rtn, rtn] # s OR Pos[c.rtf, rtf] # s; }; <> CornerInfo: TYPE ~ RECORD [set: BOOL, distance: REAL, pos: BOOL, c: REF Corner]; CornerType: TYPE ~ {lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf}; cInfo: ARRAY CornerType OF CornerInfo Ź ALL[[FALSE, 0.0, FALSE, NIL]]; IF PlaneInCube[t.plane] THEN { << [Artwork node; type 'Artwork on' to command tool] >> Face: TYPE ~ {l, r, b, t, n, f}; Code: TYPE ~ ARRAY Face OF BOOL; SetCodes: PROC RETURNS [ARRAY [0..3) OF Code] ~ { Set: PROC [t: Triple] RETURNS [r: Code] ~ {rŹ[t.xc.r,t.yc.t,t.zc.f]}; RETURN[[Set[t.p1], Set[t.p2], Set[t.p3]]]; }; NoneSet: PROC [c: Code] RETURNS [b: BOOL Ź TRUE] ~ { FOR face: Face IN Face DO IF c[face] THEN RETURN[FALSE]; ENDLOOP; }; SetInfo: PROC [c: REF Corner] RETURNS [i: CornerInfo] ~ { d: REAL Ź c.pos.x*t.plane.x+c.pos.y*t.plane.y+c.pos.z*t.plane.z+t.plane.w; RETURN[[TRUE, d, d > 0.0, c]]; }; Hit: PROC [t1, t2: CornerType] RETURNS [b: BOOL] ~ { c1: CornerInfo Ź cInfo[t1]; c2: CornerInfo Ź cInfo[t2]; b Ź G3dTriangle.InsideTriangle[InterpCorners[ }; out: ARRAY [0..3) OF Code Ź SetCodes[]; <> FOR face: Face IN Face DO IF out[0][face] AND out[1][face] AND out[2][face] THEN RETURN[FALSE]; ENDLOOP; <> FOR i: INT IN [0..3) DO IF NoneSet[out[i]] THEN RETURN[TRUE]; ENDLOOP; <> FOR i: INT IN [0..3) DO Cross: PROC [face: Face, e, e1, e2, a1, a2, amn, amx, b1, b2, bmn, bmx: REAL] RETURNS [b: BOOL Ź FALSE] ~ { IF out[i][face] # out[j][face] THEN { IF amn > amx OR bmn > bmx THEN ERROR; b Ź (a1+ }; }; p, q: Triple; j: INT Ź (i+1) MOD 3; SELECT i FROM 0 => {p Ź t.p1; q Ź t.p2}; 1 => {p Ź t.p2; q Ź t.p3}; ENDCASE => {p Ź t.p3; q Ź t.p1}; <> IF Cross[l, c.l, p.x, q.x, p.y, q.y, c.b, c.t, p.z, q.z, c.n, c.f] OR Cross[r, c.r, p.x, q.x, p.y, q.y, c.b, c.t, p.z, q.z, c.n, c.f] OR Cross[b, c.b, p.y, q.y, p.x, q.x, c.l, c.r, p.z, q.z, c.n, c.f] OR Cross[t, c.t, p.y, q.y, p.x, q.x, c.l, c.r, p.z, q.z, c.n, c.f] OR Cross[n, c.n, p.z, q.z, p.x, q.x, c.l, c.r, p.y, q.y, c.b, c.t] OR Cross[f, c.f, p.z, q.z, p.x, q.x, c.l, c.r, p.y, q.y, c.b, c.t] THEN RETURN[TRUE]; ENDLOOP; <> -- lb -- IF cInfo[lbn].pos # cInfo[lbf].pos AND Hit[lbn, lbf] THEN RETURN[TRUE]; -- ln -- IF NOT cInfo[ltn].set THEN cInfo[ltn] Ź SetInfo[c.ltn]; IF cInfo[lbn].pos # cInfo[ltn].pos AND Hit[lbn, ltn] THEN RETURN[TRUE]; -- lt -- IF NOT cInfo[ltf].set THEN cInfo[ltf] Ź SetInfo[c.ltf]; IF cInfo[ltn].pos # cInfo[ltf].pos AND Hit[ltn, ltf] THEN RETURN[TRUE]; -- lf -- IF cInfo[ltf].pos # cInfo[lbf].pos AND Hit[ltf, lbf] THEN RETURN[TRUE]; -- rb -- IF NOT cInfo[rbn].set THEN cInfo[rbn] Ź SetInfo[c.rbn]; IF NOT cInfo[rbf].set THEN cInfo[rbf] Ź SetInfo[c.rbf]; IF cInfo[rbn].pos # cInfo[rbf].pos AND Hit[rbn, rbf] THEN RETURN[TRUE]; -- rf -- IF NOT cInfo[rtf].set THEN cInfo[rtf] Ź SetInfo[c.rtf]; IF cInfo[rtf].pos # cInfo[rbf].pos AND Hit[rtf, rbf] THEN RETURN[TRUE]; -- rt -- IF NOT cInfo[rtn].set THEN cInfo[rtn] Ź SetInfo[c.rtn]; IF cInfo[rtn].pos # cInfo[rtf].pos AND Hit[rtn, rtf] THEN RETURN[TRUE]; -- rn -- IF cInfo[rtn].pos # cInfo[rbn].pos AND Hit[rtn, rbn] THEN RETURN[TRUE]; -- bn -- IF cInfo[lbn].pos # cInfo[rbn].pos AND Hit[lbn, rbn] THEN RETURN[TRUE]; -- bf -- IF cInfo[lbf].pos # cInfo[rbf].pos AND Hit[lbf, rbf] THEN RETURN[TRUE]; -- tn -- IF cInfo[ltn].pos # cInfo[rtn].pos AND Hit[ltn, rtn] THEN RETURN[TRUE]; -- tf -- IF cInfo[ltf].pos # cInfo[rtf].pos AND Hit[ltf, rtf] THEN RETURN[TRUE]; }; RETURN[FALSE]; }; TrianglesInCube: PUBLIC PROC [ candidates: REF Triangles, cube: REF Cube, scratch: REF Triangles Ź NIL] RETURNS [ret: REF Triangles] ~ { IF (ret Ź scratch) # NIL THEN ret.length Ź 0; FOR i: INT IN [0..candidates.length) DO t: REF Triangle Ź candidates[i]; IF TriangleInCube[t, cube] THEN ret Ź G3dTriangle.AddTriangle[t, ret]; ENDLOOP; }; <> Side: TYPE ~ {in, on, out}; EndState: TYPE ~ RECORD [side: Side, face: TetFace]; EventName: PUBLIC PROC [e: Event] RETURNS [r: Rope.ROPE] ~ { r Ź SELECT e.type FROM entry=> "entry", exit=> "exit", inside=> "inside", ENDCASE=> "graze"; }; scratchEvents: ARRAY [0..10) OF REF Events Ź ALL[NIL]; ObtainEvents: ENTRY PROC RETURNS [REF Events] ~ { FOR i: INT IN [0..10) DO s: REF Events Ź scratchEvents[i]; IF s # NIL THEN {scratchEvents[i] Ź NIL; s.length Ź 0; RETURN[s]}; ENDLOOP; RETURN[NEW[Events[12]]]; }; ReleaseEvents: ENTRY PROC [scratch: REF Events] ~ { FOR i: INT IN [0..10) DO IF scratchEvents[i] = NIL THEN {scratchEvents[i] Ź scratch; EXIT}; ENDLOOP; }; AddEvent: PROC [e: Event, events: REF Events] RETURNS [r: REF Events] ~ { LengthenEvents: PROC [events: REF Events] RETURNS [new: REF Events] ~ { newLength: NAT Ź MAX[Real.Ceiling[1.3*events.maxLength], 3]; new Ź NEW[Events[newLength]]; FOR i: NAT IN [0..events.length) DO new[i] Ź events[i]; ENDLOOP; new.length Ź events.length; }; IF (r Ź events) = NIL THEN r Ź NEW[Events[1]]; IF r.length = r.maxLength THEN r Ź LengthenEvents[r]; r[r.length] Ź e; r.length Ź r.length+1; }; GetEvents: PUBLIC PROC [ tet: REF Tetrahedron, segments: REF Segments, epsilon: REAL, insideToo: BOOL, scratch: REF Events Ź NIL] RETURNS [events: REF Events] ~ { Info: TYPE ~ RECORD [nInts: INT, p1, p2: Triple, f1, f2: TetFace]; GetEndState: PROC [p: Triple] RETURNS [s: EndState Ź [in, abc]] ~ { near: BOOL Ź FALSE; nearFace: TetFace; FOR l: LIST OF TetFace Ź LIST[abc, abd, acd, bcd], l.rest WHILE l # NIL DO d: REAL Ź G3dPlane.DistanceToPoint[p, tet.triangles[l.first].plane, FALSE]; IF d < -epsilon THEN RETURN[[out, abc]]; IF ABS[d] < epsilon THEN {near Ź TRUE; nearFace Ź l.first}; ENDLOOP; RETURN[IF near THEN [on, nearFace] ELSE [in, abc]]; }; end1: EndState Ź GetEndState[segments[0].p1]; IF (events Ź scratch) # NIL THEN events.length Ź 0; FOR s: INT IN [0..segments.length) DO Eq: PROC [p, q: Triple] RETURNS [b: BOOL] ~ {b Ź G3dVector.Distance[p, q] < epsilon}; Add: PROC [type: EventType, p: Triple, face: TetFace Ź abc] ~ { e: Event Ź [p, type, face, s]; IF events # NIL AND events.length > 0 AND Eq[p, events[events.length-1].p] THEN {IF events[events.length-1].type = inside THEN events[events.length-1] Ź e} ELSE events Ź AddEvent[e, events]; }; GetInfo: PROC RETURNS [i: Info] ~ { <> <> <> <> GetAlpha: PROC [f: TetFace] ~ { t: REF Triangle Ź tet.triangles[f]; d1: REAL Ź G3dPlane.DistanceToPoint[seg.p1, t.plane]; d2: REAL Ź G3dPlane.DistanceToPoint[seg.p2, t.plane]; IF (d1 > 0.0) # (d2 > 0.0) THEN { p: Triple Ź G3dVector.Interp[a2 Ź d1/(d1-d2), seg.p1, seg.p2]; IF NOT G3dTriangle.InsideTriangle[p, t] THEN RETURN; IF (i.nInts Ź i.nInts+1) = 1 THEN {a1 Ź a2; i.p1 Ź p; i.f1 Ź f} ELSE IF a2 < a1 THEN {i.p2 Ź i.p1; i.f2 Ź i.f1; i.p1 Ź p; i.f1 Ź f} ELSE {i.p2 Ź p; i.f2 Ź f}; <> IF i.nInts = 2 AND (end1.side # out OR end2.side # out) THEN i.nInts Ź 1; <<(InsideTriangle[] fudge has yielded coincident intersections along tet edge)>> }; }; a1, a2: REAL; FOR f: TetFace IN TetFace DO IF i.nInts = 2 THEN EXIT ELSE GetAlpha[f]; ENDLOOP; }; seg: Segment Ź segments[s]; end2: EndState Ź GetEndState[seg.p2]; i: Info Ź GetInfo[]; SELECT i.nInts FROM 0 => SELECT TRUE FROM end1.side = out AND end2.side = out => NULL; -- completely out end1.side = in AND end2.side = out, end1.side = out AND end2.side = in => ERROR; -- should have nInt > 0 end1.side = in AND end2.side = in, end1.side = on AND end2.side = in, end1.side = in AND end2.side = on, end1.side = on AND end2.side = on => Add[inside, seg.p1, end1.face]; end1.side = on AND end2.side = out => Add[exit, seg.p1, end1.face]; end1.side = out AND end2.side = on => Add[entry, seg.p2, end2.face]; ENDCASE; 1 => SELECT TRUE FROM end1.side = in AND end2.side = in => ERROR; -- no ints, completely in end1.side = out AND end2.side = out => Add[graze, i.p1, i.f1]; end1.side = out AND end2.side = on, end1.side = out AND end2.side = in => Add[entry, i.p1, i.f1]; end1.side = on AND end2.side = in, end1.side = in AND end2.side = on, end1.side = on AND end2.side = on => Add[inside, seg.p1, end1.face]; end1.side = on AND end2.side = out, end1.side = in AND end2.side = out => { Add[inside, seg.p1, end1.face]; Add[exit, i.p1, i.f1]; }; ENDCASE; 2 => SELECT TRUE FROM end1.side = on AND end2.side = on => Add[inside, seg.p1, end1.face]; end1.side = in AND end2.side = in, end1.side = in AND end2.side = out, end1.side = out AND end2.side = in, end1.side = on AND end2.side = in, end1.side = in AND end2.side = on => ERROR; -- couldn't have 2 ints end1.side = on AND end2.side = out => { Add[inside, seg.p1, end1.face]; Add[exit, i.p2, i.f2]; }; end1.side = out AND end2.side = on => Add[entry, i.p1, i.f1]; end1.side = out AND end2.side = out => { Add[entry, i.p1, i.f1]; Add[exit, i.p2, i.f2]; }; ENDCASE; ENDCASE; end1 Ź end2; ENDLOOP; }; BoundaryIntersectTet: PUBLIC PROC [ tet: REF Tetrahedron, segments: REF Segments, insideToo: BOOL Ź TRUE, scratch: REF TetIntersections Ź NIL] RETURNS [REF TetIntersections] ~ { <> << [Artwork node; type 'Artwork on' to command tool] >> <> << [Artwork node; type 'Artwork on' to command tool] >> Next: PROC [i: INT] RETURNS [INT] ~ {RETURN[(i+1) MOD events.length]}; maxD: REAL Ź 0.0; epsilon: REAL Ź 0.01*G3dVector.Distance[tet.a.pos, tet.b.pos]; events: REF Events Ź GetEvents[tet, segments, epsilon, insideToo, ObtainEvents[]]; mark, lastExit, start: INT Ź -1; IF scratch = NIL THEN scratch Ź NEW[TetIntersections]; scratch.length Ź 0; IF events = NIL OR events.length = 0 THEN RETURN[scratch]; <> FOR i: INT IN [0..events.length) DO type: EventType Ź events[i].type; IF type = graze OR type = exit THEN {mark Ź lastExit Ź i; EXIT}; ENDLOOP; IF mark # -1 THEN FOR i: INT Ź Next[mark], Next[i] DO CheckStart: PROC ~ { DistanceBetweenEntries: PROC [e1, e2: Event] RETURNS [d: REAL] ~ { d Ź G3dVector.Distance[e1.p, segments[e1.seg].p2]+ G3dVector.Distance[e2.p, segments[e2.seg].p1]; FOR s: INT IN (e1.seg..e2.seg) DO d Ź d+segments[s].length; ENDLOOP; }; d: REAL Ź DistanceBetweenEntries[events[lastExit], events[i]]; IF d > maxD THEN {maxD Ź d; start Ź i}; }; SELECT events[i].type FROM graze => {CheckStart[]; lastExit Ź i}; entry => CheckStart[]; exit => lastExit Ź i; ENDCASE; IF i = mark THEN EXIT; ENDLOOP; <> IF start # -1 THEN FOR i: INT IN [0..events.length) DO e: Event Ź events[start]; scratch Ź AddTetIntersection[scratch, [e.p, e.type = inside, e.face]]; start Ź Next[start]; ENDLOOP; ReleaseEvents[events]; RETURN[scratch]; }; AddTetIntersection: PROC [seq: REF TetIntersections, int: TetIntersection] RETURNS [REF TetIntersections] ~ { IF seq = NIL THEN seq Ź NEW[TetIntersections[3]]; IF seq.length = seq.maxLength THEN { old: REF TetIntersections Ź seq; seq Ź NEW[TetIntersections[MAX[Real.Ceiling[1.3*old.maxLength], 3]]]; FOR i: INT IN [0..old.length) DO seq[i] Ź old[i]; ENDLOOP; seq.length Ź old.length; }; seq[seq.length] Ź int; seq.length Ź seq.length+1; RETURN[seq]; }; <> <> <> <> <> <> <> <> <<~ {>> <> <> <> <> <> <> <> <> <> < 0.0) # (d2 > 0.0) THEN {>> <> <> <> <> <> <> <> <> <<}>> <> <> <> <> <> <> <> <> <<};>> <<};>> <<};>> <> <> <> <<1 => {>> <> <> <> <> <> <> <> <> <> <<[] Ź action[info1.p, graze, info1.f];>> <> <<};>> <> <> <> <start AND Eq[p,segments[i-1].p2])>> <> <> <> <<[] Ź action[infoStop.p, exit, infoStop.f];>> <<};>> <<2 => IF action[info1.p, entry, info1.f] THEN [] Ź action[info2.p, exit, info2.f];>> < LOOP;>> <> <> <> <<};>> <<>> TetCenter: PUBLIC PROC [tet: REF Tetrahedron] RETURNS [t: Triple] ~ { t.x Ź 0.25*(tet.a.pos.x+tet.b.pos.x+tet.c.pos.x+tet.d.pos.x); t.y Ź 0.25*(tet.a.pos.y+tet.b.pos.y+tet.c.pos.y+tet.d.pos.y); t.z Ź 0.25*(tet.a.pos.z+tet.b.pos.z+tet.c.pos.z+tet.d.pos.z); }; <<>> 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]; }; <<>> <> <> < size/2>> <> <> <<}>> <> <> <> <> <> <> <> <> <<};>> <> <> <> <> <> <> <> <> <> <<}>> <> <<};>> <<};>> <<>> scratchTetIntersections: ARRAY [0..10) OF REF TetIntersections Ź ALL[NIL]; ObtainTetIntersections: PUBLIC ENTRY PROC RETURNS [REF TetIntersections] ~ { FOR i: INT IN [0..10) DO s: REF TetIntersections Ź scratchTetIntersections[i]; IF s # NIL THEN {scratchTetIntersections[i] Ź NIL; s.length Ź 0; RETURN[s]}; ENDLOOP; RETURN[NEW[TetIntersections[12]]]; }; ReleaseTetIntersections: PUBLIC ENTRY PROC [scratch: REF TetIntersections] ~ { FOR i: INT IN [0..10) DO IF scratchTetIntersections[i] = NIL THEN {scratchTetIntersections[i] Ź scratch; EXIT}; ENDLOOP; }; END.