Track:
PUBLIC PROC [
seedSize: REAL,
surfacePoint: Triple,
valueProc: ValueProc,
threshold: REAL ¬ 1.0,
normalProc: NormalProc ¬ NIL,
statusProc: StatusProc ¬ NIL,
clientData: REF ANY ¬ NIL,
adaptLimit: NAT ¬ 0,
flatLimit: REAL ¬ 30.0,
tolerance: REAL ¬ 0.0001]
RETURNS [octree: Octree]
~ {
GetActiveFaces:
PROC [cube: Cube]
RETURNS [activeFaces: ActiveFaces] ~ {
FOR e: Edge
IN Edge
DO
corners: TwoCorners ¬ G3dOctree.EdgeCorners[cube, e];
IF corners.c0.inside # corners.c1.inside
AND
NOT corners.c0.outOfRange AND
NOT corners.c1.outOfRange
THEN {
twoFaces: TwoFaces ~ G3dOctree.EdgeFaces[e];
activeFaces[twoFaces.f0] ¬ activeFaces[twoFaces.f1] ¬ TRUE;
};
ENDLOOP;
};
ShouldSubdivideTest:
PROC [cube: Cube]
RETURNS [should:
BOOL ¬
FALSE] ~ {
crossPolygonProc: CrossPolygonProc ~ {
InnerTest:
PROC
RETURNS [shouldSubdivide:
BOOL ¬
FALSE] ~ {
SELECT
TRUE
FROM
nPolygon # 0 => RETURN[TRUE]; -- more than one polygon in cube
polygon.length >= 3 => {
Test polygon normal for collinearity with vertex normals:
center, normal, p1: Triple;
FOR n:
NAT
IN [0..polygon.length)
DO
points[n] ¬ polygon[n].point;
ENDLOOP;
points.length ¬ polygon.length;
center ¬ G3dPolygon.PolygonCenter[points];
normal ¬ ImplicitPoints.GetPointNormal[
center, valueProc, valueProc[center, clientData], delta, clientData];
FOR n:
NAT
IN [0..polygon.length)
DO
IF G3dVector.Dot[polygon[n].normal, normal] < cosFlatLimit
THEN RETURN[TRUE];
ENDLOOP;
Test polygon edges for perpendicularity to plane normal:
p1 ¬ polygon[polygon.length-1].point;
FOR n:
NAT
IN [0..polygon.length)
DO
p0: Triple ¬ p1;
v: Triple ¬ G3dVector.Sub[p1 ¬ polygon[n].point, p0];
v ¬ G3dVector.Unit[v];
IF G3dVector.Dot[normal, v] > cosPerpFlatLimit THEN RETURN[TRUE];
ENDLOOP;
};
ENDCASE;
};
IF (should ¬ InnerTest[]) THEN RETURN[FALSE];
};
IF TRUE THEN RETURN[FALSE];
IF 0.5*cube.size < sizeLimit THEN RETURN[FALSE];
IF cube.size # cubeSize THEN delta ¬ 0.001*(cubeSize ¬ cube.size);
[] ¬ ImplicitPoints.CheckAndGetCrossedEdges[
cube, valueProc, threshold, NProc, tolerance, clientData];
[] ¬ ImplicitPolygons.CrossPolygonizeCube[cube, crossPolygonProc, crossSequence, FALSE];
};
ShouldCoalesce:
PROC [parent: Cube]
RETURNS [should:
BOOL] ~ {
parent may've been coalesced by testing one of its kids, so first verify it's not terminal:
IF parent = NIL OR parent.terminal THEN RETURN[FALSE];
IF (should ¬
NOT ShouldSubdivideTest[parent])
THEN {
Check for edge-surface anomaly:
FOR edge: Edge
IN Edge
DO
twoOctants: TwoOctants ~ G3dOctree.EdgeOctants[edge];
state: BOOL ¬ parent.corners[twoOctants.o0].inside;
IF parent.corners[twoOctants.o1].inside = state
THEN {
o: Octant ¬ twoOctants.o1;
k: Cube ¬ parent.kids[twoOctants.o0];
IF k = NIL THEN {k ¬ parent.kids[twoOctants.o1]; o ¬ twoOctants.o0};
IF k # NIL AND k.corners[o].inside # state THEN RETURN[FALSE];
};
ENDLOOP;
Check for face-surface anomaly:
FOR face: Face
IN Face
DO
fourOctants: FourOctants ¬ G3dOctree.FaceOctants[face];
state: BOOL ¬ parent.corners[fourOctants.o0].inside;
IF parent.corners[fourOctants.o1].inside = state
AND
parent.corners[fourOctants.o2].inside = state AND
parent.corners[fourOctants.o3].inside = state
THEN {
o: Octant ¬ fourOctants.o2;
k: Cube ¬ parent.kids[fourOctants.o0];
IF k = NIL THEN {k ¬ parent.kids[fourOctants.o1]; o ¬ fourOctants.o3};
IF k = NIL THEN {k ¬ parent.kids[fourOctants.o2]; o ¬ fourOctants.o0};
IF k = NIL THEN {k ¬ parent.kids[fourOctants.o3]; o ¬ fourOctants.o1};
IF k # NIL AND k.corners[o].inside # state THEN RETURN[FALSE];
};
ENDLOOP;
};
};
DefaultNormalProc: NormalProc ~ {
RETURN[ImplicitPoints.GetPointNormal[point, valueProc, value, delta, clientData]];
};
AddToBottomOfStack:
PROC [cube: Cube] ~ {
IF (counter ¬ counter+1) > 10
THEN {
counter ¬ 0;
MessageWindow.Append[IO.PutFR1["\tstack size: %g", IO.int[cubeStack.size]], TRUE];
};
IF cube = NIL OR NOT cube.terminal OR cube.size < sizeLimit THEN RETURN;
ImplicitPoints.SetCornerValues[cube, valueProc, threshold, clientData];
IF NOT ImplicitPoints.IsCrossed[cube] THEN RETURN;
G3dOctree.WriteBottomOfCubeStack[cube, cubeStack];
cube.refAny ¬ $AddedToStack;
IF statusProc # NIL AND statusProc[cube] = $Abort THEN ERROR TrackAbort;
IF NOT G3dOctree.FullyPointed[G3dOctree.Root[cube]] THEN ERROR;
};
InnerSubdivide:
PROC [parent: Cube] ~ {
Subdivide parent; compute location, implicit value, and valueSet for all new child corners.
Subdivide edge neighbors if child midpoint introduces a new surface edge.
Subdivide face neighbor if child face center introduces a new surface region.
MaybeInheritCross:
PROC [edge: Edge] ~ {
kid: Corner ~ directionCorners[G3dOctree.DirectionFromEdge[edge]];
par: Corner ~ G3dOctree.EdgeCorners[parent, edge].c1;
IF par.inside = kid.inside
THEN
SELECT G3dOctree.AxisFromEdge[edge]
FROM
x => kid.lCross ¬ par.lCross;
y => kid.bCross ¬ par.bCross;
z => kid.nCross ¬ par.nCross;
ENDCASE;
};
MaybeSubdivideEdgeNeighbors:
PROC [edge: Edge] ~ {
twoCorners: TwoCorners ~ G3dOctree.EdgeCorners[parent, edge];
IF twoCorners.c0.inside = twoCorners.c1.inside
THEN {
kidCorner: Corner ~ directionCorners[G3dOctree.DirectionFromEdge[edge]];
IF twoCorners.c0.inside # kidCorner.inside
THEN {
ds: G3dOctree.ThreeDirections ¬ G3dOctree.EdgeDirections[edge];
c0: Cube ¬ neighbors[ds.d0];
c1: Cube ¬ neighbors[ds.d1];
c2: Cube ¬ neighbors[ds.d2];
IF c0 # NIL AND c0.terminal THEN InnerSubdivide[c0];
IF c1 # NIL AND c1.terminal THEN InnerSubdivide[c1];
IF c2 # NIL AND c2.terminal THEN InnerSubdivide[c2];
};
};
};
MaybeSubdivideFaceNeighbor:
PROC [face: Face] ~ {
d: FourDirections ¬ G3dOctree.FaceDirections[face];
bool: BOOL ~ directionCorners[d.d0].inside;
IF bool = directionCorners[d.d1].inside
AND
bool = directionCorners[d.d2].inside AND
bool = directionCorners[d.g3d].inside AND
bool # directionCorners[G3dOctree.DirectionFromFace[face]].inside
THEN {
neighbor: Cube ~ neighbors[G3dOctree.DirectionFromFace[face]];
IF neighbor # NIL AND neighbor.terminal THEN InnerSubdivide[neighbor];
};
};
directionCorners: DirectionCorners;
neighbors: Neighborhood ¬ G3dOctree.Neighbors[parent];
G3dOctree.Subdivide[parent];
ImplicitPoints.SetKidValues[parent, valueProc, threshold, clientData];
FOR o: Octant
IN Octant
DO
k: Cube ¬ parent.kids[o];
IF ShouldSubdivideTest[k]
THEN InnerSubdivide[k]
ELSE AddToBottomOfStack[k];
ENDLOOP;
directionCorners ¬ G3dOctree.GetDirectionCorners[parent];
FOR edge: Edge
IN Edge
DO
MaybeInheritCross[edge];
ENDLOOP;
FOR edge: Edge
IN Edge
DO
-- test for edge anomaly
MaybeSubdivideEdgeNeighbors[edge];
ENDLOOP;
FOR face: Face
IN Face
DO
-- test for face anomaly
MaybeSubdivideFaceNeighbor[face];
ENDLOOP;
};
ProcessTopOfStack:
PROC [cubeStack: CubeStack] ~ {
cube: Cube ¬ G3dOctree.ReadTopOfCubeStack[cubeStack];
IF (counter ¬ counter+1) > 10
THEN {
counter ¬ 0;
MessageWindow.Append[IO.PutFR1["\tstack size: %g", IO.int[cubeStack.size]], TRUE];
};
A cube may have been subdivided or coalesced since going on stack, so test:
IF cube.terminal
AND cube.refAny # $Coalesced
THEN {
activeFaces: ActiveFaces ¬ GetActiveFaces[cube];
FOR face: Face
IN Face
DO
In a prior pass through this loop, cube may get coalesced, so test each time:
IF cube.refAny = $Coalesced THEN RETURN;
IF activeFaces[face]
THEN {
neighbor: Cube ¬ G3dOctree.FaceNeighbor[cube, face];
IF neighbor =
NIL
THEN {
new: Cube;
New cube's parent must be fully pointed since RemoveCoalescedCubes and
ImplicitPolygons.RecursiveTrackFace presume nonterminal cubes have all kids
octree.root ¬ G3dOctree.AddCube[cube, face];
new ¬ G3dOctree.FaceNeighbor[cube, face];
IF NOT new.terminal THEN ERROR;
FOR p: Cube ¬ new.parent, p.parent
WHILE p #
NIL
DO
If new attached to coalesced part of octree, ensure all branch is coalesced;
leave new in octree (reduce future neighbor testing), don't add to stack:
IF p.refAny = $Coalesced
THEN {
Coalesce[p];
RETURN;
};
ENDLOOP;
ConsiderCube[new];
}
If neighbor was created as artifact of prior AddCube, it needs to be considered:
ELSE IF neighbor.refAny = NIL THEN ConsiderCube[neighbor];
};
ENDLOOP;
};
ELSE cube.refAny ¬ $RejectedFromStack;
};
Coalesce:
PROC [parent: Cube] ~ {
Coalesce parent.kids but keep .kids so kids' neighbors needn't re-investigate their space:
Inner: CubeProc ~ {cube.refAny ¬ $Coalesced};
G3dOctree.ApplyToKids[parent, Inner];
};
ConsiderCube:
PROC [new: Cube] ~ {
IF new.refAny = $Coalesced THEN ERROR;
new.refAny ¬ $Considered;
IF ShouldSubdivideTest[new]
THEN InnerSubdivide[new]
ELSE
IF ShouldCoalesce[new.parent]
THEN {
Coalesce[new.parent];
Parent is now effectively terminal, so, in order to be processed from stack:
new.parent.terminal ¬ TRUE;
AddToBottomOfStack[new.parent];
}
ELSE AddToBottomOfStack[new];
};
RemoveCoalescedCubes:
PROC [root: Cube] ~ {
IsCoalesced:
PROC [cube: Cube]
RETURNS [b:
BOOL ¬
TRUE] ~ {
Test: CubeProc ~ {IF cube.refAny # $Coalesced THEN ERROR};
IF (b ¬ cube.refAny = $Coalesced) THEN G3dOctree.ApplyToKids[cube, Test];
};
Inner:
PROC [cube: Cube] ~ {
cube.terminal ¬ TRUE;
FOR o: Octant
IN Octant
DO
k: Cube ¬ cube.kids[o];
IF k = NIL THEN LOOP;
IF IsCoalesced[k]
THEN cube.kids[o] ¬ NIL
ELSE {
cube.terminal ¬ FALSE;
Inner[k];
};
ENDLOOP;
};
Inner[root];
};
CheckForNeededSubdivisions:
PROC [root: Cube] ~ {
Inner: CubeProc ~ {
state: BOOL;
neighbors: Neighborhood ¬ G3dOctree.Neighbors[cube];
TestNeighborIfDivide:
PROC [d: Direction]
RETURNS [b:
BOOL] ~ {
n: Cube ¬ neighbors[d];
b ¬ n # NIL AND NOT n.terminal AND n.refAny # $Coalesced;
IF b
THEN {
c: Corner ¬ G3dOctree.GetDirectionCorner[n, G3dOctree.OppositeDirection[d]];
b ¬ c # NIL AND state # c.inside;
};
};
TestEdgeNeighborIfDivide:
PROC [d0, d1: Direction]
RETURNS [b:
BOOL] ~ {
n: Cube ¬ neighbors[d0];
b ¬ n # NIL AND NOT n.terminal AND n.refAny # $Coalesced;
IF b
THEN {
d: Direction ¬ G3dOctree.AddDirection[d1, G3dOctree.OppositeDirection[d0]];
c: Corner ¬ G3dOctree.GetDirectionCorner[n, d];
b ¬ c # NIL AND state # c.inside;
};
};
FOR edge: Edge
IN Edge
DO
twoCorners: TwoCorners ~ G3dOctree.EdgeCorners[cube, edge];
state ¬ twoCorners.c0.inside;
IF state = twoCorners.c1.inside
THEN {
ds: G3dOctree.ThreeDirections ¬ G3dOctree.EdgeDirections[edge];
IF TestEdgeNeighborIfDivide[ds.d0, ds.d1] THEN GOTO subdiv;
IF TestEdgeNeighborIfDivide[ds.d1, ds.d0] THEN GOTO subdiv;
IF TestNeighborIfDivide[ds.d2] THEN GOTO subdiv;
};
ENDLOOP;
FOR face: Face
IN Face
DO
fourCorners: FourCorners ¬ G3dOctree.FaceCorners[cube, face];
state ¬ fourCorners.c0.inside;
IF state = fourCorners.c1.inside
AND
state = fourCorners.c2.inside AND
state = fourCorners.c3.inside AND
TestNeighborIfDivide[G3dOctree.DirectionFromFace[face]] THEN GOTO subdiv;
ENDLOOP;
EXITS
subdiv => {
InnerSubdivide[cube];
FOR o: Octant IN Octant DO [] ¬ Inner[cube.kids[o]]; ENDLOOP;
};
};
G3dOctree.ApplyToTerminal[root, Inner];
};
NProc: NormalProc ¬ IF normalProc # NIL THEN normalProc ELSE DefaultNormalProc;
counter, nRecursiveTrackFaces, nRecurseFailures: NAT ¬ 0;
cubeStack: CubeStack ¬ G3dOctree.NewCubeStack[20000];
cosFlatLimit: REAL ¬ RealFns.CosDeg[flatLimit];
cosPerpFlatLimit: REAL ¬ RealFns.CosDeg[90.0-flatLimit];
delta, cubeSize: REAL ¬ 0.0;
crossSequence: CrossSequence ¬ G3dOctree.ObtainCrossSequence[];
points: TripleSequence ¬ NEW[TripleSequenceRep[12]];
sizeLimit: REAL ¬ seedSize;
FOR n: NAT IN [0..adaptLimit) DO sizeLimit ¬ 0.5*sizeLimit; ENDLOOP;
octree ¬ NEW[OctreeRep];
octree.root ¬ G3dOctree.NewCube[seedSize, surfacePoint];
ConsiderCube[octree.root];
WHILE
NOT G3dOctree.CubeStackEmpty[cubeStack]
DO
ProcessTopOfStack[cubeStack ! TrackAbort => EXIT];
ENDLOOP;
CheckForNeededSubdivisions[octree.root];
RemoveCoalescedCubes[octree.root];
IF NOT G3dOctree.FullyPointed[octree.root] THEN ERROR;
IF NOT G3dOctree.CubeOk[octree.root] THEN ERROR;
G3dOctree.ReleaseCrossSequence[crossSequence];
G3dOctree.SetOctreeFields[octree];
};