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

ImplicitMinimalImpl.mesa

Bloomenthal, November 21, 1992 4:22 pm PST

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

ImplicitTetImpl: 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 ~ 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];

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,

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

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

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

corners presumed oriented such that b, c, d appear clockwise when viewed from a

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

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

};

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

};

Sparse 3d Storage

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

Courtesy Brian Wyvill

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

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

};

Miscellany

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;

};

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

};

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

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]

~ {

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

Hash6d: PROC [i1, j1, k1, i2, j2, k2: INTEGER] RETURNS [id: INT] ~ {

IF i1 > i2 OR (i1 = i2 AND (j1 > j2 OR (j1 = j2 AND k1 > k2))) THEN {

t: INTEGER Ź i1; i1 Ź i2; i2 Ź t; t Ź j1; j1 Ź j2; j2 Ź t; t Ź k1; k1 Ź k2; k2 Ź t;

};

id Ź Basics.Int32FromF[

[Basics.HFromInt16[Hash3d[i1, j1, k1]], Basics.HFromInt16[Hash3d[i2, j2, k2]]]];

};