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

ImplicitMinimalImpl.mesa

Bloomenthal, March 10, 1993 11:45 am PST

DIRECTORY G3dBasic, G3dPolygon, G3dShape, G3dTriangle, G3dVector, ImplicitCubeTet, ImplicitDebug, ImplicitMinimal, ImplicitStorage, IO, Random, Rope;

ImplicitMinimalImpl: CEDAR PROGRAM

IMPORTS G3dBasic, G3dPolygon, G3dShape, G3dVector, ImplicitCubeTet, ImplicitDebug, ImplicitStorage, IO, Random, Rope

EXPORTS ImplicitMinimal

~ BEGIN

Types and Errors

Triple: TYPE ~ G3dBasic.Triple;

Cube: TYPE ~ ImplicitCubeTet.Cube;

CubeProc: TYPE ~ ImplicitCubeTet.CubeProc;

Corner: TYPE ~ ImplicitCubeTet.Corner;

Tetrahedron: TYPE ~ ImplicitCubeTet.Tetrahedron;

TetProc: TYPE ~ ImplicitCubeTet.TetProc;

TwoCorners: TYPE ~ ImplicitCubeTet.TwoCorners;

Function: TYPE ~ ImplicitMinimal.Function;

TriangleProc: TYPE ~ ImplicitMinimal.TriangleProc;

Centers: TYPE ~ ImplicitStorage.Centers;

IDs: TYPE ~ ImplicitStorage.IDs;

Edges: TYPE ~ ImplicitStorage.Edges;

Faces: TYPE ~ ImplicitStorage.Faces;

Vertices: TYPE ~ ImplicitStorage.Vertices;

ROPE: TYPE ~ Rope.ROPE;

Error: PUBLIC ERROR [reason: ROPE] = CODE;

Cubic Partitioning of Implicit Surface

Propagate: PUBLIC PROC [

function: Function,

level, size: REAL,

clientData: REF ANY Ź NIL,

start: Triple Ź [],

action: CubeProc]

RETURNS [centers: REF Centers]

~ {

IF action # NIL THEN {

Inner: PROC [c: REF Cube] ~ {

TestFace: PROC [i, j, k: INT, c1, c2, c3, c4: REF Corner, face: {l, r, b, t, n, f}] ~ {

Set: PROC [x, y, z: REAL, i, j, k: INT] RETURNS [c: REF Corner] ~ {

c Ź NEW[Corner Ź [[i, j, k], [x, y, z], 0.0, unset, FALSE]];

SetCornerValue[c];

};

newCube: REF Cube;

state: ImplicitCubeTet.CornerState Ź c1.state;

lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf: REF Corner;

IF c2.state = state AND c3.state = state AND c4.state = state THEN RETURN; -- do this

IF ImplicitStorage.AddCenter[centers, [i, j, k]] THEN RETURN; -- before this

SELECT face FROM

l => { -- have right (rbn, rbf, rtn, rtf); make left (lbn, lbf, ltn, ltf)

x: REAL Ź c1.pos.x-size;

i: INT Ź c1.id.i-1;

rbn Ź c1; rbf Ź c2; rtn Ź c3; rtf Ź c4;

lbn Ź Set[x, rbn.pos.y, rbn.pos.z, i, rbn.id.j, rtn.id.k];

lbf Ź Set[x, rbf.pos.y, rbf.pos.z, i, rbf.id.j, rbf.id.k];

ltn Ź Set[x, rtn.pos.y, rtn.pos.z, i, rtn.id.j, rtn.id.k];

ltf Ź Set[x, rtf.pos.y, rtf.pos.z, i, rtf.id.j, rtf.id.k];

};

r => { -- have left (lbn, lbf, ltn, ltf); make right (rbn, rbf, rtn, rtf)

x: REAL Ź c1.pos.x+size;

i: INT Ź c1.id.i+1;

lbn Ź c1; lbf Ź c2; ltn Ź c3; ltf Ź c4;

rbn Ź Set[x, lbn.pos.y, lbn.pos.z, i, lbn.id.j, lbn.id.k];

rbf Ź Set[x, lbf.pos.y, lbf.pos.z, i, lbf.id.j, lbf.id.k];

rtn Ź Set[x, ltn.pos.y, ltn.pos.z, i, ltn.id.j, ltn.id.k];

rtf Ź Set[x, ltf.pos.y, ltf.pos.z, i, ltf.id.j, ltf.id.k];

};

b => { -- have top (ltn, ltf, rtn, rtf); make bottom (lbn, lbf, rbn, rbf)

y: REAL Ź c1.pos.y-size;

j: INT Ź c1.id.j-1;

ltn Ź c1; ltf Ź c2; rtn Ź c3; rtf Ź c4;

lbn Ź Set[ltn.pos.x, y, ltn.pos.z, ltn.id.i, j, ltn.id.k];

lbf Ź Set[ltf.pos.x, y, ltf.pos.z, ltf.id.i, j, ltf.id.k];

rbn Ź Set[rtn.pos.x, y, rtn.pos.z, rtn.id.i, j, rtn.id.k];

rbf Ź Set[rtf.pos.x, y, rtf.pos.z, rtf.id.i, j, rtf.id.k];

};

t => { -- have bottom (lbn, lbf, rbn, rbf); make top (ltn, ltf, rtn, rtf)

y: REAL Ź c1.pos.y+size;

j: INT Ź c1.id.j+1;

lbn Ź c1; lbf Ź c2; rbn Ź c3; rbf Ź c4;

ltn Ź Set[lbn.pos.x, y, lbn.pos.z, lbn.id.i, j, lbn.id.k];

ltf Ź Set[lbf.pos.x, y, lbf.pos.z, lbf.id.i, j, lbf.id.k];

rtn Ź Set[rbn.pos.x, y, rbn.pos.z, rbn.id.i, j, rbn.id.k];

rtf Ź Set[rbf.pos.x, y, rbf.pos.z, rbf.id.i, j, rbf.id.k];

};

n => { -- have far (lbf, ltf, rbf, rtf); make near (lbn, ltn, rbn, rtn)

z: REAL Ź c1.pos.z-size;

k: INT Ź c1.id.k-1;

lbf Ź c1; ltf Ź c2; rbf Ź c3; rtf Ź c4;

lbn Ź Set[lbf.pos.x, lbf.pos.y, z, lbf.id.i, lbf.id.j, k];

ltn Ź Set[ltf.pos.x, ltf.pos.y, z, ltf.id.i, ltf.id.j, k];

rbn Ź Set[rbf.pos.x, rbf.pos.y, z, rbf.id.i, rbf.id.j, k];

rtn Ź Set[rtf.pos.x, rtf.pos.y, z, rtf.id.i, rtf.id.j, k];

};

f => { -- have near (lbn, ltn, rbn, rtn); make far (lbf, ltf, rbf, rtf)

z: REAL Ź c1.pos.z+size;

k: INT Ź c1.id.k+1;

lbn Ź c1; ltn Ź c2; rbn Ź c3; rtn Ź c4;

lbf Ź Set[lbn.pos.x, lbn.pos.y, z, lbn.id.i, lbn.id.j, k];

ltf Ź Set[ltn.pos.x, ltn.pos.y, z, ltn.id.i, ltn.id.j, k];

rbf Ź Set[rbn.pos.x, rbn.pos.y, z, rbn.id.i, rbn.id.j, k];

rtf Ź Set[rtn.pos.x, rtn.pos.y, z, rtn.id.i, rtn.id.j, k];

};

ENDCASE;

newCube Ź ImplicitCubeTet.SetCube[lbn.id, lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf, start];

IF NOT action[newCube, NIL] THEN RETURN;

cubes Ź CONS[newCube, cubes];

};

TestFace[c.id.i-1, c.id.j, c.id.k, c.lbn, c.lbf, c.ltn, c.ltf, l];

TestFace[c.id.i+1, c.id.j, c.id.k, c.rbn, c.rbf, c.rtn, c.rtf, r];

TestFace[c.id.i, c.id.j-1, c.id.k, c.lbn, c.lbf, c.rbn, c.rbf, b];

TestFace[c.id.i, c.id.j+1, c.id.k, c.ltn, c.ltf, c.rtn, c.rtf, t];

TestFace[c.id.i, c.id.j, c.id.k-1, c.lbn, c.ltn, c.rbn, c.rtn, n];

TestFace[c.id.i, c.id.j, c.id.k+1, c.lbf, c.ltf, c.rbf, c.rtf, f];

};

SetCornerValue: PROC [c: REF Corner] ~ {

c.state Ź IF (c.value Ź (function[c.pos, clientData]-level)) < 0.0

THEN negative ELSE positive;

c.on Ź ABS[c.value] < 0.000001; -- constant unrelated to tet size

};

pt: Triple Ź IF start = [] THEN FindStart[start, function, size, level, clientData] ELSE start;

id: G3dTriangle.IntTriple Ź ImplicitCubeTet.IDFromPoint[pt, size, pt];

c: REF Cube Ź ImplicitCubeTet.MakeCube[id, size, pt];

cubes: LIST OF REF Cube Ź LIST[c];

SetCornerValue[c.lbn]; SetCornerValue[c.lbf]; SetCornerValue[c.ltn]; SetCornerValue[c.ltf];

SetCornerValue[c.rbn]; SetCornerValue[c.rbf]; SetCornerValue[c.rtn]; SetCornerValue[c.rtf];

centers Ź ImplicitStorage.NewCenters[];

WHILE cubes # NIL DO -- need list (can't put all cubes on process stack)

c: REF Cube Ź cubes.first;

cubes Ź cubes.rest;

Inner[c];

ENDLOOP;

};

A Minimal Polygonizer

Abort: ERROR = CODE;

nOn: INT Ź 0;

Triangulate: PUBLIC PROC [

function: Function,

level, size: REAL,

clientData: REF ANY Ź NIL,

start: Triple Ź [],

triangleAction: TriangleProc]

RETURNS [error: ROPE Ź NIL]

~ {

ENABLE Abort => {error Ź "aborted"; GOTO Bad};

CubeAction: CubeProc ~ {

TetAction: TetProc ~ {

tetID Ź tetID+1;

vertices Ź ProcessTet[

function, clientData, level, tet, tetID, ids, edges, vertices, epsilon, triangleAction];

};

ImplicitCubeTet.Decompose[cube, TetAction, scratchTet];

};

ids: REF IDs Ź ImplicitStorage.NewIDs[];

edges: REF Edges Ź ImplicitStorage.NewEdges[];

tetID: INT Ź nOn Ź 0;

epsilon: REAL Ź 0.001*size;

vertices: Vertices Ź NEW[G3dShape.VertexSequenceRep[10]];

scratchTet: REF Tetrahedron Ź NEW[Tetrahedron];

ImplicitDebug.SetDebug[[edges: edges]];

vertices.valid[normal] Ź TRUE;

[] Ź Propagate[function, level, size, clientData, start, CubeAction];

EXITS Bad => NULL;

};

ProcessTet: PUBLIC PROC [

function: Function,

clientData: REF ANY,

level: REAL,

tet: REF Tetrahedron,

tetID: INT,

ids: REF IDs,

edges: REF Edges,

vertices: Vertices,

epsilon: REAL,

action: TriangleProc]

RETURNS [Vertices]

~ {

VIDfromEdge: PROC [e: ImplicitCubeTet.TetEdge] RETURNS [vid: INT Ź -1] ~ {

VIDfromCorners: PROC [c: TwoCorners] RETURNS [vid: INT Ź -1] ~ {

AddVertex: PROC [p: Triple] RETURNS [vid: INT] ~ {

v: G3dShape.Vertex Ź NEW[G3dShape.VertexRep Ź [point: p]];

v.normal Ź Normal[p, function, epsilon, clientData];

vertices Ź G3dShape.AddToVertexSequence[vertices, v];

ImplicitDebug.SetDebug[[vertices: vertices]];

ImplicitStorage.SetEdge[edges, c.a.id, c.b.id, vid Ź vertices.length-1];

};

VIDfromCorner: PROC [c: REF Corner] RETURNS [vid: INT] ~ {

IF (vid Ź ImplicitStorage.GetID[ids, c.id]) = -1

THEN ImplicitStorage.SetID[ids, c.id, vid Ź AddVertex[c.pos]];

};

SELECT TRUE FROM

corner.on is a mechanism whereby small triangles can be eliminated; the idea is to subordinate an edge id to a corner id if the corner if the distance between corner and surface is negligible. Consider, for example, the surface intersecting the tetrahedron near one of its corners. This may occur in one of three ways; the first is when a single surface/edge intersection is negligibly close to the tetrahedron corner; in effect, this simply shifts slightly the location of the intersection. The other two cases are when either two or three of the intersections are negligibly close to the corner, as shown below:

[Artwork node; type 'Artwork on' to command tool]

The transition on the left (indicated by the red arrow) shows three blue surface/edge intersections, each of which are considered negligibly close to the tetrahedron corner; this closeness is indicated by the red outline of the blue circles. After the transition, these three circles are subordinated to the corner itself, shown as a red circle. Triangle 1 has become degenerate, with all three of its vertex ids being equal to the corner vertex id. Triangle 2 is also degenerate, with two of its vertex ids being equal to the corner vertex id. Triangle 3 remains non-degenerate. The transition on the left differs in that the middle blue vertex is not considered negligibly close to the corner; thus, upon the transition, Triangle 1 becomes degenerate, with two of its vertex ids equal to the corner id. Triangles 2 and 3 remain non-degenerate.

The green circles represent surface/edge intersections from another tetrahedron; Triangle 2 on the left contains one such intersection; initially it is non-degenerate but after the transition becomes degenerate. Topological consistency is retained, however, from one tetrahedra to the next.

There are two important consequences of this scheme. The first is that whether a corner is considered `on' or not must be determined consistently; that is, if a particular corner is computed first for one tetrahedron and then later for a different tetrahedron, whether it is `on' or not must be determined equally both times (if this consistency is not maintained, then redundant vertices will be calculated and the Euler relation will not hold). The second consequence is that, if a corner is on, it must must yield an identical surface id for all those edges connecting to it.

The first consequence is simply accommodated for an implicit function: namely, the implicit function remains constant at the given corner and the test
on = function (corner) < epsilon
should yield consistent results for a given corner. Note that epsilon is not a unit of distance, and should bear no relation to the tetrahedron size; rather, it is purely a mathematical quantity.

The second consequence is somewhat subtler. It is tempting to associate the vertex id with the corner itself:
corner.vid = AddVertex[corner.position]
However, this presumes that each time the corner is visited, the same vid is available, which implies that the corner is retained in memory. This is undesirable, in that each corner includes the following information:
locationID: IntTriple,
position: Triple,
value: REAL,
state: {positive, negative},
on: BOOLEAN
Rather, we prefer to propagate the corners along with the cube propagation, computing four new corners for each new cube; the number of corners held in memory is thus a function of the number of cubes placed on the cube stack, but this is, in general, significantly less than the total number of corners created during the entire propagation. This does mean some corners are recomputed. For example, in the figure below, the initial cube, marked `1,' consists of the eight red circles. In propagating to the right and to above, the blue circles are computed. If the upper left cube propagates to the right, then the blue circles outlined in orange are re-computed; if the lower right cube propagates upward, then the blue circles outlined in yellow are re-computed.

[Artwork node; type 'Artwork on' to command tool]

The edges , however, cannot rely strictly on the six-dimensional hashing scheme to determine vertex ids that are associated with a corner, if that corner is not constantly resident in memory. Thus, an additional hash table for corner ids is required; only those corners negligibly close to the surface need store their associated vertex id in the hash table; when determining the vertex id for an intersected edge, if one of the edge corners is negligibly close, then the hash table for corner vertex ids is consulted.

c.a.on AND ids # NIL => {vid Ź VIDfromCorner[c.a]; nOn Ź nOn+1};

c.b.on AND ids # NIL => {vid Ź VIDfromCorner[c.b]; nOn Ź nOn+1};

c.a.state # c.b.state =>

IF (vid Ź ImplicitStorage.GetEdge[edges, c.a.id, c.b.id]) = -1

THEN vid Ź AddVertex[

Converge[c.a.pos, c.b.pos, c.a.value, level, function, clientData]];

ENDCASE;

};

vid Ź VIDfromCorners[ImplicitCubeTet.CornersFromEdge[e, tet]];

};

Intersect: ImplicitCubeTet.TetIntersectProc ~ {

DoIntersect: PROC [i1, i2, i3: INT] ~ {

Eq: PROC [s, t: Triple] RETURNS [b: BOOL] ~ {b Ź G3dVector.Distance[s, t] < epsilon};

IF i1 = i2 OR i1 = i3 OR i2 = i3

THEN NULL -- ERROR

ELSE {

a: Triple Ź vertices[i1].point;

b: Triple Ź vertices[i2].point;

c: Triple Ź vertices[i3].point;

IF G3dPolygon.TriangleArea[a, b, c] < minArea OR Eq[a, b] OR Eq[a, c] OR Eq[b, c]

THEN ImplicitDebug.Write[

IO.PutFR1["degenerate triangle (tetID: %g)\n", IO.int[tetID]]]

ELSE IF NOT action[i1, i2, i3, vertices] THEN ERROR Abort;

};

id0: INT Ź VIDfromEdge[e1];

id1: INT Ź VIDfromEdge[e2];

id2: INT Ź VIDfromEdge[e3];

DoIntersect[id0, id1, id2];

IF nEdges = 4 THEN DoIntersect[id0, id2, VIDfromEdge[e4]];

};

minArea: REAL Ź 0.2897354*epsilon*epsilon; -- constant determined from solid angle of tet

ImplicitDebug.SetDebug[[tet: tet]];

ImplicitCubeTet.IntersectWithSurface[tet, Intersect];

RETURN[vertices];

};

ShapeFromFunction: PUBLIC PROC [

function: Function,

level, size: REAL,

clientData: REF ANY Ź NIL,

start: Triple Ź [],

statusAction: TriangleProc Ź NIL]

RETURNS [s: G3dShape.Shape Ź NIL]

~ {

Tri: TriangleProc ~ {

triangle: G3dBasic.NatSequence Ź NEW[G3dBasic.NatSequenceRep[3]];

triangle[0] Ź i1; triangle[1] Ź i2; triangle[2] Ź i3;

triangle.length Ź 3;

triangles Ź G3dBasic.AddToSurfaceSequence[triangles, [NIL, triangle]];

points Ź vertices;

IF statusAction # NIL THEN RETURN[statusAction[i1, i2, i3, vertices]];

};

points: Vertices Ź NIL;

triangles: G3dBasic.SurfaceSequence Ź NIL;

startPt: Triple Ź IF start = [] THEN FindStart[[], function, size, level, clientData] ELSE start;

error: ROPE Ź Triangulate[function, level, size, clientData, startPt, Tri];

IF points = NIL OR points.length = 0 THEN

error Ź Rope.Cat[error, IF error # NIL THEN ", " ELSE NIL, "no vertices"];

IF triangles = NIL OR triangles.length = 0 THEN

error Ź Rope.Cat[error, IF error # NIL THEN ", " ELSE NIL, "no triangles"];

IF error # NIL

THEN Error[error]

ELSE s Ź NEW[G3dShape.ShapeRep Ź [vertices:points, surfaces:triangles, type:$Triangulated]];

};

Geometry

FindStart: PUBLIC PROC [

ballpark: Triple,

function: Function,

size, level: REAL,

clientData: REF ANY Ź NIL]

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.0) OR (sign = negative AND value < 0.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: PUBLIC PROC [

pos, neg: Triple, val, level: REAL, function: Function, clientData: REF ANY]

RETURNS [p: Triple]

~ {

p Ź G3dVector.Midpoint[pos, neg];

IF val < 0.0 THEN {c: Triple Ź neg; neg Ź pos; pos Ź c};

THROUGH [0..8) DO

IF (function[p, clientData]-level) > 0.0 THEN pos Ź p ELSE neg Ź p;

p Ź G3dVector.Midpoint[pos, neg];

ENDLOOP;

};

Normal: PUBLIC PROC [point: Triple, function: Function, delta: REAL, clientData: REF ANY]

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.