<<>> <> <> <> DIRECTORY Basics, G3dBasic, G3dShape, G3dVector, ImplicitDefs, ImplicitMinimal, Random, Real, Rope; ImplicitTetImpl: CEDAR PROGRAM IMPORTS Basics, G3dBasic, G3dShape, G3dVector, Random, Real, Rope EXPORTS ImplicitMinimal ~ BEGIN <> Triple: TYPE ~ G3dBasic.Triple; NatSequence: TYPE ~ G3dBasic.NatSequence; ValueProc: TYPE ~ ImplicitDefs.ValueProc; Vertex: TYPE ~ ImplicitTet.Vertex; Vertices: TYPE ~ ImplicitTet.Vertices; ROPE: TYPE ~ Rope.ROPE; Cube: TYPE ~ RECORD [i, j, k: INTEGER, lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf: Corner]; Corner: TYPE ~ RECORD [i, j, k: INTEGER ¬ 0, x, y, z, value: REAL ¬ 0]; Center: TYPE ~ RECORD [i, j, k: INTEGER ¬ 0]; Centers: TYPE ~ RECORD [s: SEQUENCE size: INT OF LIST OF Center]; Edge: TYPE ~ RECORD [i1, j1, k1, i2, j2, k2, id: INTEGER ¬ 0]; Edges: TYPE ~ RECORD [s: SEQUENCE size: INT OF LIST OF Edge]; <<>> <> <> <> Abort: ERROR = CODE; Polygonize: PUBLIC PROC [ function: ValueProc, level, size, delta: REAL, start: Triple ¬ [], clientData: REF ANY ¬ NIL, triangleAction: ImplicitTet.TriangleProc] RETURNS [error: ROPE ¬ NIL] ~ { ENABLE Abort => {error ¬ "aborted"; GOTO Failed}; Poly: PROC [cube: Cube] ~ { VertexID: PROC [c1, c2: Corner] RETURNS [id: INTEGER] ~ { <> v: Vertex; IF (c1.value < 0.0) = (c2.value < 0.0) THEN RETURN[-1]; -- no vertex IF (id ¬ GetEdge[edges, c1.i, c1.j, c1.k, c2.i, c2.j, c2.k]) # -1 THEN RETURN[id]; v.position ¬ Converge[[c1.x, c1.y, c1.z], [c2.x, c2.y, c2.z], c1.value, level, function, clientData]; v.normal ¬ Normal[v.position, function, clientData, delta]; vertices ¬ AddToVertices[vertices, v]; SetEdge[edges, c1.i, c1.j, c1.k, c2.i, c2.j, c2.k, id ¬ vertices.length-1]; }; TriangulateTetrahedron: PROC [a, b, c, d: Corner] ~ { <> Output: PROC [i0, i1, i2: INTEGER] ~ INLINE { IF NOT triangleAction[i0, i1, i2, vertices] THEN ERROR Abort; }; aPos, bPos, cPos, dPos: BOOL; index, e1, e2, e3, e4, e5, e6: INT ¬ 0; IF (aPos ¬ a.value > 0) THEN index ¬ index+8; IF (bPos ¬ b.value > 0) THEN index ¬ index+4; IF (cPos ¬ c.value > 0) THEN index ¬ index+2; IF (dPos ¬ d.value > 0) THEN index ¬ index+1; IF aPos # bPos THEN e1 ¬ VertexID[a, b]; IF aPos # cPos THEN e2 ¬ VertexID[a, c]; IF aPos # dPos THEN e3 ¬ VertexID[a, d]; IF bPos # cPos THEN e4 ¬ VertexID[b, c]; IF bPos # dPos THEN e5 ¬ VertexID[b, d]; IF cPos # dPos THEN e6 ¬ VertexID[c, d]; SELECT index FROM 1 => Output[e5, e6, e3]; 2 => Output[e2, e6, e4]; 3 => {Output[e3, e5, e4]; Output[e3, e4, e2]}; 4 => Output[e1, e4, e5]; 5 => {Output[e3, e1, e4]; Output[e3, e4, e6]}; 6 => {Output[e1, e2, e6]; Output[e1, e6, e5]}; 7 => Output[e1, e2, e3]; 8 => Output[e1, e3, e2]; 9 => {Output[e1, e5, e6]; Output[e1, e6, e2]}; 10 => {Output[e1, e3, e6]; Output[e1, e6, e4]}; 11 => Output[e1, e5, e4]; 12 => {Output[e3, e2, e4]; Output[e3, e4, e5]}; 13 => Output[e6, e2, e4]; 14 => Output[e5, e3, e6]; ENDCASE; -- 0, 15 }; IF (cube.i+cube.j+cube.k) MOD 2 # 0 THEN { TriangulateTetrahedron[cube.lbn, cube.lbf, cube.ltn, cube.rbn]; TriangulateTetrahedron[cube.ltf, cube.ltn, cube.lbf, cube.rtf]; TriangulateTetrahedron[cube.rbf, cube.rbn, cube.rtf, cube.lbf]; TriangulateTetrahedron[cube.rtn, cube.rtf, cube.rbn, cube.ltn]; TriangulateTetrahedron[cube.lbf, cube.ltn, cube.rbn, cube.rtf]; } ELSE { TriangulateTetrahedron[cube.lbf, cube.rbf, cube.ltf, cube.lbn]; TriangulateTetrahedron[cube.ltn, cube.ltf, cube.rtn, cube.lbn]; TriangulateTetrahedron[cube.rbn, cube.rbf, cube.lbn, cube.rtn]; TriangulateTetrahedron[cube.rtf, cube.rtn, cube.ltf, cube.rbf]; TriangulateTetrahedron[cube.lbn, cube.ltf, cube.rtn, cube.rbf]; }; }; SetCorner: PROC [x, y, z: REAL, i, j, k: INTEGER] RETURNS [Corner] ~ { RETURN[[i, j, k, x, y, z, function[[x, y, z], clientData]-level]]; }; TestFace: PROC [i, j, k: INTEGER, c1, c2, c3, c4: Corner, face: {l, r, b, t, n, f}] ~ { cube: Cube; sign: BOOL ¬ c1.value > 0; lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf: Corner; IF (c2.value > 0) = sign AND (c3.value > 0) = sign AND (c4.value > 0) = sign THEN RETURN; IF SetCenter[centers, i, j, k] THEN RETURN; SELECT face FROM l => { -- given right face (rbn, rbf, rtn, rtf); compute left (lbn, lbf, ltn, ltf) x: REAL ¬ c1.x-size; i: INTEGER ¬ c1.i-1; rbn ¬ c1; rbf ¬ c2; rtn ¬ c3; rtf ¬ c4; lbn ¬ SetCorner[x, rbn.y, rbn.z, i, rbn.j, rtn.k]; lbf ¬ SetCorner[x, rbf.y, rbf.z, i, rbf.j, rbf.k]; ltn ¬ SetCorner[x, rtn.y, rtn.z, i, rtn.j, rtn.k]; ltf ¬ SetCorner[x, rtf.y, rtf.z, i, rtf.j, rtf.k]; }; r => { -- given left face (lbn, lbf, ltn, ltf); compute right (rbn, rbf, rtn, rtf) x: REAL ¬ c1.x+size; i: INTEGER ¬ c1.i+1; lbn ¬ c1; lbf ¬ c2; ltn ¬ c3; ltf ¬ c4; rbn ¬ SetCorner[x, lbn.y, lbn.z, i, lbn.j, lbn.k]; rbf ¬ SetCorner[x, lbf.y, lbf.z, i, lbf.j, lbf.k]; rtn ¬ SetCorner[x, ltn.y, ltn.z, i, ltn.j, ltn.k]; rtf ¬ SetCorner[x, ltf.y, ltf.z, i, ltf.j, ltf.k]; }; b => { -- given top face (ltn, ltf, rtn, rtf); compute bottom (lbn, lbf, rbn, rbf) face y: REAL ¬ c1.y-size; j: INTEGER ¬ c1.j-1; ltn ¬ c1; ltf ¬ c2; rtn ¬ c3; rtf ¬ c4; lbn ¬ SetCorner[ltn.x, y, ltn.z, ltn.i, j, ltn.k]; lbf ¬ SetCorner[ltf.x, y, ltf.z, ltf.i, j, ltf.k]; rbn ¬ SetCorner[rtn.x, y, rtn.z, rtn.i, j, rtn.k]; rbf ¬ SetCorner[rtf.x, y, rtf.z, rtf.i, j, rtf.k]; }; t => { -- given bottom face (lbn, lbf, rbn, rbf); compute top (ltn, ltf, rtn, rtf) y: REAL ¬ c1.y+size; j: INTEGER ¬ c1.j+1; lbn ¬ c1; lbf ¬ c2; rbn ¬ c3; rbf ¬ c4; ltn ¬ SetCorner[lbn.x, y, lbn.z, lbn.i, j, lbn.k]; ltf ¬ SetCorner[lbf.x, y, lbf.z, lbf.i, j, lbf.k]; rtn ¬ SetCorner[rbn.x, y, rbn.z, rbn.i, j, rbn.k]; rtf ¬ SetCorner[rbf.x, y, rbf.z, rbf.i, j, rbf.k]; }; n => { -- given far face (lbf, ltf, rbf, rtf); compute near (lbn, ltn, rbn, rtn) z: REAL ¬ c1.z-size; k: INTEGER ¬ c1.k-1; lbf ¬ c1; ltf ¬ c2; rbf ¬ c3; rtf ¬ c4; lbn ¬ SetCorner[lbf.x, lbf.y, z, lbf.i, lbf.j, k]; ltn ¬ SetCorner[ltf.x, ltf.y, z, ltf.i, ltf.j, k]; rbn ¬ SetCorner[rbf.x, rbf.y, z, rbf.i, rbf.j, k]; rtn ¬ SetCorner[rtf.x, rtf.y, z, rtf.i, rtf.j, k]; }; f => { -- given near face (lbn, ltn, rbn, rtn); compute far (lbf, ltf, rbf, rtf) z: REAL ¬ c1.z+size; k: INTEGER ¬ c1.k+1; lbn ¬ c1; ltn ¬ c2; rbn ¬ c3; rtn ¬ c4; lbf ¬ SetCorner[lbn.x, lbn.y, z, lbn.i, lbn.j, k]; ltf ¬ SetCorner[ltn.x, ltn.y, z, ltn.i, ltn.j, k]; rbf ¬ SetCorner[rbn.x, rbn.y, z, rbn.i, rbn.j, k]; rtf ¬ SetCorner[rtn.x, rtn.y, z, rtn.i, rtn.j, k]; }; ENDCASE; cube ¬ [i, j, k, lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf]; cubes ¬ IF cubes = NIL THEN LIST[cube] ELSE CONS[cube, cubes]; }; Propagate: PROC [c: Cube] ~ { <> TestFace[c.i-1, c.j, c.k, c.lbn, c.lbf, c.ltn, c.ltf, l]; TestFace[c.i+1, c.j, c.k, c.rbn, c.rbf, c.rtn, c.rtf, r]; TestFace[c.i, c.j-1, c.k, c.lbn, c.lbf, c.rbn, c.rbf, b]; TestFace[c.i, c.j+1, c.k, c.ltn, c.ltf, c.rtn, c.rtf, t]; TestFace[c.i, c.j, c.k-1, c.lbn, c.ltn, c.rbn, c.rtn, n]; TestFace[c.i, c.j, c.k+1, c.lbf, c.ltf, c.rbf, c.rtf, f]; }; cubes: LIST OF Cube; vertices: Vertices ¬ NIL; lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf: Corner; edges: REF Edges ¬ NEW[Edges[2*hashSize]]; centers: REF Centers ¬ NEW[Centers[hashSize]]; cen: Triple ¬ FindStart[start, function, clientData, size, level ! Error => {error ¬ reason; GOTO Failed}]; FOR i: INT IN [0..centers.size) DO centers[i] ¬ NIL; ENDLOOP; FOR i: INT IN [0..edges.size) DO edges[i] ¬ NIL; ENDLOOP; lbn ¬ SetCorner[cen.x-0.5*size, cen.y-0.5*size, cen.z-0.5*size, 0, 0, 0]; lbf ¬ SetCorner[cen.x-0.5*size, cen.y-0.5*size, cen.z+0.5*size, 0, 0, 1]; ltn ¬ SetCorner[cen.x-0.5*size, cen.y+0.5*size, cen.z-0.5*size, 0, 1, 0]; ltf ¬ SetCorner[cen.x-0.5*size, cen.y+0.5*size, cen.z+0.5*size, 0, 1, 1]; rbn ¬ SetCorner[cen.x+0.5*size, cen.y-0.5*size, cen.z-0.5*size, 1, 0, 0]; rbf ¬ SetCorner[cen.x+0.5*size, cen.y-0.5*size, cen.z+0.5*size, 1, 0, 1]; rtn ¬ SetCorner[cen.x+0.5*size, cen.y+0.5*size, cen.z-0.5*size, 1, 1, 0]; rtf ¬ SetCorner[cen.x+0.5*size, cen.y+0.5*size, cen.z+0.5*size, 1, 1, 1]; Propagate[[0, 0, 0, lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf]]; WHILE cubes # NIL DO c: Cube ¬ cubes.first; Poly[c]; cubes ¬ cubes.rest; Propagate[c]; ENDLOOP; EXITS Failed => RETURN[error]; }; <> hashNBits: INTEGER ~ 5; hashMask: INTEGER ~ Basics.BITSHIFT[1, hashNBits]-1; hashSize: INTEGER ~ Basics.BITSHIFT[1, 3*hashNBits]; <<>> HashAnd: PROC [i: INTEGER] RETURNS [j: WORD] ~INLINE {j¬Basics.BITAND[ABS[i],hashMask]}; HashShift: PROC [i: INTEGER] RETURNS [j: WORD] ~ INLINE {j¬Basics.BITSHIFT[i, hashNBits]}; Hash: PROC [i, j, k: INTEGER] RETURNS [id: INTEGER] ~ { <> id ¬ Basics.BITOR[ HashShift[ Basics.BITOR[ HashShift[HashAnd[i]], HashAnd[j]]], HashAnd[k]]; }; SetCenter: PROC [table: REF Centers, i, j, k: INTEGER] RETURNS [b: BOOL ¬ FALSE] ~ { <> index: INTEGER ¬ Hash[i, j, k]; list: LIST OF Center ¬ table[index]; FOR l: LIST OF Center ¬ list, l.rest WHILE l # NIL DO IF l.first.i = i AND l.first.j = j AND l.first.k = k THEN RETURN[TRUE]; ENDLOOP; table[index] ¬ CONS[[i, j, k], list]; }; Order: PROC [i1, j1, k1, i2, j2, k2: INTEGER] RETURNS [a, b, c, d, e, f: INTEGER] ~ INLINE { IF i1 > i2 OR (i1 = i2 AND (j1 > j2 OR (j1 = j2 AND k1 > k2))) THEN RETURN[i2, j2, k2, i1, j1, k1] ELSE RETURN[i1, j1, k1, i2, j2, k2]; }; SetEdge: PROC [table: REF Edges, i1, j1, k1, i2, j2, k2, id: INTEGER] ~ { index: INTEGER; [i1, j1, k1, i2, j2, k2] ¬ Order[i1, j1, k1, i2, j2, k2]; index ¬ Hash[i1, j1, k1]+Hash[i2, j2, k2]; table[index] ¬ CONS[[i1, j1, k1, i2, j2, k2, id], table[index]]; }; GetEdge: PROC [table: REF Edges, i1, j1, k1, i2, j2, k2: INTEGER] RETURNS [i: INTEGER ¬ -1] ~ { [i1, j1, k1, i2, j2, k2] ¬ Order[i1, j1, k1, i2, j2, k2]; FOR l: LIST OF Edge ¬ table[Hash[i1, j1, k1]+Hash[i2, j2, k2]], l.rest WHILE l # NIL DO IF l.first.i1 = i1 AND l.first.j1 = j1 AND l.first.k1 = k1 AND l.first.i2 = i2 AND l.first.j2 = j2 AND l.first.k2 = k2 THEN RETURN[l.first.id]; ENDLOOP; }; <> AddToVertices: PROC [vertices: Vertices, v: Vertex] RETURNS [new: Vertices] ~ { new ¬ IF vertices = NIL THEN NEW[ImplicitTet.VerticesRep[1]] ELSE vertices; IF new.length = new.maxLength THEN { old: Vertices ¬ new; new ¬ NEW[ImplicitTet.VerticesRep[MAX[Real.Ceiling[1.5*old.maxLength], 3]]]; FOR i: NAT IN [0..new.length ¬ old.length) DO new[i] ¬ old[i]; ENDLOOP; }; new[new.length] ¬ v; new.length ¬ new.length+1; }; Error: PUBLIC ERROR [reason: ROPE] = CODE; ShapeFromFunction: PUBLIC PROC [ function: ValueProc, level, size, delta: REAL, start: Triple ¬ [], clientData: REF ANY ¬ NIL, status: ImplicitTet.StatusProc ¬ NIL] RETURNS [s: G3dShape.Shape ¬ NIL] ~ { Tri: ImplicitTet.TriangleProc ~ { triangle: NatSequence ¬ NEW[G3dBasic.NatSequenceRep[3]]; triangle[0] ¬ i0; triangle[1] ¬ i1; triangle[2] ¬ i2; triangle.length ¬ 3; points ¬ vertices; polys ¬ G3dBasic.AddToSurfaceSequence[polys, [NIL, triangle]]; IF status # NIL THEN RETURN[status[i0, i1, i2, vertices]]; }; points: ImplicitTet.Vertices ¬ NIL; polys: G3dBasic.SurfaceSequence ¬ NIL; error: ROPE ¬ Polygonize[function, level, size, delta, start, clientData, Tri]; IF points = NIL OR points.length = 0 THEN error ¬ Rope.Cat[error, IF error # NIL THEN ", " ELSE NIL, "no vertices"]; IF polys = NIL OR polys.length = 0 THEN error ¬ Rope.Cat[error, IF error # NIL THEN ", " ELSE NIL, "no polygons"]; IF error # NIL THEN Error[error] ELSE { vertices: G3dBasic.TripleSequence ¬ NEW[G3dBasic.TripleSequenceRep[points.length]]; normals: G3dBasic.TripleSequence ¬ NEW[G3dBasic.TripleSequenceRep[points.length]]; FOR n: NAT IN [0..vertices.length ¬ normals.length ¬ points.length) DO vertices[n] ¬ points[n].position; normals[n] ¬ points[n].normal; ENDLOOP; s ¬ G3dShape.ShapeFromData[NIL, polys, vertices, normals,,,, FALSE,, $Triangulated]; s.selected ¬ TRUE; }; }; <> FindStart: PUBLIC PROC [ ballpark: Triple, function: ValueProc, clientData: REF ANY, size, level: REAL] RETURNS [c: Triple] ~ { Found: TYPE ~ RECORD [bad: BOOL, p: Triple]; Find: PROC [sign: {negative, positive}] RETURNS [found: Found ¬ [FALSE, ballpark]] ~ { value: REAL; FOR i: NAT IN [0..5000) DO found.p.x ¬ found.p.x+size*(Random.NextInt[rs]/REAL[LAST[INT]]-0.5); found.p.y ¬ found.p.y+size*(Random.NextInt[rs]/REAL[LAST[INT]]-0.5); found.p.z ¬ found.p.z+size*(Random.NextInt[rs]/REAL[LAST[INT]]-0.5); value ¬ function[found.p, clientData]-level; IF (sign = positive AND value > 0) OR (sign = negative AND value < 0) THEN RETURN; size ¬ size*1.005; ENDLOOP; found.bad ¬ TRUE; }; rs: Random.RandomStream ¬ Random.Create[]; pos: Found ¬ Find[positive]; neg: Found ¬ Find[negative]; IF pos.bad OR neg.bad THEN Error[Rope.Cat["can't find (", IF pos.bad THEN "+" ELSE "-", ") starting point"]] ELSE c ¬ Converge[pos.p, neg.p, 1.0, level, function, clientData]; }; Converge: PROC [pos, neg: Triple, val, level: REAL, function: ValueProc, clientData: REF ANY] RETURNS [p: Triple] ~ { nTries: NAT ¬ 0; p ¬ pos; IF val < 0 THEN {c: Triple ¬ neg; neg ¬ pos; pos ¬ c}; DO p ¬ G3dVector.Midpoint[pos, neg]; IF (nTries ¬ nTries+1) = 10 THEN EXIT; IF (function[p, clientData]-level) > 0.0 THEN pos ¬ p ELSE neg ¬ p; ENDLOOP; }; Normal: PROC [point: Triple, function: ValueProc, clientData: REF ANY, delta: REAL] RETURNS [n: Triple] ~ { <> v: REAL ¬ function[point, clientData]; x: REAL ¬ v-function[[point.x+delta, point.y, point.z], clientData]; y: REAL ¬ v-function[[point.x, point.y+delta, point.z], clientData]; z: REAL ¬ v-function[[point.x, point.y, point.z+delta], clientData]; n ¬ G3dVector.Unit[[x, y, z]]; }; END. .. <> <> <> < i2 OR (i1 = i2 AND (j1 > j2 OR (j1 = j2 AND k1 > k2))) THEN {>> <> <<};>> <> <<[Basics.HFromInt16[Hash3d[i1, j1, k1]], Basics.HFromInt16[Hash3d[i2, j2, k2]]]];>> <<};>> <<>>