[_CDCSL_93-16_]<1>Cedar>release>ImplicitSurfaces>ImplicitMinimalOld2Impl.mesa

ImplicitMinimalImpl.mesa

Bloomenthal, August 11, 1992 4:02 pm PDT

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

ImplicitMinimalImpl: CEDAR PROGRAM

IMPORTS Basics, G3dBasic, G3dShape, G3dVector, Random, Real, Rope

EXPORTS ImplicitMinimal

~ BEGIN

Types

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];

IntegerArray: TYPE ~ ImplicitMinimal.IntegerArray;

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

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, dir: BOOL] RETURNS [Face] ~ {

return face about which negative to positive direction

along edge is clockwise if dir or counter-clockwise if not dir.

Sign: PROC [value: REAL] RETURNS [b: BOOL] ~ {

b Ź (dir AND value > 0.0) OR (NOT dir AND value < 0.0);

};

RETURN[SELECT edge FROM

lb => IF Sign[c.lbn.value] THEN l ELSE b,

lt => IF Sign[c.ltn.value] THEN t ELSE l,

ln => IF Sign[c.lbn.value] THEN n ELSE l,

lf => IF Sign[c.lbf.value] THEN l ELSE f,

rb => IF Sign[c.rbn.value] THEN b ELSE r,

rt => IF Sign[c.rtn.value] THEN r ELSE t,

rn => IF Sign[c.rbn.value] THEN r ELSE n,

rf => IF Sign[c.rbf.value] THEN f ELSE r,

bn => IF Sign[c.lbn.value] THEN b ELSE n,

bf => IF Sign[c.lbf.value] THEN f ELSE b,

tn => IF Sign[c.ltn.value] THEN n ELSE t,

ENDCASE => IF Sign[c.ltf.value] THEN t ELSE f];

};

VertexID: PROC [e: Edge] RETURNS [id: INTEGER] ~ {

Return -1 if no vertex on edge; otherwise, return vertex (compute & save if unknown)

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 ABS[p1.value] < tiny AND ABS[p2.value] < tiny THEN RETURN[-1]; -- no vertex

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: IntegerArray] ~ {

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: IntegerArray;

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 [nVerts: INTEGER, ids: IntegerArray] 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..nVerts) DO

IF G3dVector.Dot[aveNormal, vertices[ids[i]].normal] < adaptive THEN RETURN[TRUE];

ENDLOOP;

RETURN[FALSE];

};

tiny: REAL Ź 1.0e-3;

id, nVertices, nPos: INTEGER Ź 0;

ids: IntegerArray;

done: ARRAY Edge OF BOOL Ź ALL[FALSE];

FOR l: LIST OF Corner Ź LIST[c.lbn, c.lbf, c.ltn, c.ltf, c.rbn, c.rbf, c.rtn, c.rtf],

l.rest WHILE l # NIL DO

IF l.first.value > 0.0 THEN nPos Ź nPos+1;

ENDLOOP;

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, nPos > 4];

nVertices Ź 1;

done[e Ź NextCWEdge[e, face]] Ź TRUE;

IF e = start THEN EXIT;

IF (id Ź VertexID[e]) # -1 THEN {

IF nVertices = 6 THEN ERROR; -- RETURN;

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] ~ {

simple, due to cube symmetry

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;

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

hashNBits: INTEGER ~ 5;

hashSize: INTEGER ~ Basics.BITSHIFT[1, 3*hashNBits];

Hash: PROC [i, j, k: INTEGER] RETURNS [id: INTEGER] ~ INLINE {

Courtesy Brian Wyvill

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] ~ {

Return true if already set; otherwise, set and return 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] ~ {

Get integer at i, j, k; id = -1 if [i, j, k] not set.

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;

};

Miscellany

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;

};

Topology

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];

};

Geometry

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};

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]

~ {

opposite usual gradient so surface normals point outwards

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.

Debug and/or Old

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;

ERROR startError;

};