DIRECTORY Basics, G3dBasic, G3dShape, G3dVector, ImplicitDefs, ImplicitMinimal, Random, Real, Rope;

ImplicitMinimalImpl: 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 ~ ImplicitMinimal.Vertex;
Vertices:		TYPE ~ ImplicitMinimal.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];
Face:			TYPE ~ {l, r, b, t, n, f, unknown};
Edge:			TYPE ~ {lb, lt, ln, lf, rb, rt, rn, rf, bn, bf, tn, tf};
Entry:		TYPE ~ RECORD [i, j, k, id: INTEGER];
Table:			TYPE ~ RECORD [s: SEQUENCE length: CARDINAL OF LIST OF Entry];
Abort: ERROR = CODE;

Polygonize: PUBLIC PROC [
function: ValueProc,
start: Triple ฌ [],
level: REAL,
size: REAL,
delta: REAL,
clientData: REF ANY ฌ NIL,
polygonAction: ImplicitMinimal.PolygonProc ฌ NIL,
triangleAction: ImplicitMinimal.TriangleProc ฌ NIL,
adaptive: REAL ฌ 0.0]
RETURNS [error: ROPE ฌ NIL]
~ {
ENABLE Abort => {error ฌ "aborted"; GOTO Failed};
Poly: PROC [c: Cube] ~ {
DirectedFace: PROC [edge: Edge] RETURNS [Face] ~ {
RETURN[SELECT edge FROM
lb			=> IF c.lbn.value < 0.0 THEN l	ELSE b,
lt			=> IF c.ltn.value	< 0.0 THEN t	ELSE l,
ln			=> IF c.lbn.value < 0.0 THEN n ELSE l,
lf			=> IF c.lbf.value	< 0.0 THEN l	ELSE f,
rb			=> IF c.rbn.value < 0.0 THEN b ELSE r,
rt			=> IF c.rtn.value	< 0.0 THEN r	ELSE t,
rn			=> IF c.rbn.value < 0.0 THEN r	ELSE n,
rf			=> IF c.rbf.value	< 0.0 THEN f	ELSE r,
bn			=> IF c.lbn.value < 0.0 THEN b	ELSE n,
bf			=> IF c.lbf.value	< 0.0 THEN f	ELSE b,
tn			=> IF c.ltn.value	< 0.0 THEN n	ELSE t,
ENDCASE	=> IF c.ltf.value	< 0.0 THEN t	ELSE f];
};
VertexID: PROC [e: Edge] RETURNS [id: INTEGER] ~ {
v: Vertex;
table: REF Table;
dir: Face; -- from p1 to p2; e.g., if p1 = lbn and p2 = lbf, direction = f
p1, p2: Corner;
SELECT e FROM
lb	=> {table ฌ zEdges;	dir ฌ f; p1 ฌ c.lbn;	p2 ฌ c.lbf};
lt	=> {table ฌ zEdges;	dir ฌ f; p1 ฌ c.ltn;	p2 ฌ c.ltf};
ln	=> {table ฌ yEdges;	dir ฌ t; p1 ฌ c.lbn;	p2 ฌ c.ltn};
lf	=> {table ฌ yEdges;	dir ฌ t; p1 ฌ c.lbf;	p2 ฌ c.ltf};
rb	=> {table ฌ zEdges;	dir ฌ f; p1 ฌ c.rbn;	p2 ฌ c.rbf};
rt	=> {table ฌ zEdges;	dir ฌ f; p1 ฌ c.rtn;	p2 ฌ c.rtf};
rn	=> {table ฌ yEdges;	dir ฌ t; p1 ฌ c.rbn;	p2 ฌ c.rtn};
rf	=> {table ฌ yEdges;	dir ฌ t; p1 ฌ c.rbf;	p2 ฌ c.rtf};
bn	=> {table ฌ xEdges;	dir ฌ r; p1 ฌ c.lbn;	p2 ฌ c.rbn};
bf	=> {table ฌ xEdges;	dir ฌ r; p1 ฌ c.lbf;	p2 ฌ c.rbf};
tn	=> {table ฌ xEdges;	dir ฌ r; p1 ฌ c.ltn;	p2 ฌ c.rtn};
tf	=> {table ฌ xEdges;	dir ฌ r; p1 ฌ c.ltf;	p2 ฌ c.rtf};
ENDCASE;
IF (p1.value < 0.0) = (p2.value < 0.0) THEN RETURN[-1]; -- no vertex
IF (id ฌ GetIJK[table, p1.i, p1.j, p1.k]) # -1 THEN RETURN[id];
v.position ฌ Converge[[p1.x, p1.y, p1.z], [p2.x, p2.y, p2.z], p1.value, level, function, clientData, dir];
v.normal ฌ Normal[v.position, function, clientData, delta];
vertices ฌ AddToVertices[vertices, v];
[] ฌ SetIJK[table, p1.i, p1.j, p1.k, id ฌ vertices.length-1];
};
AddPolyCenter: PROC [level, nVertices: INTEGER, ids: ARRAY[0..6) OF INTEGER] ~ {
NewVertex: PROC [p: Triple] ~ {
vertices ฌ AddToVertices[vertices, [p, Normal[p, function, clientData, delta]]];
};
Opp: PROC [a,b: REAL] RETURNS [c: BOOL]~{cฌ(a<0 AND b>0) OR (a>0 AND b<0)};
Add: PROC [a,b: Triple] RETURNS [c: Triple] ~ {cฌ[a.x+b.x, a.y+b.y, a.z+b.z]};
Mul: PROC [a: Triple, s: REAL] RETURNS [c: Triple] ~ {cฌ[s*a.x, s*a.y, s*a.z]};
val: REAL;
p0, p1, normal: Triple ฌ [0, 0, 0];
vertexID: INTEGER ฌ vertices.length;
tris: ARRAY [0..6) OF INTEGER;
FOR i: NAT IN [0..nVertices) DO p0 ฌ Add[p0, vertices[ids[i]].position]; ENDLOOP;
val ฌ function[p0 ฌ Mul[p0, 1.0/REAL[nVertices]], clientData];
normal ฌ Normal[p0, function, clientData, delta];
FOR i: NAT IN [0..10) DO
IF Opp[val, function[p1 ฌ Add[p0, Mul[normal, step*i]], clientData]] OR
Opp[val, function[p1 ฌ Add[p0, Mul[normal, -step*i]], clientData]] THEN {
NewVertex[Converge[p0, p1, val, level, function, clientData]];
EXIT;
};
REPEAT FINISHED => NewVertex[p0];
ENDLOOP;
tris[0] ฌ vertexID;
IF triangleAction # NIL THEN
FOR i: NAT IN [0..nVertices) DO
IF NOT triangleAction[vertexID, ids[i], ids[(i+1) MOD nVertices], vertices]
THEN ERROR Abort;
IF level = 0 THEN {
tris[1] ฌ ids[i];
tris[2] ฌ ids[(i+1) MOD nVertices];
IF Subdivide[3, tris] THEN AddPolyCenter[1, 3, tris];
};
ENDLOOP;
IF polygonAction # NIL THEN {
FOR i: NAT IN [0..nVertices) DO
tris[1] ฌ ids[i];
tris[2] ฌ ids[(i+1) MOD nVertices];
IF NOT polygonAction[[3, tris], vertices] THEN ERROR Abort;
IF level = 0 AND Subdivide[3, tris] THEN AddPolyCenter[1, 3, tris];
ENDLOOP;
};
};
Subdivide: PROC [nVertices: INTEGER, ids: ARRAY [0..6) OF INTEGER] RETURNS [BOOL] ~ {
aveNormal: Triple ฌ [0., 0., 0.];
FOR i: NAT IN [0..nVertices) DO
aveNormal ฌ G3dVector.Add[aveNormal, vertices[ids[i]].normal];
ENDLOOP;
aveNormal ฌ G3dVector.Unit[aveNormal];
FOR i: NAT IN [0..nVertices) DO
IF G3dVector.Dot[aveNormal, vertices[ids[i]].normal] < adaptive THEN RETURN[TRUE];
ENDLOOP;
RETURN[FALSE];
};
id, nVertices: INTEGER;
ids: ARRAY [0..6) OF INTEGER;
done: ARRAY Edge OF BOOL ฌ ALL[FALSE];
FOR edge: Edge IN Edge DO
IF NOT done[edge] AND (ids[0] ฌ VertexID[edge]) # -1 THEN {
e, start: Edge ฌ edge;
face: Face ฌ DirectedFace[e];
nVertices ฌ 1;
DO
done[e ฌ NextCWEdge[e, face]] ฌ TRUE;
IF e = start THEN EXIT;
IF (id ฌ VertexID[e]) # -1 THEN {
IF nVertices = 6 THEN RETURN; -- ERROR
ids[nVertices] ฌ id;
nVertices ฌ nVertices+1;
face ฌ OtherFace[e, face];
};
ENDLOOP;
IF nVertices < 3 THEN RETURN; -- ERROR
IF triangleAction # NIL THEN AddPolyCenter[0, nVertices, ids];
IF polygonAction # NIL THEN {
IF adaptive # 0.0 AND triangleAction = NIL AND Subdivide[nVertices, ids]
THEN AddPolyCenter[0, nVertices, ids]
ELSE IF NOT polygonAction[[nVertices, ids], vertices] THEN ERROR Abort;
};
};
ENDLOOP;
};
Set: 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: Face] ~ {
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 SetIJK[centers, i, j, k, 1] 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 ฌ Set[x, rbn.y, rbn.z, i, rbn.j, rtn.k];
lbf ฌ Set[x, rbf.y, rbf.z, i, rbf.j, rbf.k];
ltn ฌ Set[x, rtn.y, rtn.z, i, rtn.j, rtn.k];
ltf ฌ Set[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 ฌ Set[x, lbn.y, lbn.z, i, lbn.j, lbn.k];
rbf ฌ Set[x, lbf.y, lbf.z, i, lbf.j, lbf.k];
rtn ฌ Set[x, ltn.y, ltn.z, i, ltn.j, ltn.k];
rtf ฌ Set[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 ฌ Set[ltn.x, y, ltn.z, ltn.i, j, ltn.k];
lbf ฌ Set[ltf.x, y, ltf.z, ltf.i, j, ltf.k];
rbn ฌ Set[rtn.x, y, rtn.z, rtn.i, j, rtn.k];
rbf ฌ Set[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 ฌ Set[lbn.x, y, lbn.z, lbn.i, j, lbn.k];
ltf ฌ Set[lbf.x, y, lbf.z, lbf.i, j, lbf.k];
rtn ฌ Set[rbn.x, y, rbn.z, rbn.i, j, rbn.k];
rtf ฌ Set[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 ฌ Set[lbf.x, lbf.y, z, lbf.i, lbf.j, k];
ltn ฌ Set[ltf.x, ltf.y, z, ltf.i, ltf.j, k];
rbn ฌ Set[rbf.x, rbf.y, z, rbf.i, rbf.j, k];
rtn ฌ Set[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 ฌ Set[lbn.x, lbn.y, z, lbn.i, lbn.j, k];
ltf ฌ Set[ltn.x, ltn.y, z, ltn.i, ltn.j, k];
rbf ฌ Set[rbn.x, rbn.y, z, rbn.i, rbn.j, k];
rtf ฌ Set[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];
};
vertices: Vertices ฌ NIL;
xEdges: REF Table ฌ NEW[Table[hashSize]];
yEdges: REF Table ฌ NEW[Table[hashSize]];
zEdges: REF Table ฌ NEW[Table[hashSize]];
centers: REF Table ฌ NEW[Table[hashSize]];
lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf: Corner;
cubes: LIST OF Cube;
step: REAL ~ 0.1*size;
cen: Triple ฌ FindStart[start, function, clientData, size, level
! Error => {error ฌ reason; GOTO Failed}];
lbn ฌ Set[cen.x-0.5*size, cen.y-0.5*size, cen.z-0.5*size, 0, 0, 0];
lbf ฌ Set[cen.x-0.5*size, cen.y-0.5*size, cen.z+0.5*size, 0, 0, 1];
ltn ฌ Set[cen.x-0.5*size, cen.y+0.5*size, cen.z-0.5*size, 0, 1, 0];
ltf ฌ Set[cen.x-0.5*size, cen.y+0.5*size, cen.z+0.5*size, 0, 1, 1];
rbn ฌ Set[cen.x+0.5*size, cen.y-0.5*size, cen.z-0.5*size, 1, 0, 0];
rbf ฌ Set[cen.x+0.5*size, cen.y-0.5*size, cen.z+0.5*size, 1, 0, 1];
rtn ฌ Set[cen.x+0.5*size, cen.y+0.5*size, cen.z-0.5*size, 1, 1, 0];
rtf ฌ Set[cen.x+0.5*size, cen.y+0.5*size, cen.z+0.5*size, 1, 1, 1];
cubes ฌ LIST[[0, 0, 0, lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf]];
WHILE cubes # NIL DO
c: Cube ฌ cubes.first;
debug ฌ c.i=2 AND c.j=2 AND c.k=4;
Poly[c];
cubes ฌ cubes.rest;
Propagate[c];
ENDLOOP;
EXITS Failed => RETURN[error];
};

debug: BOOL ฌ FALSE;

hashNBits: INTEGER ~ 5;
hashSize:	INTEGER ~ Basics.BITSHIFT[1, 3*hashNBits];
Hash: PROC [i, j, k: INTEGER] RETURNS [id: INTEGER] ~ INLINE {
mask: INTEGER ~ Basics.BITSHIFT[1, hashNBits]-1;
id ฌ Basics.BITOR[
Basics.BITSHIFT[
Basics.BITOR[
Basics.BITSHIFT[Basics.BITAND[ABS[i], mask], hashNBits],
Basics.BITAND[ABS[j], mask]],
hashNBits],
Basics.BITAND[ABS[k], mask]];
};

SetIJK: PROC [table: REF Table, i, j, k, id: INTEGER] RETURNS [b: BOOL ฌ FALSE] ~ {
index: INTEGER ฌ Hash[i, j, k];
list: LIST OF Entry ฌ table[index];
IF list = NIL THEN {table[index] ฌ LIST[[i, j, k, id]]; RETURN};
FOR l: LIST OF Entry ฌ list, l.rest WHILE l # NIL DO
IF l.first.i = i AND l.first.j = j AND l.first.k = k THEN {l.first.id ฌ id; RETURN[TRUE]};
ENDLOOP;
table[index] ฌ CONS[[i, j, k, id], list];
};

GetIJK: PROC [table: REF Table, i, j, k: INTEGER] RETURNS [id: INTEGER ฌ -1] ~ {
FOR l: LIST OF Entry ฌ table[Hash[i, j, k]], l.rest WHILE l # NIL DO
IF l.first.i = i AND l.first.j = j AND l.first.k = k THEN RETURN[l.first.id];
ENDLOOP;
};
AddToVertices: PROC [vertices: Vertices, v: Vertex] RETURNS [new: Vertices] ~ {
new ฌ IF vertices = NIL THEN NEW[ImplicitMinimal.VerticesRep[1]] ELSE vertices;
IF new.length = new.maxLength THEN {
old: Vertices ฌ new;
new ฌ NEW[ImplicitMinimal.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,
start: Triple ฌ [],
level, size, delta: REAL,
clientData: REF ANY ฌ NIL,
status: ImplicitMinimal.StatusProc ฌ NIL,
triangulate: BOOL ฌ FALSE,
adaptive: REAL ฌ 0.0]
RETURNS [s: G3dShape.Shape ฌ NIL]
~ {
Tri: ImplicitMinimal.TriangleProc ~ {
triangle: NatSequence ฌ NEW[G3dBasic.NatSequenceRep[3]];
triangle[0] ฌ i2; triangle[1] ฌ i1; triangle[2] ฌ i0; -- reverse triangle direction
triangle.length ฌ 3;
RETURN[SavePoly[triangle, vertices]];
};
Poly: ImplicitMinimal.PolygonProc ~ {
poly: NatSequence ฌ NEW[G3dBasic.NatSequenceRep[polygon.nVertices]];
FOR i: NAT IN [0..poly.length ฌ polygon.nVertices) DO
poly[i] ฌ polygon.ids[poly.length-1-i]; -- reverse polygon direction
ENDLOOP;
RETURN[SavePoly[poly, vertices]];
};
SavePoly: PROC [polygon: NatSequence, vertices: Vertices] RETURNS [cont: BOOL ฌ TRUE] ~ {
points ฌ vertices;
polys ฌ G3dBasic.AddToSurfaceSequence[polys, [NIL, polygon]];
IF status # NIL AND NOT status[polygon, vertices] THEN RETURN[FALSE];
};
points: ImplicitMinimal.Vertices ฌ NIL;
polys: G3dBasic.SurfaceSequence ฌ NIL;
tri: ImplicitMinimal.TriangleProc ฌ IF triangulate THEN Tri ELSE NIL;
poly: ImplicitMinimal.PolygonProc ฌ IF triangulate THEN NIL ELSE Poly;
error: ROPE ฌ Polygonize[function, start, level, size, delta, clientData, poly, tri, adaptive];
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 {
type: ATOM ฌ IF triangulate THEN $Triangulated ELSE $ConvexPolygon;
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,, type];
s.selected ฌ TRUE;
};
};

Triangulate: PUBLIC PROC [
shape: G3dShape.Shape,		-- polygonized shape to be triangulated
function: ValueProc,
level, size, delta: REAL,
clientData: REF ANY ฌ NIL]
RETURNS [s: G3dShape.Shape]
~ {
Opp: PROC [a, b: REAL] RETURNS [c: BOOL] ~ {c ฌ (a<0 AND b>0) OR (a>0 AND b<0)};
nTriangles: NAT ฌ 0;
val, step: REAL ฌ 0.1*size;
IF shape = NIL THEN RETURN[NIL];
FOR i: NAT IN [0..shape.surfaces.length) DO
nTriangles ฌ nTriangles+shape.surfaces[i].vertices.length;
ENDLOOP;
s ฌ NEW[G3dShape.ShapeRep ฌ [vertices: G3dShape.CopyVertexSequence[shape.vertices]]];
s.surfaces ฌ NEW[G3dBasic.SurfaceSequenceRep[nTriangles]];
FOR i: NAT IN [0..shape.surfaces.length) DO
p, q, n: Triple ฌ [0, 0, 0];
poly: NatSequence ฌ shape.surfaces[i].vertices;
FOR i: NAT IN [0..poly.length) DO
p ฌ G3dVector.Add[p, shape.vertices[poly[i]].point];
ENDLOOP;
val ฌ function[p ฌ G3dVector.Mul[p, 1.0/REAL[poly.length]], clientData];
n ฌ Normal[p, function, clientData, delta];
FOR j: NAT IN [0..10) DO
IF Opp[val, function[q ฌ G3dVector.ScaleRay[[p, n], -step*i], clientData]-level] OR
Opp[val, function[q ฌ G3dVector.ScaleRay[[p, n], step*i], clientData]-level] THEN {
position: Triple ฌ Converge[p, q, val, level, function, clientData, unknown];
normal: Triple ฌ Normal[position, function, clientData, delta];
v: G3dShape.Vertex ฌ NEW[G3dShape.VertexRep ฌ [position, normal]];
p0: NAT ฌ shape.vertices.length;
shape.vertices ฌ G3dShape.AddToVertexSequence[shape.vertices, v];
FOR p: NAT IN [0..poly.length) DO
nats: NatSequence ฌ NEW[G3dBasic.NatSequenceRep[3]];
nats.length ฌ 3;
nats[0] ฌ p0; nats[1] ฌ poly[p]; nats[2] ฌ poly[(p+1) MOD poly.length];
s.surfaces[s.surfaces.length].vertices ฌ nats;
s.surfaces.length ฌ s.surfaces.length+1;
ENDLOOP;
EXIT;
};
ENDLOOP;
ENDLOOP;
};
NextCWEdge: PROC [edge: Edge, face: Face] RETURNS [Edge] ~ {
RETURN[SELECT edge FROM
lb			=> IF face = l	THEN lf	ELSE bn,
lt			=> IF face = l	THEN ln	ELSE tf,
ln			=> IF face = l	THEN lb	ELSE tn,
lf			=> IF face = l	THEN lt	ELSE bf,
rb			=> IF face = r	THEN rn	ELSE bf,
rt			=> IF face = r	THEN rf	ELSE tn,
rn			=> IF face = r	THEN rt	ELSE bn,
rf			=> IF face = r	THEN rb	ELSE tf,
bn			=> IF face = b	THEN rb	ELSE ln,
bf			=> IF face = b	THEN lb	ELSE rf,
tn			=> IF face = t	THEN lt	ELSE rn,
ENDCASE	=> IF face = t	THEN rt	ELSE lf];
};

OtherFace: PROC [edge: Edge, face: Face] RETURNS [Face] ~ {
RETURN[SELECT edge FROM
lb			=> IF face = b	THEN l	ELSE b,
lt			=> IF face = l	THEN t	ELSE l,
ln			=> IF face = l	THEN n	ELSE l,
lf			=> IF face = f	THEN l	ELSE f,
rb			=> IF face = r	THEN b	ELSE r,
rt			=> IF face = t	THEN r	ELSE t,
rn			=> IF face = n	THEN r	ELSE n,
rf			=> IF face = r	THEN f	ELSE r,
bn			=> IF face = n	THEN b	ELSE n,
bf			=> IF face = b	THEN f	ELSE b,
tn			=> IF face = t	THEN n	ELSE t,
ENDCASE	=> IF face = f	THEN t	ELSE f];
};
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, unknown];
};

Converge: PROC [
pos, neg: Triple,
vpos, level: REAL,
function: ValueProc,
clientData: REF ANY,
direction: Face ฌ unknown]
RETURNS [p: Triple]
~ {
positive: BOOL;
nTries: NAT ฌ 0;
p ฌ pos;
IF vpos < 0 THEN {c: Triple ฌ neg; neg ฌ pos; pos ฌ c};
DO
SELECT direction FROM
n, f =>  p.z ฌ 0.5*(pos.z+neg.z);
b, t =>  p.y ฌ 0.5*(pos.y+neg.y);
l, r =>  p.x ฌ 0.5*(pos.x+neg.x);
ENDCASE => p ฌ G3dVector.Midpoint[pos, neg];
IF (nTries ฌ nTries+1) = 10 THEN EXIT;
positive ฌ (function[p, clientData]-level) > 0.0;
SELECT direction FROM
n, f => IF positive THEN pos.z ฌ p.z ELSE neg.z ฌ p.z;
b, t => IF positive THEN pos.y ฌ p.y ELSE neg.y ฌ p.y;
l, r => IF positive THEN pos.x ฌ p.x ELSE neg.x ฌ p.x;
ENDCASE => IF positive 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.
PrintCorner: PROC [name: ROPE, c: Corner] ~ {
TerminalIO.PutF["%g=[%g,%g,%g],value=%g\n",IO.rope[name],IO.real[c.x],IO.real[c.y],IO.real[c.z],IO.real[c.value]];
};

PrintCube: PROC [lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf: Corner] ~ {
FOR l: LIST OF Corner ฌ LIST[lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf], l.rest WHILE l # NIL DO
TerminalIO.PutF["[%g,%g,%g] ", IO.real[l.first.x], IO.real[l.first.y], IO.real[l.first.z]];
ENDLOOP;
TerminalIO.PutF["\n"];
};

Find: PROC [sign: {negative, positive}] RETURNS [c: Triple] ~ {
value, stepSize: REAL ฌ 0.1;
FOR n: NAT IN [0..5) DO
FOR longitude: REAL ฌ 0.0, longitude+30.0 WHILE longitude <= 360.0 DO
FOR lattitude: REAL ฌ 0.0, lattitude+30.0 WHILE lattitude <= 360.0 DO
c ฌ G3dVector.CartesianFromPolar[longitude, lattitude, n*stepSize];
value ฌ function[c, clientData]-level;
IF sign = negative AND value <= 0.0 THEN RETURN;
IF sign = positive AND value >= 0.0 THEN RETURN;
ENDLOOP;
ENDLOOP;
ENDLOOP;
ERROR startError;
};

๔
ImplicitMinimalImpl.mesa
Copyright c 1991 by Xerox Corporation.  All rights reserved.
Bloomenthal, August 11, 1992 4:02 pm PDT
Types

l(left): -x,-i; r(right): +x,+i; b(bottom): -y,-j; t(top): +y,+j; n(near): -z,-k; f(far): +z,+k
the lbn corner of cube (i, j, k) corresponds with location (i, j, k).
A `Minimal' Polygonizer
negative to positive along edge is cw about face
Return -1 if no vertex on edge; otherwise, return vertex (compute & save if unknown)
IF ABS[p1.value] < tiny AND ABS[p2.value] < tiny THEN RETURN[-1]; -- no vertex
tiny: REAL ฌ 1.0e-3;
simple, due to cube symmetry
PrintCube: PROC [c: Cube, title: ROPE] ~ {IF debug THEN TerminalIO.PutF["%g i, j, k= %g, %g, %g\n", IO.rope[title], IO.int[c.i], IO.int[c.j], IO.int[c.k]]};
Sparse 3d Storage

Courtesy Brian Wyvill
Return true if already set; otherwise, set and return false.
Get integer at i, j, k; id = -1 if [i, j, k] not set.
Miscellany
Topology
Geometry
opposite usual gradient so surface normals point outwards
Debug and/or Old
ส‚•NewlineDelimiter

™š
œ™Jšœ
ฯm
œ1™<J™(—J˜Jšฯk	œZ˜cJ˜šัblnœžœ
ž˜"Jšžœ:˜AJšžœ˜—J˜Jšœž˜
headš
ฯl™Jšœ	žœ˜ Jšœ
žœ˜)Jšœžœ˜)Jšœ	žœ˜'Jšœžœ˜+šžœžœžœ
˜J˜—Jšœžœžœžœฯs
œก
œก
œก
œก
œก
œก
œก
œ˜YJšœ	žœžœก
œก
œก
žก
œ
ก
œก
œก
œก
œก
œก
žก
œ
ก
œ˜HJšœžœ˜+Jšœžœ4˜@Jšœžœžœžœ˜-Jšœ	žœžœžœ	žœ
žœ
žœ
žœ˜GJ™J™_Jš
œE™E—š
 ™šฯnœžœžœ
˜J˜—šข
œžœ
žœ˜Jš
œ˜J˜Jšœžœ
˜Jšœžœ
˜Jšœžœ
˜Jšœžœ
žœžœ
˜Jšœ-žœ
˜1Jšœ/žœ
˜3Jšœ
žœ˜Jšžœ	žœžœ
˜Jš
œ˜Jšžœžœ	˜1š
ขœ
ก
žก
œ	ก
œ
ก
œ
˜šขœžœžœ˜2J™0šžœ
žœž˜Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšœžœžœžœ˜+Jšžœžœžœžœ˜/—J˜—šขœžœžœžœ˜2J™TJš
œ
˜
Jšœžœ˜Jšœฯc?˜JJ˜šžœž˜
J˜8J˜8J˜8J˜8J˜8J˜8J˜8J˜8J˜8J˜8J˜8J˜8Jšžœ
˜—Jš
žœ
žœžœ
žœžœ
žœฃ™NJšžœ%žœ
žœฃ˜DJšžœ-žœ
žœ˜?Jšœ
ก
œก
œก
œก
œก
œก
œก
œ	ก
œก
œ	ก
œก
œ˜jJ˜;J˜&J˜=J˜—šข
œžœžœžœžœ
žœ˜Pšข	œžœ˜J˜PJ˜—Jšขœ
ก
žก
œก
žœ
ก
žก
œก
žœ	ก
žก
œก
žก
œก
žก
œ˜KJšขœžœžœ/˜NJšขœžœžœžœ)˜OJšœžœ
˜
J˜#Jšœ
žœ˜$Jšœžœžœ
žœ
˜Jšžก
œก
žœ
žœžœ*žœ
˜QJšœ žœ˜>J˜1šžœžœ
žœ	ž˜š	žก
œก
œก
œ"ก
ž˜Gš	œก
œก
œก
œžœ˜IJš	œก
œก
œก
œก
œ˜>Jšžœ
˜J˜——Jšžœ
žœ˜!Jšžœ
˜—J˜šžœžœ
ž˜šžœžœ
žœž˜šžก
žก
œก
œก
œ	ก
žก
œก
œ	˜KJšžก
žก
œ˜—šžœžœ˜J˜Jšœžœ˜#Jšžœžœ
ฯb
œ
˜5J˜—Jšžœ
˜——šžœžœ
žœ˜šžœžœ
žœž˜J˜Jšœžœ˜#Jšžก
žก
œก
œ
žœ
žœ˜;Jšžœžœžœ
ค
œ
˜CJšžœ
˜—J˜—J˜—šข	œ
ก
žก
œก
žœ
ก
œก
žœก
žก
žœ
ก
žก
œ
žœ
ก
œ
ก
œ
˜UJ˜!šžœžœ
žœž˜J˜>Jšžœ
˜—J˜&šžœžœ
žœž˜Jšžœ>žœ
žœ
žœ˜RJšžœ
˜—Jšžœ
žœ˜J˜—Jšœžœ
™Jšœžœ
˜Jšœžœžœ
žœ
˜Jšœžœžœ
žœžœ
žœ˜&šžœžœž˜šžœ
žœžœ žœ˜;J˜J˜J˜š
ž˜Jšœ žœ
˜%Jšžœžœ
žœ
˜šžœžœ˜!Jšžœžœ
žœฃะcs˜&J˜J˜J˜J˜—Jšžœ
˜—Jšžœžœ
žœฃฅ˜&Jšžœžœ
žœ"˜>šžœžœ
žœ˜šžœžœžœ
žœ˜HJšžœ!˜%Jš
žœ
žœ
žœ+žœ
žœ˜G—J˜—J˜—Jšžœ
˜—J˜—š
ขœžœžœžœžœ
˜@Jšžœ<˜BJ˜—šขœ
ก
žก
œก
œก
œก
žœ
ก
œก
œก
œก
œก
œก
œก
œก
œ
˜IJ˜Jšœžœ˜Jš
œ/˜/Jš(žก
œ	ก
œ
ก
œก
œ
ก
œก
žก
œ	ก
œ
ก
œก
œ
ก
œก
žก
œ	ก
œ
ก
œก
œ
ก
œก
žก
žœ˜ZJšžœžœ
žœ
˜+šžœž˜šœฃK˜RJšœžœ
˜Jšœžœ
˜J˜'J˜,J˜,J˜,J˜,J˜—šœฃK˜RJšœžœ
˜Jšœžœ
˜J˜'J˜,J˜,J˜,J˜,Jš
œ˜—šœฃP˜WJšœžœ
˜Jšœžœ
˜J˜'J˜,J˜,J˜,J˜,J˜—šœฃK˜RJšœžœ
˜Jšœžœ
˜J˜'J˜,J˜,J˜,J˜,J˜—šœฃI˜PJšœžœ
˜Jšœžœ
˜J˜'J˜,J˜,J˜,J˜,J˜—šœฃI˜PJšœžœ
˜Jšœžœ
˜J˜'J˜,J˜,J˜,J˜,J˜—Jšžœ
˜—J˜9Jš
œžœ	žœ
žœ
žœžœ
žœ˜>J˜—šข	œžœ˜J™Jšคœ1˜9Jšคœ1˜9Jšคœ1˜9Jšคœ1˜9Jšคœ1˜9Jšคœ1˜9J˜—Jšœžœ
˜Jšœžœ	žœ˜)Jšœžœ	žœ˜)Jšœžœ	žœ˜)Jšœ	žœ	žœ˜*Jš
œ/˜/Jšœžœ
žœ˜Jšœžœ˜˜@Jšœžœ
˜*—Jšœ)ก
œ˜CJšœ)ก
œ˜CJ˜CJ˜CJšœ)ก
œ˜CJšœ)ก
œ˜CJ˜CJ˜CJšœžœ4˜@šžœ	žœ
ž˜J˜Jšœžœžœ˜"Jš
œ˜J˜Jš
œ
˜
Jšžœ
˜—Jšžœžœ˜J˜J˜—šœžœžœ
˜J˜—Jšข	œžœžœžœ-žœžœžœžœ™œ
—š
 ™Jšœ
ก
žœ˜Jšœ
žœ
žœ˜4J™šขœžœžœžœžœžœ˜>J™Jšœžœ
žœ˜0šœžœ
˜šœžœ
˜šœžœ
˜
Jšœžœžœ
žœ˜8Jšœžœ
žœ˜—Jš
œ˜—Jšœžœ
žœ˜—J˜J˜—šขœžœ	žœžœžœžœžœ˜SJ™<Jšœžœ˜Jšœžœ
žœ˜#Jš
žœžœ
žœžœžœ˜@šžœžœ
žœžœžœ
ž˜4Jšžœžœžœžœžœ
žœ˜ZJšžœ
˜—Jšœžœ˜)Jš
œ˜J˜—šขœžœ	žœžœžœžœ
˜PJ™5šžœžœ
žœ&žœžœ
ž˜DJš
žœžœžœžœ
žœ
˜MJšžœ
˜—J˜——š
 
™

codešข
œžœ!žœ˜OLšœžœžœ
žœ
žœ!žœ
˜Ošžœžœ˜$L˜Jšœžœžœ'˜PJš
žœžœ
žœžœžœ
˜GL˜—L˜L˜L˜L˜—š
ขœžœ
žœ
žœžœ
˜*J˜—šขœžœ
žœ˜ Lš
œ˜L˜Lšœžœ
˜Lšœžœ
žœžœ
˜Lšœ%žœ
˜)Lšœ
žœžœ
˜Lšœ
žœ˜Lšžœžœ
˜!Lš
œ˜šขœ"˜%Jšœžœ˜8Jšœ6ฃ˜SJ˜Jšžœ˜%Jš
œ˜—šขœ!˜%Jšœžœ-˜Dšžก
œก
žก
žก
œ$ก
ž˜5Jšœก
œ
ก
œฃ˜DJšžœ
˜—Jšžœ˜!Jš
œ˜—šขœžœ+ก
žก
œก
žก
œ
ก
žœ
ก
œ
ก
œ
˜YJ˜Jšœ.žœ˜=Jšžœ
žœ
žœ
žœžœ
žœ
žœ˜EJ˜—Jšœ#žœ
˜'Jšœ"žœ
˜&Jš	œ$žœ
žœžœ
žœ
˜EJš	œ$žœ
žœ
žœ
žœ˜FJšœžœT˜_šžœ
žœ
žœž˜)Jšœžœ	žœ
žœžœ
žœ˜K—šžœ	žœ
žœž˜'Jšœžœ	žœ
žœžœ
žœ˜K—šžœ	ž˜Jšžœ
˜šžœ˜Jš	œžœžœ
žœžœ˜CJšœ	ก
œก
œ
ก
žœ,˜SJšœก
œก
œ
ก
žœ,˜Ršžœžœ
žœ7ž˜FJ˜!J˜Jšžœ
˜—Jšœžœ
ก
œก
œ	ก
œžœ	˜KJšœ
žœ
˜J˜——J˜J˜—šขœžœ
žœ˜Lšœฃ'˜?Lš
œ˜Lšœžœ
˜Lšœžœ
žœžœ
˜Lšžœ˜Lš
œ˜Jšขœ
ก
žก
œก
žœ
ก
žก
œก
žœ
ก
žก
œก
žก
œก
žก
œ˜PJšœžœ˜Jšœžœ˜Jš
žœ	žœ
žœ
žœ
žœ˜ šžœžœ
žœž˜+L˜:Lšžœ
˜—JšœžœN˜ULšœ
žœ*˜:šžœžœ
žœž˜+J˜L˜/š	žก
œก
žœ
žœž˜!J˜4Jšžœ
˜—Jšœ(žœ˜HJ˜+šžœžœ
žœ	ž˜š	žก
œก
œ#ก
œ ก
ž˜Sšœก
œ
ก
œ
ก
œก
œก
œก
œก
žก
œ
˜SJš
œก
œ
ก
œก
œก
œก
œก
œ˜MJ˜?Jšœžœ*˜BJšœžœ˜ J˜Ašžœžœ
žœž˜!Lšœžœ˜4L˜Lšœ6žœ˜GL˜.L˜(Jšžœ
˜—Jšžœ
˜J˜——Lšžœ
˜—Lšžœ
˜—J˜——š
 ™šข
œžœžœ˜<šžœ
žœž˜Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšœžœ
žœžœ˜$Jšžœžœ
žœžœ˜(—J˜J˜—šข	œžœžœ˜;šžœ
žœž˜Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšœžœ
žœžœ˜"Jšžœžœ
žœžœ˜&—J˜——š
 ™šข	œžœ˜Jš
œ˜Jš
œ˜Jšœžœ
žœ
˜Jšœ
žœ
˜Jšžœ˜Jš
œ˜Jšœžœžœžœ
˜,š
ขœžก
œก
œ
ก
œžœžœ˜VJšœžœ
˜šžœžœ
žœž˜Jšœ/žœ
žœ
žœ˜DJšœ/žœ
žœ
žœ˜DJšœ/žœ
žœ
žœ˜DJ˜,Jšžœžœžœžœžœ
žœ
˜RJ˜Jšžœ
˜—Jšœžœ
˜J˜—J˜*J˜J˜šžœ	žœ˜Jšžœ žœ	žœžœ˜VJšžœG˜K—J˜J˜—šขœ
ก
žก
œ
˜Jšœก
œก
œ˜Jšœก
žœ
˜Jšœ	ก
œ
˜Jšœžœ
žœ
˜Jšœ
ก
œ˜Jšžœ˜Jš
œ˜Jšœ
žœ
˜Jšœžœ˜J˜Jšžก
œก
œ
ก
œ
ก
žก
œก
œก
œ
ก
œก
œ
ก
œ	˜7š
ž˜šžœž˜J˜!J˜!J˜!Jšžœ%˜,—Jšžœžœ
žœ
˜&J˜1šžœž˜Jšœžœ
žœ
žœ
˜6Jšœžœ
žœ
žœ
˜6Jšœžœ
žœ
žœ
˜6Jšžœžœ
žœ	žœ	˜1—Jšžœ
˜—J˜J˜—š
ขœžœ2žœ
žœ	žœ
˜SJšžœ˜Jš
œ˜Jš
œ9™9Jšœžœ˜&Jšœžœ=˜DJšœžœ=˜DJšœžœ=˜DJ˜J˜——J˜Jšžœ
˜š
 ™šขœžœžœ˜-Jšœ+žœžœžœžœžœ˜rJ˜J˜—šข	œžœ5˜Dš
žœžœ
žœ
žœ1žœžœ
ž˜]Jšœžœžœžœ˜[Jšžœ
˜—J˜J˜J˜—šขœžก
œก
œ
ก
œžœ˜?Jšœžœ˜šžœžœ
žœž˜šžœžœžœž˜Ešžœžœžœž˜EJ˜CJ˜&Jšžœžœžœ
žœ
˜0Jšžœžœžœ
žœ
˜0Jšžœ
˜—Jšžœ
˜—Jšžœ
˜—Jšžœ˜J˜——J˜—…—LาkH