<<>>
<<G3dOctreeImpl.mesa>>
<<Copyright Ó 1985, 1992 by Xerox Corporation.  All rights reserved.>>
<<Bloomenthal, November 21, 1992 2:30 pm PST>>

DIRECTORY Basics, G3dBasic, G3dBox, G3dIO, G3dOctree, G3dPlane, FS, FileNames, IO, Real, RefText, Rope;

G3dOctreeImpl: CEDAR MONITOR
IMPORTS Basics, G3dBox, G3dIO, G3dOctree, G3dPlane, FS, FileNames, IO, Real, RefText, Rope
EXPORTS G3dOctree

~ BEGIN
<<Types>>
ROPE:                        TYPE ~ Rope.ROPE;
STREAM:                    TYPE ~ IO.STREAM;
Box:                        TYPE ~ G3dBasic.Box;
Pair:                        TYPE ~ G3dBasic.Pair;
Triple:                     TYPE ~ G3dBasic.Triple;
Plane:                        TYPE ~ G3dPlane.Plane;

ActiveFaces:                TYPE ~ G3dOctree.ActiveFaces;
Axis:                        TYPE ~ G3dOctree.Axis;
Corner:                    TYPE ~ G3dOctree.Corner;
CornerDataProc:            TYPE ~ G3dOctree.CornerDataProc;
CornerProc:                TYPE ~ G3dOctree.CornerProc;
CornerRep:                TYPE ~ G3dOctree.CornerRep;
Corners:                    TYPE ~ G3dOctree.Corners;
Cross:                        TYPE ~ G3dOctree.Cross;
CrossArray:                TYPE ~ G3dOctree.CrossArray;
CrossPolygonProc:        TYPE ~ G3dOctree.CrossPolygonProc;
CrossRep:                    TYPE ~ G3dOctree.CrossRep;
CrossSequence:            TYPE ~ G3dOctree.CrossSequence;
CrossSequenceRep:        TYPE ~ G3dOctree.CrossSequenceRep;
CrossedEdge:                TYPE ~ G3dOctree.CrossedEdge;
CrossedEdges:            TYPE ~ G3dOctree.CrossedEdges;
Cube:                        TYPE ~ G3dOctree.Cube;
CubeObject:                TYPE ~ G3dOctree.CubeObject;
CubeProc:                TYPE ~ G3dOctree.CubeProc;
CubeRep:                    TYPE ~ G3dOctree.CubeRep;
CubeStack:                TYPE ~ G3dOctree.CubeStack;
CubeStackRep:            TYPE ~ G3dOctree.CubeStackRep;
CubeSequence:            TYPE ~ G3dOctree.CubeSequence;
CubeSequenceRep:        TYPE ~ G3dOctree.CubeSequenceRep;
Direction:                    TYPE ~ G3dOctree.Direction;
DirectionCorners:        TYPE ~ G3dOctree.DirectionCorners;
DirectionPairs:            TYPE ~ G3dOctree.DirectionPairs;
DirectionSelect:            TYPE ~ G3dOctree.DirectionSelect;
DirectionType:            TYPE ~ G3dOctree.DirectionType;
Edge:                        TYPE ~ G3dOctree.Edge;
Face:                        TYPE ~ G3dOctree.Face;
FaceNeighborInfo:        TYPE ~ G3dOctree.FaceNeighborInfo;
Intersection:                TYPE ~ G3dOctree.Intersection;
IntersectionList:            TYPE ~ G3dOctree.IntersectionList;
IntersectionType:        TYPE ~ G3dOctree.IntersectionType;
Kids:                        TYPE ~ G3dOctree.Kids;
Neighborhood:            TYPE ~ G3dOctree.Neighborhood;
Octant:                    TYPE ~ G3dOctree.Octant;
OctantPairs:                TYPE ~ G3dOctree.OctantPairs;
Octree:                    TYPE ~ G3dOctree.Octree;
OctreeMode:                TYPE ~ G3dOctree.OctreeMode;
OctreeRep:                TYPE ~ G3dOctree.OctreeRep;
OctreeType:                TYPE ~ G3dOctree.OctreeType;
FourDirections:            TYPE ~ G3dOctree.FourDirections;
FourOctants:                TYPE ~ G3dOctree.FourOctants;
FourCorners:                TYPE ~ G3dOctree.FourCorners;
ThreeCubes:                TYPE ~ G3dOctree.ThreeCubes;
ThreeDirections:            TYPE ~ G3dOctree.ThreeDirections;
ThreeEdges:                TYPE ~ G3dOctree.ThreeEdges;
TwoCorners:                TYPE ~ G3dOctree.TwoCorners;
TwoCubes:                TYPE ~ G3dOctree.TwoCubes;
TwoDirections:            TYPE ~ G3dOctree.TwoDirections;
TwoFaces:                TYPE ~ G3dOctree.TwoFaces;
TwoOctants:                TYPE ~ G3dOctree.TwoOctants;
<<Allocation>>
NewCube: PUBLIC PROC [
size: REAL,
center: Triple Ź [],
cornerDataProc: CornerDataProc Ź NIL]
RETURNS [Cube]
~ {
RETURN[NewKid[NIL, size, lbn, 0, NewCorners[size, center, cornerDataProc]]];
};

NewKid: PROC [
parent: Cube,
size: REAL,
octant: Octant,
level: INT,
corners: Corners]
RETURNS [cube: Cube]
~ {
cube Ź NEW[CubeRep];
cube.parent Ź parent;
cube.size Ź size;
cube.octant Ź octant;
cube.level Ź level;
cube.terminal Ź TRUE;
cube.corners Ź corners;
};

NewCorner: PROC [
point: Triple,
cornerDataProc: CornerDataProc Ź NIL]
RETURNS [corner: Corner]
~ {
corner Ź NEW[CornerRep Ź [point: point]];
};

NewCorners: PUBLIC PROC [
size: REAL,
center: Triple,
cornerDataProc: CornerDataProc Ź NIL]
RETURNS [corners: Corners]
~ {
d: REAL ~ 0.5*size;
l: REAL ~ center.x-d;
r: REAL ~ center.x+d;
b: REAL ~ center.y-d;
t: REAL ~ center.y+d;
n: REAL ~ center.z-d;
f: REAL ~ center.z+d;
corners Ź [
NewCorner[[l, b, n], cornerDataProc],
NewCorner[[l, b, f], cornerDataProc],
NewCorner[[l, t, n], cornerDataProc],
NewCorner[[l, t, f], cornerDataProc],
NewCorner[[r, b, n], cornerDataProc],
NewCorner[[r, b, f], cornerDataProc],
NewCorner[[r, t, n], cornerDataProc],
NewCorner[[r, t, f], cornerDataProc]
];
};
<<Attributes>>
Root: PUBLIC PROC [cube: Cube] RETURNS [root: Cube] ~ {
root Ź cube;
IF cube # NIL THEN WHILE root.parent # NIL DO root Ź root.parent; ENDLOOP;
};
<<>>
CubeOk: PUBLIC PROC [cube: Cube] RETURNS [BOOL] ~ {
k: BOOL Ź FALSE;
IF cube = NIL THEN RETURN[FALSE];
FOR o: Octant IN Octant DO IF cube.kids[o] # NIL THEN {k Ź TRUE; EXIT}; ENDLOOP;
IF k = cube.terminal THEN RETURN[FALSE];
FOR o: Octant IN Octant DO
k: Cube Ź cube.kids[o];
IF k # NIL THEN {
IF k.parent # cube THEN RETURN[FALSE];
IF NOT CubeOk[k] THEN RETURN[FALSE];
};
ENDLOOP;
RETURN[TRUE];
};

FullyPointed: PUBLIC PROC [root: Cube] RETURNS [b: BOOL Ź TRUE] ~ {
bad: ERROR = CODE;
Inner: PROC [cube: Cube] ~ {
terminal: BOOL;
FOR o: Octant IN Octant DO
k: Cube Ź cube.kids[o];
IF o = Octant.FIRST
THEN terminal Ź k = NIL
ELSE IF (k # NIL) AND terminal THEN ERROR bad;
IF k # NIL THEN Inner[k];
ENDLOOP;
};
Inner[root ! bad => {b Ź FALSE; CONTINUE}];
};

NCubes: PUBLIC PROC [root: Cube] RETURNS [nCubes: INT Ź 0] ~ {
cubeProc: CubeProc ~ {nCubes Ź nCubes+1};
Apply[root, cubeProc];
};

NTerminalCubes: PUBLIC PROC [root: Cube] RETURNS [nCubes: INT Ź 0] ~ {
cubeProc: CubeProc ~ {nCubes Ź nCubes+1};
ApplyToTerminal[root, cubeProc];
};

Center: PUBLIC PROC [cube: Cube] RETURNS [Triple] ~ {
IF cube # NIL THEN {
p0: Triple ~ cube.corners[lbn].point;
p1: Triple ~ cube.corners[rtf].point;
RETURN[[0.5*(p0.x+p1.x), 0.5*(p0.y+p1.y), 0.5*(p0.z+p1.z)]];
};
RETURN[[0.0, 0.0, 0.0]];
};

Size: PUBLIC PROC [cube: Cube] RETURNS [REAL] ~ {
RETURN[IF cube = NIL THEN 0.0 ELSE cube.size];
};

MinSize: PUBLIC PROC [root: Cube] RETURNS [min: REAL] ~ {
MinProc: CubeProc ~ {IF cube.size < min THEN min Ź cube.size};
min Ź 100000000.0;
ApplyToTerminal[root, MinProc];
};

BoxOfCube: PUBLIC PROC [cube: Cube] RETURNS [Box] ~ {
RETURN[[cube.corners[lbn].point, cube.corners[rtf].point]];
};

AnyKids: PUBLIC PROC [cube: Cube] RETURNS [BOOL] ~ {
FOR o: Octant IN Octant DO IF cube.kids[o] # NIL THEN RETURN[TRUE]; ENDLOOP;
RETURN[FALSE];
};

DepthOf: PUBLIC PROC [cube: Cube] RETURNS [d: INT Ź 0] ~ {
Inner: PROC [cube: Cube, level: INT] ~ {
ll: INT Ź level+1;
FOR o: Octant IN Octant DO
k: Cube Ź cube.kids[o];
IF k # NIL THEN Inner[k, ll];
REPEAT
FINISHED => IF level > d THEN d Ź level;
ENDLOOP;
};
Inner[cube, 0];
};
<<Subdivision>>
Subdivide: PUBLIC PROC [cube: Cube] ~ {
IF cube # NIL AND cube.terminal THEN {
FaceTry: PROC [d: Direction] RETURNS [corner: Corner Ź NIL] ~ {
KidsTry: PROC [o0a, o0b, o1a, o1b: Octant] RETURNS [corner: Corner Ź NIL] ~ {
IF c.kids[o0a] # NIL THEN RETURN[c.kids[o0a].corners[o0b]];
IF c.kids[o0b] # NIL THEN RETURN[c.kids[o0b].corners[o0a]];
IF c.kids[o1a] # NIL THEN RETURN[c.kids[o1a].corners[o1b]];
IF c.kids[o1b] # NIL THEN RETURN[c.kids[o1b].corners[o1a]];
};
c: Cube ~ neighbors[d];
IF c # NIL AND NOT c.terminal
THEN corner Ź SELECT d FROM
l    => KidsTry[rbn, rtf, rbf, rtn],
r    => KidsTry[lbn, ltf, lbf, ltn],
b    => KidsTry[ltn, rtf, ltf, rtn],
t    => KidsTry[lbn, rbf, lbf, rbn],
n    => KidsTry[lbf, rtf, ltf, rbf],
f    => KidsTry[lbn, rtn, ltn, rbn]
ENDCASE => NIL;
};
EdgeTry: PROC [d0, d1, d2: Direction, e0, e1, e2: Edge] RETURNS [corner: Corner] ~ {
NeighborTry: PROC [d: Direction, e: Edge] RETURNS [corner: Corner Ź NIL] ~ {
KidsTry: PROC [o0, o1: Octant] RETURNS [corner: Corner Ź NIL] ~ {
IF c.kids[o0] # NIL THEN RETURN[c.kids[o0].corners[o1]];
IF c.kids[o1] # NIL THEN RETURN[c.kids[o1].corners[o0]];
};
c: Cube ~ neighbors[d];
IF c # NIL AND NOT c.terminal
THEN corner Ź SELECT e FROM
lb    => KidsTry[lbn, lbf],
lt    => KidsTry[ltn, ltf],
ln    => KidsTry[lbn, ltn],
lf    => KidsTry[lbf, ltf],
rb    => KidsTry[rbn, rbf],
rt    => KidsTry[rtn, rtf],
rn    => KidsTry[rbn, rtn],
rf    => KidsTry[rbf, rtf],
bn    => KidsTry[lbn, rbn],
bf    => KidsTry[lbf, rbf],
tn    => KidsTry[ltn, rtn],
tf    => KidsTry[ltf, rtf],
ENDCASE => NIL;
};
corner Ź NeighborTry[d0, e0];
IF corner = NIL THEN corner Ź NeighborTry[d1, e1];
IF corner = NIL THEN corner Ź NeighborTry[d2, e2];
};
SubCube: PROC [octant: Octant, corners: Corners] ~ {
kid: Cube ~ cube.kids[octant] Ź NewKid[cube, kSize, octant, kidLevel, corners];
kid.terminal Ź TRUE;
kid.size Ź kSize;
};
kSize: REAL ~ 0.5*Size[cube];
kidLevel: INT ~ cube.level+1;
c: Triple ~ Center[cube];
neighbors: Neighborhood ~ Neighbors[cube];
cornerDataProc: CornerDataProc Ź NIL;
ccc:    Corner ~ NewCorner[c, cornerDataProc];
lbn:    Corner ~ cube.corners[lbn];
lbf:    Corner ~ cube.corners[lbf];
ltn:    Corner ~ cube.corners[ltn];
ltf:    Corner ~ cube.corners[ltf];
rbn:    Corner ~ cube.corners[rbn];
rbf:    Corner ~ cube.corners[rbf];
rtn:    Corner ~ cube.corners[rtn];
rtf:    Corner ~ cube.corners[rtf];
lb:        Corner Ź EdgeTry[l, lb, b,        rb, rt, lt];
lt:        Corner Ź EdgeTry[l, lt, t,        rt, rb, lb];
ln:        Corner Ź EdgeTry[l, ln, n,    rn, rf, lf];
lf:        Corner Ź EdgeTry[l, lf, f,        rf, rn, ln];
rb:        Corner Ź EdgeTry[r, rb, b,    lb, lt, rt];
rt:        Corner Ź EdgeTry[r, rt, t,        lt, lb, rb];
rn:    Corner Ź EdgeTry[r, rn, n,    ln, lf, rf];
rf:        Corner Ź EdgeTry[r, rf, f,        lf, ln, rn];
bn:    Corner Ź EdgeTry[b, bn, n,    tn, tf, bf];
bf:        Corner Ź EdgeTry[b, bf, f,    tf, tn, bn];
tn:        Corner Ź EdgeTry[t, tn, n,    bn, bf, tf];
tf:        Corner Ź EdgeTry[t, tf, f,        bf, bn, tn];
l:        Corner Ź FaceTry[l];
r:        Corner Ź FaceTry[r];
b:        Corner Ź FaceTry[b];
t:        Corner Ź FaceTry[t];
n:        Corner Ź FaceTry[n];
f:        Corner Ź FaceTry[f];
IF lb = NIL THEN lb Ź    NewCorner[[c.x-kSize,    c.y-kSize,    c.z], cornerDataProc];
IF lt = NIL THEN lt Ź    NewCorner[[c.x-kSize,    c.y+kSize,    c.z], cornerDataProc];
IF ln = NIL THEN ln Ź    NewCorner[[c.x-kSize,    c.y,            c.z-kSize], cornerDataProc];
IF lf = NIL THEN lf Ź    NewCorner[[c.x-kSize,    c.y,            c.z+kSize], cornerDataProc];
IF rb = NIL THEN rb Ź    NewCorner[[c.x+kSize,    c.y-kSize,    c.z], cornerDataProc];
IF rt = NIL THEN rt Ź    NewCorner[[c.x+kSize,    c.y+kSize,    c.z], cornerDataProc];
IF rn = NIL THEN rn Ź    NewCorner[[c.x+kSize,    c.y,            c.z-kSize], cornerDataProc];
IF rf = NIL THEN rf Ź    NewCorner[[c.x+kSize,    c.y,            c.z+kSize], cornerDataProc];
IF bn = NIL THEN bn Ź NewCorner[[c.x,            c.y-kSize,    c.z-kSize], cornerDataProc];
IF bf = NIL THEN bf Ź    NewCorner[[c.x,            c.y-kSize,    c.z+kSize], cornerDataProc];
IF tn = NIL THEN tn Ź    NewCorner[[c.x,            c.y+kSize,    c.z-kSize], cornerDataProc];
IF tf = NIL THEN tf Ź    NewCorner[[c.x,            c.y+kSize,    c.z+kSize], cornerDataProc];
IF l = NIL THEN    l Ź    NewCorner[[c.x-kSize,    c.y,            c.z], cornerDataProc];
IF r = NIL THEN    r Ź    NewCorner[[c.x+kSize,    c.y,            c.z], cornerDataProc];
IF b = NIL THEN b Ź    NewCorner[[c.x,            c.y-kSize,    c.z], cornerDataProc];
IF t = NIL THEN t Ź        NewCorner[[c.x,            c.y+kSize,    c.z], cornerDataProc];
IF n = NIL THEN n Ź    NewCorner[[c.x,            c.y,            c.z-kSize], cornerDataProc];
IF f = NIL THEN f Ź    NewCorner[[c.x,            c.y,            c.z+kSize], cornerDataProc];
cube.terminal Ź FALSE;
SubCube[lbn, [lbn, lb, ln, l, bn, b, n, ccc]];
SubCube[lbf, [lb, lbf, l, lf, b, bf, ccc, f]];
SubCube[ltn, [ln, l, ltn, lt, n, ccc, tn, t]];
SubCube[ltf, [l, lf, lt, ltf, ccc, f, t, tf]];
SubCube[rbn, [bn, b, n, ccc, rbn, rb, rn, r]];
SubCube[rbf, [b, bf, ccc, f, rb, rbf, r, rf]];
SubCube[rtn, [n, ccc, tn, t, rn, r, rtn, rt]];
SubCube[rtf, [ccc, f, t, tf, r, rf, rt, rtf]];
};
};

SubdivideTerminal: PUBLIC PROC [cube: Cube] ~ {
cubeProc: CubeProc ~ {Subdivide[cube]};
ApplyToTerminal[cube, cubeProc];
};

RecursivelySubdivide: PUBLIC PROC [cube: Cube, nLevels: INT] ~ {
THROUGH [0..nLevels) DO SubdivideTerminal[cube]; ENDLOOP;
};
<<Modification>>
TestCube: PROC [cube: Cube] ~ {
tol: REAL ~ 0.01*Size[cube];
InnerTest: PROC [cube: Cube] ~ {
Test: PROC [o0, o1: Octant, s: {x, y, z}] ~ {
IF s # x AND ABS[c[o0].point.x-c[o1].point.x] > tol THEN ERROR;
IF s # y AND ABS[c[o0].point.y-c[o1].point.y] > tol THEN ERROR;
IF s # z AND ABS[c[o0].point.z-c[o1].point.z] > tol THEN ERROR;
};
c: Corners Ź cube.corners;
Test[lbn, ltn, y]; Test[lbf, ltf, y]; Test[lbn, lbf, z]; Test[ltn, ltf, z];
Test[rbn, rtn, y]; Test[rbf, rtf, y]; Test[rbn, rbf, z]; Test[rtn, rtf, z];
Test[lbn, rbn, x]; Test[lbf, rbf, x]; Test[ltn, rtn, x]; Test[ltf, rtf, x];
};
InnerTest[cube];
};

AddCube: PUBLIC PROC [cube: Cube, face: Face, fullyPointed: BOOL Ź TRUE]
RETURNS [root: Cube]
~ {
IF cube = NIL THEN RETURN[NIL];
TestCube[cube];
IF FaceNeighbor[cube, face] # NIL THEN RETURN[Root[cube]];
IF cube.parent = NIL THEN {
<<Create cube's parent, using cube's corner[o] as the center of parent:>>
cubeProc: CubeProc ~ {cube.level Ź cube.level+1};
NewParent: PUBLIC PROC [cube: Cube, o: Octant] RETURNS [root: Cube] ~ {
oo: Octant ~ G3dOctree.OppositeOctant[o];
d: REAL ~ 2.0*Size[cube];
l, r, b, t, n, f: REAL;
c: Triple ~ cube.corners[o].point;
pp: Triple ~ cube.corners[oo].point;
SELECT oo FROM
lbn        => {l Ź pp.x;    b Ź pp.y;    n Ź pp.z;    r Ź pp.x+d;    t Ź pp.y+d;    f Ź pp.z+d};
lbf        => {l Ź pp.x;    b Ź pp.y;    f Ź pp.z;    r Ź pp.x+d;    t Ź pp.y+d;    n Ź pp.z-d};
ltn        => {l Ź pp.x;    t Ź pp.y;    n Ź pp.z;    r Ź pp.x+d;    b Ź pp.y-d;    f Ź pp.z+d};
ltf            => {l Ź pp.x;    t Ź pp.y;    f Ź pp.z;    r Ź pp.x+d;    b Ź pp.y-d;    n Ź pp.z-d};
rbn        => {r Ź pp.x;    b Ź pp.y;    n Ź pp.z;    l Ź pp.x-d;    t Ź pp.y+d;    f Ź pp.z+d}; 
rbf        => {r Ź pp.x;    b Ź pp.y;    f Ź pp.z;    l Ź pp.x-d;    t Ź pp.y+d;    n Ź pp.z-d};
rtn        => {r Ź pp.x;    t Ź pp.y;    n Ź pp.z;    l Ź pp.x-d;    b Ź pp.y-d;    f Ź pp.z+d};
ENDCASE    => {r Ź pp.x;    t Ź pp.y;    f Ź pp.z;    l Ź pp.x-d;    b Ź pp.y-d;    n Ź pp.z-d};
root Ź NEW[CubeRep Ź [size: d]];
root.corners[oo] Ź cube.corners[oo];
IF oo # lbn THEN root.corners[lbn] Ź NewCorner[[l, b, n]];
IF oo # lbf THEN root.corners[lbf] Ź NewCorner[[l, b, f]];
IF oo # ltn THEN root.corners[ltn] Ź NewCorner[[l, t, n]];
IF oo # ltf THEN root.corners[ltf] Ź NewCorner[[l, t, f]];
IF oo # rbn THEN root.corners[rbn] Ź NewCorner[[r, b, n]];
IF oo # rbf THEN root.corners[rbf] Ź NewCorner[[r, b, f]];
IF oo # rtn THEN root.corners[rtn] Ź NewCorner[[r, t, n]];
IF oo # rtf THEN root.corners[rtf] Ź NewCorner[[r, t, f]];
};
o: Octant ~ SELECT face FROM l=> lbn, r=> rtf, b=> rbn, t=> ltf, n=> ltn, ENDCASE=> rbf;
opp: Octant ~ G3dOctree.OppositeOctant[o];
cube.parent Ź root Ź NewParent[cube, o];
cube.octant Ź opp;
root.corners[opp] Ź cube.corners[opp];
<<Notice that the root has a kid, but is not itself truly subdivided:>>
root.kids[opp] Ź cube;
Apply[root, cubeProc];
root.level Ź 0;
TestCube[root];
};
{
AddCubeToParent: PROC [parent: Cube, octant: Octant, size: REAL, level: INT] ~ {
new: CubeŹ NEW[CubeRep Ź [size: size, octant: octant, terminal: TRUE, parent: parent]];
c: Corner Ź parent.corners[octant];
l, r: REAL Ź parent.corners[octant].point.x;
b, t: REAL Ź parent.corners[octant].point.y;
n, f: REAL Ź parent.corners[octant].point.z;
SELECT octant FROM
lbn    => {r Ź l+size;    t Ź b+size;    f Ź n+size};
lbf    => {r Ź l+size;    t Ź b+size;    n Ź f-size};
ltn    => {r Ź l+size;    b Ź t-size;    f Ź n+size};
ltf        => {r Ź l+size;    b Ź t-size;    n Ź f-size};
rbn    => {l Ź r-size;    t Ź b+size;    f Ź n+size};
rbf    => {l Ź r-size;    t Ź b+size;    n Ź f-size};
rtn    => {l Ź r-size;    b Ź t-size;    f Ź n+size};
rtf        => {l Ź r-size;    b Ź t-size;    n Ź f-size};
ENDCASE;
parent.kids[octant] Ź new;
new.level Ź level;
FOR face: Face IN Face DO
neighbor: Cube ~ FaceNeighbor[new, face];
IF neighbor # NIL THEN {
n: Corners ~ neighbor.corners;
SELECT face FROM
l => {
new.corners[lbn] Ź n[rbn]; new.corners[lbf] Ź n[rbf];
new.corners[ltn] Ź n[rtn];    new.corners[ltf] Ź n[rtf];
};
r => {
new.corners[rbn] Ź n[lbn]; new.corners[rbf] Ź n[lbf];
new.corners[rtn] Ź n[ltn];    new.corners[rtf] Ź n[ltf];
};
b => {
new.corners[lbn] Ź n[ltn];    new.corners[lbf] Ź n[ltf];
new.corners[rbn] Ź n[rtn];    new.corners[rbf] Ź n[rtf];
};
t => {
new.corners[ltn] Ź n[lbn];    new.corners[ltf] Ź n[lbf];
new.corners[rtn] Ź n[rbn];    new.corners[rtf] Ź n[rbf];
};
n => {
new.corners[lbn] Ź n[lbf];    new.corners[ltn] Ź n[ltf];
new.corners[rbn] Ź n[rbf];    new.corners[rtn] Ź n[rtf];
};
ENDCASE => {
new.corners[lbf] Ź n[lbn];    new.corners[ltf] Ź n[ltn];
new.corners[rbf] Ź n[rbn];    new.corners[rtf] Ź n[rtn];
};
};
ENDLOOP;
new.corners[octant] Ź parent.corners[octant];
FOR o: Octant IN Octant DO
IF new.corners[o] = NIL THEN {
new.corners[o] Ź NEW[CornerRep];
new.corners[o].point Ź SELECT o FROM
lbn        => [l, b, n],
lbf        => [l, b, f],
ltn        => [l, t, n],
ltf            => [l, t, f],
rbn        => [r, b, n],
rbf        => [r, b, f],
rtn        => [r, t, n],
ENDCASE    => [r, t, f];
};
ENDLOOP;
TestCube[new];
};
parent: Cube Ź cube.parent;
info: FaceNeighborInfo ~ GetFaceNeighborInfo[cube, face];
IF info.recurse THEN {
<<This may change cube's level:>>
neighborOfParent: Cube Ź FaceNeighbor[cube.parent, face];
IF neighborOfParent = NIL THEN root Ź AddCube[cube.parent, face, fullyPointed];
IF (neighborOfParent Ź FaceNeighbor[cube.parent, face]) = NIL THEN ERROR;
parent Ź neighborOfParent;
};
parent.terminal Ź FALSE;
AddCubeToParent[parent, info.nOctant, cube.size, cube.level];
IF fullyPointed THEN FOR o: Octant IN Octant DO
IF parent.kids[o] = NIL THEN AddCubeToParent[parent, o, cube.size, cube.level];
ENDLOOP;
};
RETURN[Root[cube]];
};

DeleteCube: PUBLIC PROC [cube: Cube] ~ {
IF cube # NIL AND cube.parent # NIL THEN cube.parent.kids[cube.octant] Ź NIL;
};

MakeRestrictedOctree: PUBLIC PROC [root: Cube] ~ {
Inner: CubeProc ~ {
Test: PROC [c: Cube] RETURNS [b: BOOL] ~ {b Ź c # NIL AND NOT c.terminal};
IF cube # NIL THEN {
neighbors: Neighborhood Ź Neighbors[cube];
FOR face: Face IN Face DO
neighbor: Cube Ź neighbors[DirectionFromFace[face]];
IF neighbor # NIL AND NOT neighbor.terminal THEN {
fourOcts: FourOctants Ź FaceOctants[G3dOctree.OppositeFace[face]];
IF Test[neighbor.kids[fourOcts.o0]] THEN GOTO subdivide;
IF Test[neighbor.kids[fourOcts.o1]] THEN GOTO subdivide;
IF Test[neighbor.kids[fourOcts.o2]] THEN GOTO subdivide;
IF Test[neighbor.kids[fourOcts.o3]] THEN GOTO subdivide;
};
ENDLOOP;
FOR edge: Edge IN Edge DO
neighbor: Cube Ź neighbors[DirectionFromEdge[edge]];
IF neighbor # NIL AND NOT neighbor.terminal THEN {
twoOcts: TwoOctants Ź EdgeOctants[DiagonalEdge[edge]];
IF Test[neighbor.kids[twoOcts.o0]] THEN GOTO subdivide;
IF Test[neighbor.kids[twoOcts.o1]] THEN GOTO subdivide;
};
ENDLOOP;
FOR octant: Octant IN Octant DO
neighbor: Cube Ź neighbors[DirectionFromOctant[octant]];
IF neighbor # NIL AND NOT neighbor.terminal
AND Test[neighbor.kids[G3dOctree.OppositeOctant[octant]]] THEN GOTO subdivide;
ENDLOOP;
};
EXITS
subdivide => {
Subdivide[cube];
IF cube.parent # NIL THEN {
grandNeighbors: Neighborhood Ź Neighbors[cube.parent];
threeEdges: ThreeEdges Ź EdgesFromOctant[cube.octant];
threeDirections: ThreeDirections Ź DirectionsFromOctant[cube.octant];
[] Ź Inner[grandNeighbors[DirectionFromOctant[cube.octant]]];
[] Ź Inner[grandNeighbors[DirectionFromEdge[threeEdges.e0]]];
[] Ź Inner[grandNeighbors[DirectionFromEdge[threeEdges.e1]]];
[] Ź Inner[grandNeighbors[DirectionFromEdge[threeEdges.e2]]];
[] Ź Inner[grandNeighbors[threeDirections.d0]];
[] Ź Inner[grandNeighbors[threeDirections.d1]];
[] Ź Inner[grandNeighbors[threeDirections.d2]];
};
};
};
ApplyToTerminal[root, Inner];
};
<<Searching>>
PointInCube: PUBLIC PROC [point: Triple, cube: Cube, fudge: REAL Ź 0.] RETURNS [BOOL] ~ {
IF cube = NIL
THEN RETURN[FALSE]
ELSE {
extra: REAL ~ fudge*Size[cube];
RETURN[
point.x >= cube.corners[lbn].point.x-extra AND
point.x <= cube.corners[rbn].point.x+extra AND
point.y >= cube.corners[lbn].point.y-extra AND
point.y <= cube.corners[ltn].point.y+extra AND
point.z >= cube.corners[lbn].point.z-extra AND
point.z <= cube.corners[lbf].point.z+extra];
};
};

WhichCube: PUBLIC PROC [point: Triple, root: Cube] RETURNS [Cube] ~ {
Test: PROC [cube: Cube] RETURNS [Cube] ~ {
FOR octant: Octant IN Octant DO
kid: Cube Ź cube.kids[octant];
IF kid # NIL AND PointInCube[point, kid] THEN RETURN[Test[kid]];
ENDLOOP;
RETURN[cube];
};
RETURN[IF PointInCube[point, root] THEN Test[root] ELSE NIL];
};

FirstKid: PUBLIC PROC [cube: Cube] RETURNS [Cube] ~ {
IF cube = NIL THEN RETURN[NIL];
FOR o: Octant IN Octant DO IF cube.kids[o] # NIL THEN RETURN[cube.kids[o]]; ENDLOOP;
RETURN[NIL];
};

FirstTerminalKid: PUBLIC PROC [cube: Cube] RETURNS [Cube] ~ {
kid: Cube Ź cube;
WHILE (kid # NIL AND NOT kid.terminal) DO kid Ź FirstKid[kid]; ENDLOOP;
RETURN[kid];
};
<<Callback>>
Apply: PUBLIC PROC [cube: Cube, cubeProc: CubeProc] ~ {
Inner: PROC [cube: Cube] ~ {
IF cube # NIL AND cubeProc[cube]
THEN FOR o: Octant IN Octant DO Inner[cube.kids[o]]; ENDLOOP;
};
Inner[cube];
};

ApplyToKids: PUBLIC PROC [cube: Cube, cubeProc: CubeProc] ~ {
Inner: PROC [cube: Cube] ~ {
IF cube # NIL AND cubeProc[cube]
THEN FOR o: Octant IN Octant DO Inner[cube.kids[o]]; ENDLOOP;
};
IF cube # NIL THEN FOR o: Octant IN Octant DO Inner[cube.kids[o]]; ENDLOOP;
};

ApplyToTerminal: PUBLIC PROC [cube: Cube, cubeProc: CubeProc] ~ {
Inner: PROC [cube: Cube] ~ {
IF cube = NIL THEN RETURN;
IF cube.terminal
THEN {IF NOT cubeProc[cube] THEN RETURN}
ELSE FOR o: Octant IN Octant DO Inner[cube.kids[o]]; ENDLOOP;
};
Inner[cube];
};

ApplyToLevel: PUBLIC PROC [cube: Cube, cubeProc: CubeProc, level: INT] ~ {
Inner: PROC [cube: Cube] ~ {
IF cube = NIL OR cube.level > level THEN RETURN;
IF cube.level = level
THEN {IF NOT cubeProc[cube] THEN RETURN}
ELSE FOR o: Octant IN Octant DO Inner[cube.kids[o]]; ENDLOOP;
};
Inner[cube];
};

ApplyToTerminalCorners: PUBLIC PROC [cube: Cube, cornerProc: CornerProc] ~ {
abort: ERROR = CODE;
cubeProc: CubeProc ~ {
Test: PROC [d: Direction, o: Octant] ~ {
IF neighbors[d] = NIL THEN IF NOT cornerProc[cube.corners[o]] THEN ERROR abort;
};
neighbors: Neighborhood ~ Neighbors[cube];
IF NOT cornerProc[cube.corners[lbn]] THEN ERROR abort;
Test[f, lbf]; Test[rf, rbf]; Test[r, rbn]; Test[t, ltn]; Test[tf, ltf]; Test[rtf, rtf]; Test[rt, rtn];
};
ApplyToTerminal[cube, cubeProc ! abort => CONTINUE];
};

ApplyToNonParentKidCorners: PUBLIC PROC [cube: Cube, cornerProc: CornerProc] ~ {
Apply: PROC [d: Direction] RETURNS [BOOL] ~ {RETURN[cornerProc[directionCorners[d]]]};
directionCorners: DirectionCorners ~ GetDirectionCorners[cube];
IF NOT Apply[ln] THEN RETURN;            -- midpoint of ln edge
IF NOT Apply[lf] THEN RETURN;            -- midpoint of lf edge
IF NOT Apply[rn] THEN RETURN;        -- midpoint of rn edge
IF NOT Apply[rf] THEN RETURN;            -- midpoint of rf edge
IF NOT Apply[lb] THEN RETURN;            -- midpoint of lb edge
IF NOT Apply[lt] THEN RETURN;            -- midpoint of lt edge
IF NOT Apply[rb] THEN RETURN;            -- midpoint of rb edge
IF NOT Apply[rt] THEN RETURN;            -- midpoint of rt edge
IF NOT Apply[bn] THEN RETURN;        -- midpoint of bn edge
IF NOT Apply[tn] THEN RETURN;            -- midpoint of tn edge
IF NOT Apply[bf] THEN RETURN;            -- midpoint of bf edge
IF NOT Apply[tf] THEN RETURN;            -- midpoint of tf edge
IF NOT Apply[l] THEN RETURN;            -- center of l face
IF NOT Apply[r] THEN RETURN;            -- center of r face
IF NOT Apply[b] THEN RETURN;            -- center of b face
IF NOT Apply[t] THEN RETURN;            -- center of t face
IF NOT Apply[n] THEN RETURN;            -- center of n face
IF NOT Apply[f] THEN RETURN;            -- center of f face
IF NOT Apply[c] THEN RETURN;            -- center
};
<<Directions>>
Dir: TYPE ~ RECORD [x --(l-r)--, y --(b-t)--, z --(n-f)--: [-1..1], none: BOOL Ź FALSE];

opposites:        ARRAY Direction OF Direction ~ InitOpposites[];
sums:            ARRAY Direction OF ARRAY Direction OF Direction ~ InitSums[];
fromOctants:    ARRAY Octant OF ARRAY Octant OF Direction ~ InitFromOctants[];

decode: ARRAY Direction OF Dir ~ [
none:    [0,        0,        0,    TRUE],
c:        [0,        0,        0],
l:        [-1,    0,        0],
r:        [1,        0,        0],
b:        [0,        -1,    0],
t:        [0,        1,        0],
n:        [0,        0,        -1],
f:        [0,        0,        1],
lb:        [-1,    -1,    0],
lt:        [-1,    1,        0],
ln:        [-1,    0,        -1],
lf:        [-1,    0,        1],
rb:        [1,        -1,    0],
rt:        [1,        1,        0],
rn:    [1,        0,        -1],
rf:        [1,        0,        1],
bn:    [0,        -1,    -1],
bf:        [0,        -1,    1],
tn:        [0,        1,        -1],
tf:        [0,        1,        1],
lbn:    [-1,    -1,    -1],
lbf:    [-1,    -1,    1],
ltn:    [-1,    1,        -1],
ltf:    [-1,    1,        1],
rbn:    [1,        -1,    -1],
rbf:    [1,        -1,    1],
rtn:    [1,        1,        -1],
rtf:    [1,        1,        1]
];

Encode: PROC [x, y, z: INT] RETURNS [Direction] ~ {
FOR d: Direction IN Direction DO
dir: Dir ~ decode[d];
IF NOT dir.none AND dir.x = x AND dir.y = y AND dir.z = z THEN RETURN [d]
ENDLOOP;
RETURN [none]
};

InitOpposites: PROC RETURNS [a: ARRAY Direction OF Direction] ~ INLINE {
Opposite: PROC [d: Direction] RETURNS [Direction] ~ {
a: Dir ~ decode[d];
RETURN [Encode[x: -a.x, y: -a.y, z: -a.z]]
};
FOR d: Direction IN Direction DO
a[d] Ź Opposite[d];
ENDLOOP;
};

InitSums: PROC RETURNS [a: ARRAY Direction OF ARRAY Direction OF Direction] ~ INLINE {
Sum: PROC [d0, d1: Direction] RETURNS [Direction] ~ {
d0Dir: Dir ~ decode[d0];
d1Dir: Dir ~ decode[d1];
RETURN[SELECT TRUE FROM
d0Dir.none => d1,
d1Dir.none => d0,
ENDCASE => Encode[d0Dir.x+d1Dir.x, d0Dir.y+d1Dir.y, d0Dir.z+d1Dir.z]];
};
FOR d0: Direction IN Direction DO
FOR d1: Direction IN Direction DO
a[d0][d1] Ź Sum[d0, d1];
ENDLOOP;
ENDLOOP;
};

InitFromOctants: PROC RETURNS [a: ARRAY Octant OF ARRAY Octant OF Direction]~INLINE {
FromOctants: PROC [octant0, octant1: Octant] RETURNS [Direction] ~ {
Dif: PROC [b0, b1: [-1..1]] RETURNS [[-1..1]] ~ {
RETURN[SELECT b1 FROM b0 => 0, > b0 => 1, ENDCASE => -1];
};
d0: Dir ~ decode[DirectionFromOctant[octant0]];
d1: Dir ~ decode[DirectionFromOctant[octant1]];
RETURN[Encode[Dif[d0.x, d1.x], Dif[d0.y, d1.y], Dif[d0.z, d1.z]]];
};
FOR o0: Octant IN Octant DO
FOR o1: Octant IN Octant DO a[o0][o1] Ź FromOctants[o0, o1]; ENDLOOP;
ENDLOOP;
};

DirectionTypeFromDirection: PUBLIC PROC [direction: Direction]
RETURNS [DirectionType]
~ {
RETURN[SELECT direction FROM
l, r, b, t, n, f => face,
lb, lt, ln, lf, rb, rt, rn, rf, bn, bf, tn, tf => edge,
lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf => corner,
ENDCASE => none];
};

DirectionFromFace: PUBLIC PROC [face: Face] RETURNS [Direction] ~ {
RETURN[SELECT face FROM l => l, r => r, b => b, t => t, n => n, f => f, ENDCASE => none];
};

DirectionFromEdge: PUBLIC PROC [edge: Edge] RETURNS [Direction] ~ {
RETURN[SELECT edge FROM
lb=>lb, lt=>lt, rb=>rb, rt=>rt, ln=>ln, lf=>lf, rn=>rn, rf=>rf, bn=>bn, bf=>bf, tn=>tn, tf=>tf,
ENDCASE => none];
};

DirectionFromOctant: PUBLIC PROC [octant: Octant] RETURNS [Direction] ~ {
RETURN[SELECT octant FROM
lbn => lbn, lbf => lbf, ltn => ltn, ltf => ltf, rbn => rbn, rbf => rbf, rtn => rtn, rtf => rtf,
ENDCASE => none];
};

DirectionFromOctants: PUBLIC PROC [octant0, octant1: Octant] RETURNS [Direction] ~ {
RETURN[fromOctants[octant0][octant1]];
};

DirectionFromPoints: PUBLIC PROC [point0, point1: Triple] RETURNS [Direction] ~ {
xSame: BOOL ~ point0.x = point1.x;
ySame: BOOL ~ point0.y = point1.y;
zSame: BOOL ~ point0.z = point1.z;
IF xSame THEN {
IF ySame THEN RETURN[IF point1.z > point0.z THEN f ELSE n];
IF zSame THEN RETURN[IF point1.y > point0.y THEN t ELSE b];
}
ELSE IF ySame THEN {
IF xSame THEN RETURN[IF point1.z > point0.z THEN f ELSE n];
IF zSame THEN RETURN[IF point1.x > point0.x THEN r ELSE l];
}
ELSE IF zSame THEN {
IF xSame THEN RETURN[IF point1.y > point0.y THEN t ELSE b];
IF ySame THEN RETURN[IF point1.x > point0.x THEN r ELSE l];
};
RETURN[none];
};

OppositeFace: PUBLIC PROC [face: Face] RETURNS [Face] ~ {
RETURN[SELECT face FROM l => r, r => l, b => t, t => b, n => f, ENDCASE => n];
};

OppositeDirection: PUBLIC PROC [direction: Direction] RETURNS [Direction] ~ {
RETURN[opposites[direction]];
};

AddDirection: PUBLIC PROC [d0, d1: Direction] RETURNS [Direction] ~ {
RETURN[sums[d0][d1]];
};

EdgeFromDirection: PUBLIC PROC [d: Direction] RETURNS [Edge] ~ {
RETURN[SELECT d FROM
ln    => ln, lf    => lf, rn    => rn, rf    => rf, bn    => bn, bf    => bf, tn    => tn, tf    => tf,
ENDCASE => ERROR];
};
<<>>
OctantFromDirection: PUBLIC PROC [d: Direction] RETURNS [Octant] ~ {
RETURN[SELECT d FROM
lbn => lbn, lbf => lbf, ltn => ltn, ltf => ltf, rbn => rbn, rbf => rbf, rtn => rtn, rtf => rtf,
ENDCASE => ERROR];
};
<<>>
OctantFromThreeDirections: PUBLIC PROC [d0, d1, d2: Direction] RETURNS [o: Octant] ~ {
o Ź OctantFromDirection[AddDirection[AddDirection[d0, d1], d2]];
};

FaceFromDirection: PUBLIC PROC [d: Direction] RETURNS [Face] ~ {
RETURN[SELECT d FROM
l => l, r => r, b => b, t => t, n => n, f => f,
ENDCASE => ERROR];
};

FaceDirections: PUBLIC PROC [face: Face] RETURNS [fd: FourDirections] ~ {
fd Ź SELECT face FROM
l => [lb, lt, ln, lf], r => [rb, rt, rn, rf], b => [lb, rb, bn, bf],
t => [lt, rt, tn, tf], n => [ln, rn, bn, tn], ENDCASE => [lf, rf, bf, tf];
};

NormalFromFace: PUBLIC PROC [face: Face] RETURNS [Triple] ~ {
RETURN[SELECT face FROM
l => [-1.0, 0.0, 0.0], r => [1.0, 0.0, 0.0], b => [0.0, -1.0, 0.0],
t => [0.0, 1.0, 0.0], n => [0.0, 0.0, -1.0], ENDCASE => [0.0, 0.0, 1.0]];
};
<<Edge-Face-Octant>>
PlaneFromCubeFace: PUBLIC PROC [cube: Cube, face: Face] RETURNS [p: Plane] ~ {
o: Octant Ź SELECT face FROM l=>lbn, r=>rbn, b=>lbn, t=>ltn, n=>lbn, ENDCASE=>lbf;
p Ź G3dPlane.FromPointAndNormal[cube.corners[o].point, NormalFromFace[face]];
};

EdgesFromOctant: PUBLIC PROC [octant: Octant] RETURNS [ThreeEdges] ~ {
RETURN[SELECT octant FROM
lbn => [lb, ln, bn], lbf => [lb, lf, bf], ltn => [lt, ln, tn], ltf => [lt, lf, tf],
rbn => [rb, rn, bn], rbf => [rb, rf, bf], rtn => [rt, rn, tn], ENDCASE => [rt, rf, tf]];
};
<<>>
DirectionsFromOctant: PUBLIC PROC [octant: Octant] RETURNS [ThreeDirections] ~ {
RETURN[SELECT octant FROM
lbn => [l, b, n], lbf => [l, b, f], ltn => [l, t, n], ltf => [l, t, f],
rbn => [r, b, n], rbf => [r, b, f], rtn => [r, t, n], ENDCASE => [r, t, f]];
};

FaceFromEdgeOctant: PUBLIC PROC [edge: Edge, octant: Octant] RETURNS [Face] ~ {
RETURN[SELECT edge FROM
lb =>            IF octant = lbn    THEN l    ELSE b,
lt =>            IF octant = ltn    THEN t    ELSE l,
ln =>            IF octant = lbn    THEN n    ELSE l,
lf =>            IF octant = lbf    THEN l    ELSE f,
rb =>            IF octant = rbn    THEN b    ELSE r,
rt =>            IF octant = rtn    THEN r    ELSE t,
rn =>            IF octant = rbn    THEN r    ELSE n,
rf =>            IF octant = rbf    THEN f    ELSE r,
bn =>            IF octant = lbn    THEN b    ELSE n,
bf =>            IF octant = lbf    THEN f    ELSE b,
tn =>            IF octant = ltn    THEN n    ELSE t,
ENDCASE =>    IF octant = ltf    THEN t    ELSE f];
};

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

NextCCWEdge: PUBLIC PROC [edge: Edge, face: Face] RETURNS [Edge] ~ {
RETURN[SELECT edge FROM
lb =>            IF face = l    THEN ln    ELSE bf,
lt =>            IF face = l    THEN lf    ELSE tn,
ln =>            IF face = l    THEN lt    ELSE bn,
lf =>            IF face = l    THEN lb    ELSE tf,
rb =>            IF face = r    THEN rf    ELSE bn,
rt =>            IF face = r    THEN rn    ELSE tf,
rn =>            IF face = r    THEN rb    ELSE tn,
rf =>            IF face = r    THEN rt    ELSE bf,
bn =>            IF face = b    THEN lb    ELSE rn,
bf =>            IF face = b    THEN rb    ELSE lf,
tn =>            IF face = t    THEN rt    ELSE ln,
ENDCASE =>    IF face = t    THEN lt    ELSE rf];
};

EdgeDirections: PUBLIC PROC [edge: Edge] RETURNS [ThreeDirections] ~ {
RETURN[SELECT edge FROM
lb => [l, b, lb], lt => [l, t, lt], ln => [l, n, ln], lf => [l, f, lf], rb => [r, b, rb],
rt => [r, t, rt], rn => [r, n, rn], rf => [r, f, rf], bn => [b, n, bn],
bf => [b, f, bf], tn => [t, n, tn], tf => [t, f, tf], ENDCASE => [none, none, none]];
};

EdgeFaces: PUBLIC PROC [edge: Edge] RETURNS [TwoFaces] ~ {
RETURN[SELECT edge FROM
lb => [l, b], lt => [l, t], ln => [l, n], lf => [l, f], rb => [r, b], rt => [r, t],
rn => [r, n], rf => [r, f], bn => [b, n], bf => [b, f], tn => [t, n], ENDCASE => [t, f]];
};

EdgeOctants: PUBLIC PROC [edge: Edge] RETURNS [TwoOctants] ~ {
RETURN[SELECT edge FROM
lb => [lbn, lbf], lt => [ltn, ltf], ln => [lbn, ltn], lf => [lbf, ltf],
rb => [rbn, rbf], rt => [rtn, rtf], rn => [rbn, rtn], rf => [rbf, rtf],
bn => [lbn, rbn], bf => [lbf, rbf], tn => [ltn, rtn], ENDCASE    =>    [ltf, rtf]];
};

FaceOctants: PUBLIC PROC [face: Face] RETURNS [FourOctants] ~ {
RETURN[SELECT face FROM
l            => [lbf, ltf, ltn, lbn],
r            => [rbn, rtn, rtf, rbf],
b            => [lbf, lbn, rbn, rbf],
t            => [ltn, ltf, rtf, rtn],
n            => [lbn, ltn, rtn, rbn],
ENDCASE    => [rbf, rtf, ltf, lbf]];
};

EdgeCorners: PUBLIC PROC [cube: Cube, edge: Edge] RETURNS [TwoCorners] ~ {
octants: TwoOctants ~ EdgeOctants[edge];
RETURN[[cube.corners[octants.o0], cube.corners[octants.o1]]];
};

FaceCorners: PUBLIC PROC [cube: Cube, face: Face] RETURNS [FourCorners] ~ {
o: FourOctants ~ FaceOctants[face];
RETURN[[cube.corners[o.o0], cube.corners[o.o1], cube.corners[o.o2], cube.corners[o.o3]]];
};

OppositeEdge: PUBLIC PROC [edge: Edge, face: Face] RETURNS [Edge] ~ {
RETURN[SELECT edge FROM
lb =>            IF face = l    THEN lt    ELSE rb,
lt =>            IF face = l    THEN lb    ELSE rt,
ln =>            IF face = l    THEN lf    ELSE rn,
lf =>            IF face = l    THEN ln    ELSE rf,
rb =>            IF face = r    THEN rt    ELSE lb,
rt =>            IF face = r    THEN rb    ELSE lt,
rn =>            IF face = r    THEN rf    ELSE ln,
rf =>            IF face = r    THEN rn    ELSE lf,
bn =>            IF face = b    THEN bf    ELSE tn,
bf =>            IF face = b    THEN bn    ELSE tf,
tn =>            IF face = t    THEN tf    ELSE bn,
ENDCASE =>    IF face = t    THEN tn    ELSE bf];
};

InverseEdge: PUBLIC PROC [edge: Edge, face: Face] RETURNS [Edge] ~ {
RETURN[SELECT edge FROM
lb =>            IF face = l    THEN rb    ELSE lt,
lt =>            IF face = l    THEN rt    ELSE lb,
ln =>            IF face = l    THEN rn    ELSE lf,
lf =>            IF face = l    THEN rf    ELSE ln,
rb =>            IF face = r    THEN lb    ELSE rt,
rt =>            IF face = r    THEN lt    ELSE rb,
rn =>            IF face = r    THEN ln    ELSE rf,
rf =>            IF face = r    THEN lf    ELSE rn,
bn =>            IF face = b    THEN tn    ELSE bf,
bf =>            IF face = b    THEN tf    ELSE bn,
tn =>            IF face = t    THEN bn    ELSE tf,
ENDCASE =>    IF face = t    THEN bf    ELSE tn];
};

DiagonalEdge: PUBLIC PROC [edge: Edge] RETURNS [Edge] ~ {
RETURN[SELECT edge FROM
lb => rt, lt => rb, ln => rf, lf => rn, rb => lt, rt => lb, rn => lf,
rf => ln, bn => tf, bf => tn, tn => bf, ENDCASE => bn];
};

OtherFace: PUBLIC PROC [edge: Edge, face: Face] RETURNS [Face] ~ {
twoFaces: TwoFaces Ź EdgeFaces[edge];
RETURN[IF face = twoFaces.f0 THEN twoFaces.f1 ELSE twoFaces.f1];
};

AxisFromEdge: PUBLIC PROC [edge: Edge] RETURNS [Axis] ~ {
RETURN[SELECT edge FROM bn, tn, bf, tf => x, ln, lf, rn, rf => y, ENDCASE => z];
};
<<Neighbors>>
GetFaceNeighborInfo: PROC [cube: Cube, face: Face] RETURNS [FaceNeighborInfo] ~ {
octantWord: WORD ~ ORD[cube.octant];
faceBit: WORD ~ SELECT face FROM l, r => 4, b, t => 2, ENDCASE => 1;
faceWord: WORD ~ Basics.BITAND[octantWord, faceBit];
RETURN[[
VAL[Basics.BITXOR[octantWord, faceBit]],
SELECT face FROM l, b, n => faceWord = 0, ENDCASE => faceWord # 0
]];
};

FaceNeighbor: PUBLIC PROC [cube: Cube, face: Face] RETURNS [Cube] ~ {
IF cube = NIL OR cube.parent = NIL
THEN RETURN[NIL]
ELSE {
faceNeighborInfo: FaceNeighborInfo ~ GetFaceNeighborInfo[cube, face];
IF faceNeighborInfo.recurse
THEN {
parentNeighbor: Cube Ź FaceNeighbor[cube.parent, face];
IF parentNeighbor = NIL
THEN RETURN[NIL]
ELSE RETURN[parentNeighbor.kids[faceNeighborInfo.nOctant]];
}
ELSE RETURN[cube.parent.kids[faceNeighborInfo.nOctant]];
};
};

EdgeNeighbors: PUBLIC PROC [cube: Cube, edge: Edge] RETURNS [cubes: ThreeCubes] ~ {
f1, f2: Face;
SELECT edge FROM
lb    => {f1 Ź l;    f2 Ź b};
lt    => {f1 Ź l;    f2 Ź t};
rb    => {f1 Ź r;    f2 Ź b};
rt    => {f1 Ź r;    f2 Ź t};
ln    => {f1 Ź l;    f2 Ź n};
lf    => {f1 Ź l;    f2 Ź f};
rn    => {f1 Ź r;    f2 Ź n};
rf    => {f1 Ź r;    f2 Ź f};
bn    => {f1 Ź b;    f2 Ź n};
bf    => {f1 Ź b;    f2 Ź f};
tn    => {f1 Ź t;    f2 Ź n};
tf    => {f1 Ź t;    f2 Ź f};
ENDCASE;
cubes.c0 Ź FaceNeighbor[cube, f1];
cubes.c1 Ź FaceNeighbor[cube, f2];
cubes.c2 Ź FaceNeighbor[cubes.c0, f2];
};

Neighbor: PUBLIC PROC [cube: Cube, direction: Direction] RETURNS [x: Cube] ~ {
IF cube = NIL THEN RETURN[NIL];
SELECT direction FROM
<<six faces:>>
l => x Ź FaceNeighbor[cube, l];
r => x Ź FaceNeighbor[cube, r];
b => x Ź FaceNeighbor[cube, b];
t => x Ź FaceNeighbor[cube, t];
n => x Ź FaceNeighbor[cube, n];
f => x Ź FaceNeighbor[cube, f];
<<twelve edges:>>
lb => {
l: Cube Ź FaceNeighbor[cube, l];
x Ź IF l # NIL THEN FaceNeighbor[l, b] ELSE FaceNeighbor[FaceNeighbor[cube, b], l];
};
lt => {
l: Cube Ź FaceNeighbor[cube, l];
x Ź IF l # NIL THEN FaceNeighbor[l, t] ELSE FaceNeighbor[FaceNeighbor[cube, t], l];
};
ln => {
l: Cube Ź FaceNeighbor[cube, l];
x Ź IF l # NIL THEN FaceNeighbor[l, n] ELSE FaceNeighbor[FaceNeighbor[cube, n], l];
};
lf => {
l: Cube Ź FaceNeighbor[cube, l];
x Ź IF l # NIL THEN FaceNeighbor[l, f] ELSE FaceNeighbor[FaceNeighbor[cube, f], l];
};
rb => {
r: Cube Ź FaceNeighbor[cube, r];
x Ź IF r # NIL THEN FaceNeighbor[r, b] ELSE FaceNeighbor[FaceNeighbor[cube, b], r];
};
rt => {
r: Cube Ź FaceNeighbor[cube, r];
x Ź IF r # NIL THEN FaceNeighbor[r, t] ELSE FaceNeighbor[FaceNeighbor[cube, t], r];
};
rn => {
r: Cube Ź FaceNeighbor[cube, r];
x Ź IF r # NIL THEN FaceNeighbor[r, n] ELSE FaceNeighbor[FaceNeighbor[cube, n], r];
};
rf => {
r: Cube Ź FaceNeighbor[cube, r];
x Ź IF r # NIL THEN FaceNeighbor[r, f] ELSE FaceNeighbor[FaceNeighbor[cube, f], r];
};
bn => {
b: Cube Ź FaceNeighbor[cube, b];
x Ź IF b # NIL THEN FaceNeighbor[b, n] ELSE FaceNeighbor[FaceNeighbor[cube, n],b];
};
bf => {
b: Cube Ź FaceNeighbor[cube, b];
x Ź IF b # NIL THEN FaceNeighbor[b, f] ELSE FaceNeighbor[FaceNeighbor[cube, f],b ];
};
tn => {
t: Cube Ź FaceNeighbor[cube, t];
x Ź IF t # NIL THEN FaceNeighbor[t, n] ELSE FaceNeighbor[FaceNeighbor[cube, n], t];
};
tf => {
t: Cube Ź FaceNeighbor[cube, t];
x Ź IF t # NIL THEN FaceNeighbor[t, f] ELSE FaceNeighbor[FaceNeighbor[cube, f], t];
};
<<eight corners:>>
lbn => {
temp: Cube Ź Neighbor[cube, lb];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, n]];
temp Ź Neighbor[cube, ln];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, b]];
RETURN[FaceNeighbor[Neighbor[cube, bn], l]];
};
lbf => {
temp: Cube Ź Neighbor[cube, lb];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, f]];
temp Ź Neighbor[cube, lf];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, b]];
RETURN[FaceNeighbor[Neighbor[cube, bf], l]];
};
ltn => {
temp: Cube Ź Neighbor[cube, lt];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, n]];
temp Ź Neighbor[cube, ln];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, t]];
RETURN[FaceNeighbor[Neighbor[cube, tn], l]];
};
ltf => {
temp: Cube Ź Neighbor[cube, lt];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, f]];
temp Ź Neighbor[cube, lf];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, t]];
RETURN[FaceNeighbor[Neighbor[cube, tf], l]];
};
rbn => {
temp: Cube Ź Neighbor[cube, rb];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, n]];
temp Ź Neighbor[cube, rn];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, b]];
RETURN[FaceNeighbor[Neighbor[cube, bn], r]];
};
rbf => {
temp: Cube Ź Neighbor[cube, rb];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, f]];
temp Ź Neighbor[cube, rf];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, b]];
RETURN[FaceNeighbor[Neighbor[cube, bf], r]];
};
rtn => {
temp: Cube Ź Neighbor[cube, rt];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, n]];
temp Ź Neighbor[cube, rn];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, t]];
RETURN[FaceNeighbor[Neighbor[cube, tn], r]];
};
rtf => {
temp: Cube Ź Neighbor[cube, rt];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, f]];
temp Ź Neighbor[cube, rf];
IF temp # NIL THEN RETURN[FaceNeighbor[temp, t]];
RETURN[FaceNeighbor[Neighbor[cube, tf], r]];
};
ENDCASE    => x Ź cube;
};

Neighbors: PUBLIC PROC [cube: Cube] RETURNS [neighbors: Neighborhood] ~ {
IF cube # NIL THEN
FOR direction: Direction IN Direction DO
neighbors[direction] Ź Neighbor[cube, direction];
ENDLOOP;
};

Move: PUBLIC PROC [cube: Cube, direction: Direction] RETURNS [Neighborhood] ~ {
RETURN[Neighbors[Neighbor[cube, direction]]];
};
<<Children>>
GetDirectionCorner: PUBLIC PROC [cube: Cube, d: Direction] RETURNS [c: Corner] ~ {
c Ź SELECT d FROM
lbn    => cube.corners[lbn],                    -- lbn corner
lbf    => cube.corners[lbf],                    -- lbf corner
ltn    => cube.corners[ltn],                    -- ltn corner
ltf        => cube.corners[ltf],                        -- ltf corner
rbn    => cube.corners[rbn],                    -- rbn corner
rbf    => cube.corners[rbf],                    -- rbf corner
rtn    => cube.corners[rtn],                    -- rtn corner
rtf        => cube.corners[rtf],                    -- rtf corner
ln        => cube.kids[ltn].corners[lbn],            -- midpoint of ln edge
lf        => cube.kids[ltf].corners[lbf],            -- midpoint of lf edge
rn        => cube.kids[rtn].corners[rbn],            -- midpoint of rn edge
rf        => cube.kids[rtf].corners[rbf],            -- midpoint of rf edge
lt        => cube.kids[ltn].corners[ltf],            -- midpoint of lt edge
lb        => cube.kids[lbn].corners[lbf],            -- midpoint of lb edge
rt        => cube.kids[rtn].corners[rtf],            -- midpoint of rt edge
rb        => cube.kids[rbn].corners[rbf],            -- midpoint of rb edge
tn        => cube.kids[ltn].corners[rtn],            -- midpoint of tn edge
bn        => cube.kids[lbn].corners[rbn],        -- midpoint of bn edge
tf        => cube.kids[ltf].corners[rtf],            -- midpoint of tf edge
bf        => cube.kids[lbf].corners[rbf],            -- midpoint of bf edge
l        => cube.kids[lbn].corners[ltf],            -- center of l face
r        => cube.kids[rbn].corners[rtf],            -- center of r face
b        => cube.kids[lbn].corners[rbf],            -- center of b face
t        => cube.kids[ltn].corners[rtf],            -- center of t face
n        => cube.kids[lbn].corners[rtn],            -- center of n face
f        => cube.kids[lbf].corners[rtf],            -- center of f face
c        => cube.kids[lbn].corners[rtf],            -- center of cube
ENDCASE => NIL;
};

GetDirectionCorners: PUBLIC PROC [cube: Cube]
RETURNS [directionCorners: DirectionCorners]
~ {
directionCorners[lbn]    Ź cube.corners[lbn];                    -- lbn corner
directionCorners[lbf]    Ź cube.corners[lbf];                        -- lbf corner
directionCorners[ltn]    Ź cube.corners[ltn];                        -- ltn corner
directionCorners[ltf]    Ź cube.corners[ltf];                        -- ltf corner
directionCorners[rbn]    Ź cube.corners[rbn];                    -- rbn corner
directionCorners[rbf]    Ź cube.corners[rbf];                        -- rbf corner
directionCorners[rtn]    Ź cube.corners[rtn];                        -- rtn corner
directionCorners[rtf]    Ź cube.corners[rtf];                        -- rtf corner
directionCorners[ln]        Ź cube.kids[ltn].corners[lbn];            -- midpoint of ln edge
directionCorners[lf]        Ź cube.kids[ltf].corners[lbf];            -- midpoint of lf edge
directionCorners[rn]    Ź cube.kids[rtn].corners[rbn];            -- midpoint of rn edge
directionCorners[rf]        Ź cube.kids[rtf].corners[rbf];            -- midpoint of rf edge
directionCorners[lt]        Ź cube.kids[ltn].corners[ltf];            -- midpoint of lt edge
directionCorners[lb]        Ź cube.kids[lbn].corners[lbf];            -- midpoint of lb edge
directionCorners[rt]        Ź cube.kids[rtn].corners[rtf];            -- midpoint of rt edge
directionCorners[rb]        Ź cube.kids[rbn].corners[rbf];            -- midpoint of rb edge
directionCorners[tn]        Ź cube.kids[ltn].corners[rtn];            -- midpoint of tn edge
directionCorners[bn]    Ź cube.kids[lbn].corners[rbn];            -- midpoint of bn edge
directionCorners[tf]        Ź cube.kids[ltf].corners[rtf];            -- midpoint of tf edge
directionCorners[bf]        Ź cube.kids[lbf].corners[rbf];            -- midpoint of bf edge
directionCorners[l]        Ź cube.kids[lbn].corners[ltf];            -- center of l face
directionCorners[r]        Ź cube.kids[rbn].corners[rtf];            -- center of r face
directionCorners[b]        Ź cube.kids[lbn].corners[rbf];            -- center of b face
directionCorners[t]        Ź cube.kids[ltn].corners[rtf];            -- center of t face
directionCorners[n]        Ź cube.kids[lbn].corners[rtn];            -- center of n face
directionCorners[f]        Ź cube.kids[lbf].corners[rtf];            -- center of f face
directionCorners[c]        Ź cube.kids[lbn].corners[rtf];            -- center of cube
};
<<Miscellany>>
SetTerminalField: PUBLIC PROC [root: Cube] ~ {
cubeProc: CubeProc ~ {
FOR o: Octant IN Octant DO
IF cube.kids[o] # NIL THEN {
cube.terminal Ź FALSE;
RETURN;
};
ENDLOOP;
cube.terminal Ź TRUE;
};
Apply[root, cubeProc];
};

SetLevels: PUBLIC PROC [root: Cube] ~ {
SetKidsLevel: PROC [cube: Cube, level: INT] ~ {
nextLevel: INT ~ level+1;
FOR o: Octant IN Octant DO
kid: Cube Ź cube.kids[o];
IF kid = NIL THEN LOOP;
kid.level Ź level;
SetKidsLevel[kid, nextLevel];
ENDLOOP;
};
root.level Ź 0;
SetKidsLevel[root, 1];
};
    GetClientData: PUBLIC PROC [cube: Cube] RETURNS [REF ANY] ~ {
        RETURN [cube.refAny];
        };
    SetClientData: PUBLIC PROC [cube: Cube, r: REF ANY]  ~ {
        IF cube # NIL THEN cube.refAny Ź r;
        };

<<Object Query>>
<<By Ken Fishkin.>>
<<These operations are performed using the refAny field of the cube to store a pointer>>
<<to a record containing a list of all objects which (at least partially) are within this cube.>>
<<This is NIL for non-terminal cubes.>>
ObjectTrack:        TYPE ~ REF ObjectTrackRep;
ObjectTrackRep:    TYPE ~ RECORD [objects: LIST OF CubeObject  _ NIL];

    <<PrintBox: PROC [b: Box, s: IO.STREAM, hdr: ROPE Ź NIL]  ~ {>>
        <<IF hdr # NIL THEN IO.PutRope[s,hdr];>>
        <<G3dIO.PrintTriple[s,b.min,"min:"];>>
        <<G3dIO.PrintTriple[s,b.max,"max:"];>>
        <<};>>

    FromObjects: PUBLIC PROC [objs: LIST OF CubeObject] RETURNS [Cube] ~ {
        bbox: Box Ź objs.first.bbox;
        size: REAL;
        ans: Cube Ź NIL;
        track: ObjectTrack Ź NEW[ObjectTrackRep];

        <<first, go through and find the bounding box of all the objects: needed to create the Cube.>>
        FOR l: LIST OF CubeObject _ objs, l.rest WHILE l # NIL DO
            bbox Ź G3dBox.BoxUnionBox[bbox, l.first.bbox];
            ENDLOOP;
        <<cube size is the max side length of the bbox:>>
        size Ź MAX[bbox.max.x - bbox.min.x, bbox.max.y - bbox.min.y];
        size Ź MAX[size, bbox.max.z - bbox.min.z];
        ans Ź NewCube[size, G3dBox.BoxCenter[bbox]];
        track.objects Ź objs;
        SetClientData[ans, track];
        RETURN[ans];
        };

    AddObject: PUBLIC PROC [cube: Cube, obj: CubeObject] ~ {
        track: ObjectTrack Ź NIL;
        intersection: Box;
        hasIntersection: BOOL;
        <<sanity checks....>>
        IF cube = NIL OR obj = NIL THEN RETURN;

        IF cube.refAny = NIL AND cube.terminal THEN {
            <<this can happen after subdividing. Set up an empty track, and fall into regular case.>>
            track Ź NEW[ObjectTrackRep];
            track.objects Ź NIL;
            cube.refAny Ź track;
            };
        IF NOT ISTYPE[cube.refAny,ObjectTrack] THEN RETURN;
        track Ź NARROW[cube.refAny];
        <<does this object belong within me?>>
        [hasIntersection, intersection] Ź G3dBox.Intersect[BoxOfCube[cube], obj.bbox];
        IF NOT hasIntersection THEN RETURN;
        IF cube.terminal 
            THEN {
                <<I am a terminal: add to me directly.>>
                myCo: CubeObject Ź NEW[G3dOctree.CubeObjectRep Ź
                    [val: obj.val, bbox: intersection, query: obj.query]];
                track.objects Ź CONS[myCo, track.objects];
                }
            ELSE {
                <<add to the kids. Since this guy can cross cube boundaries, might be added to more>>
                <<than one child: go through all of them.>>
                FOR octant: Octant IN Octant DO
                    AddObject[cube.kids[octant], obj];
                    ENDLOOP;
                };
        };

    SubdivideGivenObjects: PUBLIC PROC [cube: Cube, maxPerCube: INT, maxDepth: INT]  ~ {
        <<Subdivide the given cube (which must have been created with FromObjects/AddObject) until every sub-cube has either
        1) <= (maxPerCube) objects in it
        2) a depth = (maxDepth)>>
        <<>>
        track: ObjectTrack;

        IF cube = NIL OR cube.refAny = NIL OR NOT ISTYPE[cube.refAny,ObjectTrack] THEN RETURN;
        track Ź NARROW[cube.refAny];
        IF cube.level >= maxDepth THEN RETURN;
        IF cube.terminal 
            THEN {
                nObjs: INT Ź 0;
                FOR l: LIST OF CubeObject _ track.objects, l.rest WHILE l # NIL DO
                    nObjs Ź nObjs + 1;
                    ENDLOOP;
                IF nObjs <= maxPerCube THEN {
                    RETURN;
                    };
                <<can't weasel out of it anymore: have to split myself up, and hand off the objects to my kids. Since my objects may cross more 
                than one of my kids, will do this by adding each object to each kid: the kid will check and see if it lies outside 
                his boundaries.>>
                <<PrintBox[BoxOfCube[cube],TerminalIO.TOS[], IO.PutFR["%g: %g objs: Subdividing:", IO.int[cube.level], IO.int[nObjs]]];>>
                Subdivide[cube];
                FOR l: LIST OF CubeObject _ track.objects, l.rest WHILE l # NIL DO
                    FOR octant: Octant IN Octant DO
                        AddObject[cube.kids[octant], l.first];
                        ENDLOOP;
                    ENDLOOP;
                track.objects Ź NIL;
                <<and now fall down to non-terminal case (I'm not a terminal anymore):>>
                };
            <<need to make sure the kids aren't overloaded. Note that I don't have to change my bounding box, as the objects are unchanged and 
            hence have the same bounding boxes.>>
            FOR octant: Octant IN Octant DO
                SubdivideGivenObjects[cube.kids[octant], maxPerCube, maxDepth];
                ENDLOOP;

        };

    QueryObjects: PUBLIC PROC [cube: Cube, pt: Triple] RETURNS [obj: REF ANY] ~ {
        track: ObjectTrack;
        IF cube = NIL THEN RETURN [NIL];
        <<G3dIO.PrintTriple[TerminalIO.TOS[], pt, "QueryObjects:pt:"];>>
        <<PrintBox[BoxOfCube[cube], TerminalIO.TOS[], "checking against cube:"];>>
        IF cube.refAny = NIL THEN RETURN [NIL];
        IF NOT ISTYPE[cube.refAny,ObjectTrack] THEN {
            RETURN [NIL];
            };
        track Ź NARROW[cube.refAny];
        IF NOT G3dBox.PointInBox[BoxOfCube[cube], pt] THEN {
            RETURN [NIL];
            };
        IF cube.terminal 
            THEN {
                FOR l: LIST OF CubeObject _ track.objects, l.rest WHILE l # NIL DO
                    co: CubeObject Ź l.first;
                    IF G3dBox.PointInBox[co.bbox,pt] THEN {
                        IF co.query[co.val, pt] THEN RETURN [co.val];
                        };
                    ENDLOOP;
                }
            ELSE {
                <<Figure out which octant could contain this point:>>
                d0, d1, d2: Direction;
                centerPt: Triple Ź Center[cube];
                octant: Octant;
                val: REF ANY;

                d0 Ź IF pt.x > centerPt.x THEN r ELSE l;
                d1 Ź IF pt.y > centerPt.y THEN t ELSE b;
                d2 Ź IF pt.z > centerPt.z THEN f ELSE n;
                octant Ź OctantFromThreeDirections[d0,d1,d2];
                val Ź QueryObjects[cube.kids[octant], pt];
                IF val # NIL THEN RETURN [val];
                };
            RETURN [NIL];
        };
<<Miscellany>>
SetOctreeFields: PUBLIC PROC [octree: Octree] ~ {
IF octree = NIL OR octree.root = NIL THEN RETURN;
octree.depth Ź DepthOf[octree.root];
octree.radius Ź 0.5*(octree.diameter Ź Size[octree.root]);
octree.terminalRadius Ź 0.5*(octree.terminalDiameter Ź MinSize[octree.root]);
octree.nTerminalCubes Ź NTerminalCubes[octree.root];
octree.nCubes Ź NCubes[octree.root];
};
<<IO>>
WriteCubesToFile: PUBLIC PROC [fileName: ROPE, cube: Cube, miscInfo: ROPE Ź NIL] ~ {
IF cube # NIL AND fileName # NIL THEN {
Inner: PROC [cube: Cube] ~ {
AddChar: PROC [c: CHAR] ~ INLINE {
text[text.length] Ź c;
text.length Ź text.length+1;
};
text.length Ź 6;
FOR n: INT IN [0..cube.level) DO
AddChar['\t];
ENDLOOP;
FOR o: Octant IN Octant DO
AddChar[IF cube.kids[o] = NIL THEN '0 ELSE '1];
ENDLOOP;
AddChar['\n];
IO.PutText[out, text];
FOR o: Octant IN Octant DO
IF cube.kids[o] # NIL THEN Inner[cube.kids[o]];
ENDLOOP;
};
text: REF TEXT Ź RefText.ObtainScratch[100];
out: STREAM Ź FS.StreamOpen[fileName, $create];
center: Triple Ź Center[cube];
IO.PutFL[out, "Volume~ size: %g\tcenter: %g\t%g\t%g\n\n",
LIST[IO.real[cube.size], IO.real[center.x], IO.real[center.y], IO.real[center.z]]];
IF miscInfo # NIL THEN IO.PutF1[out, "MiscInfo~ %g\n\n", IO.rope[miscInfo]];
text Ź RefText.AppendRope[text, "cube:\t"];
Inner[cube];
IO.Close[out];
RefText.ReleaseScratch[text];
};
};

ReadCubesFromFile: PUBLIC PROC [fileName: ROPE] RETURNS [root: Cube] ~ {
IF fileName # NIL THEN {
InnerReadCube: PROC [parent: Cube] ~ {
CubeLineOK: PROC RETURNS [ok: BOOL Ź TRUE] ~ {
DO
text Ź IO.GetLine[in, text ! IO.EndOfStream => GOTO eof];
IF Rope.Run[id, 0, RefText.TrustTextAsRope[text], 0, TRUE] = 6 THEN RETURN;
REPEAT
eof => ok Ź FALSE;
ENDLOOP;
};
IF CubeLineOK[] THEN {
n: INT Ź 6;
WHILE text[n] = '\t DO n Ź n+1; ENDLOOP;
FOR nn: INT IN [n..n+8) DO
IF text[nn] = '1 THEN EXIT;
REPEAT
FINISHED => RETURN;
ENDLOOP;
Subdivide[parent];
FOR o: Octant IN Octant DO
IF text[n] = '0 THEN parent.kids[o] Ź NIL;
n Ź n+1;
ENDLOOP;
FOR o: Octant IN Octant DO
kid: Cube Ź parent.kids[o];
IF kid # NIL THEN InnerReadCube[kid];
ENDLOOP;
};
};
in: STREAM Ź FS.StreamOpen[FileNames.ResolveRelativePath[fileName] !
FS.Error => GOTO noOpen];
text: REF TEXT Ź RefText.ObtainScratch[100];
id: ROPE Ź "cube:\t";
line: G3dIO.Line Ź
G3dIO.FindKeyWord[in, "Octree",,, G3dIO.ObtainLine[]];
skip1: ROPE Ź G3dIO.GetWord[line];
size: REAL Ź G3dIO.GetReal[line];
skip2: ROPE Ź G3dIO.GetWord[line];
center: Triple Ź G3dIO.GetTriple[line];
root Ź NewCube[size, center];
InnerReadCube[root];
G3dIO.ReleaseLine[line];
RefText.ReleaseScratch[text];
IO.Close[in];
EXITS
noOpen => RETURN;
};
};
<<Cube Sequence Procedures>>
AddToCubeSequence: PUBLIC PROC [cube: Cube, cubes: CubeSequence Ź NIL]
RETURNS [CubeSequence]
~ {
IF cubes = NIL THEN cubes Ź NEW[CubeSequenceRep[10]];
IF cubes.length = cubes.maxLength THEN cubes Ź LengthenCubeSequence[cubes];
cubes[cubes.length] Ź cube;
cubes.length Ź cubes.length+1;
RETURN[cubes];
};

LengthenCubeSequence: PUBLIC PROC [cubes: CubeSequence, amount: REAL Ź 1.3]
RETURNS [new: CubeSequence]
~ {
new Ź NEW[CubeSequenceRep[Real.Round[1.3*cubes.length]]];
FOR i: INT IN [0..cubes.length) DO new[i] Ź cubes[i]; ENDLOOP;
new.length Ź cubes.length;
};
<<Cube Stack Support Procedures>>
StackOverflow:    PUBLIC ERROR = CODE;
StackUnderflow:    PUBLIC ERROR = CODE;

NewCubeStack: PUBLIC PROC [length: INT] RETURNS [CubeStack] ~ {
RETURN[NEW[CubeStackRep[length]]];
};

WriteBottomOfCubeStack: PUBLIC PROC [cube: Cube, cubeStack: CubeStack] ~ {
newBottom: INT ~ (cubeStack.bottom+1) MOD cubeStack.maxLength;
cubeStack[cubeStack.bottom] Ź cube;
IF newBottom = cubeStack.top THEN ERROR StackOverflow;
cubeStack.bottom Ź newBottom;
cubeStack.size Ź cubeStack.size+1;
cubeStack.maxSize Ź MAX[cubeStack.maxSize, cubeStack.size];
};

ReadTopOfCubeStack: PUBLIC PROC [cubeStack: CubeStack] RETURNS [cube: Cube] ~ {
IF cubeStack.bottom = cubeStack.top THEN ERROR StackUnderflow;
cube Ź cubeStack[cubeStack.top];
cubeStack.top Ź (cubeStack.top+1) MOD cubeStack.maxLength;
cubeStack.size Ź cubeStack.size-1;
};

CubeStackSize: PUBLIC PROC [cubeStack: CubeStack] RETURNS [INT] ~ {
RETURN[cubeStack.size];
<<RETURN[IF cubeStack.top > cubeStack.bottom>>
<<THEN 2*cubeStack.maxLength-cubeStack.top-cubeStack.bottom>>
<<ELSE cubeStack.bottom-cubeStack.top];>>
};

MaxCubeStackSize: PUBLIC PROC [cubeStack: CubeStack] RETURNS [INT] ~ {
RETURN[cubeStack.maxSize];
};

CubeStackEmpty: PUBLIC PROC [cubeStack: CubeStack] RETURNS [BOOL] ~ {
RETURN[cubeStack.bottom = cubeStack.top];
};

LengthenCubeStack: PUBLIC PROC [cubeStack: CubeStack, amount: REAL Ź 1.3]
RETURNS [newStack: CubeStack]
~ {
newStack Ź NEW[CubeStackRep[Real.Round[amount*cubeStack.maxLength]]];
FOR i: INT IN [0..cubeStack.maxLength) DO newStack[i] Ź cubeStack[i]; ENDLOOP;
newStack.top Ź cubeStack.top;
newStack.bottom Ź cubeStack.bottom;
newStack.size Ź cubeStack.size;
};
<<Cross Support Procedures>>
nScratchCrossSequences: INT ~ 6;
scratchCrossSequences: ARRAY [0..nScratchCrossSequences) OF CrossSequence Ź ALL[NIL];

ObtainCrossSequence: PUBLIC ENTRY PROC RETURNS [crossSequence: CrossSequence] ~ {
FOR i: INT IN [0..nScratchCrossSequences) DO
crossSequence Ź scratchCrossSequences[i];
IF crossSequence = NIL THEN LOOP;
crossSequence.length Ź 0;
scratchCrossSequences[i] Ź NIL;
RETURN;
ENDLOOP;
crossSequence Ź NEW[CrossSequenceRep[100]];
};

ReleaseCrossSequence: PUBLIC ENTRY PROC [crossSequence: CrossSequence] ~ {
FOR i: INT IN [0..nScratchCrossSequences) DO
IF scratchCrossSequences[i] # NIL THEN LOOP;
scratchCrossSequences[i] Ź crossSequence;
RETURN;
ENDLOOP;
};
<<Ropes>>
NoSuchOctant:    PUBLIC ERROR [msg: ROPE] ~ CODE;
NoSuchFace:        PUBLIC ERROR [msg: ROPE] ~ CODE;

RopeFromOctant: PUBLIC PROC [octant: Octant] RETURNS [ROPE] ~ {
RETURN[Rope.Substr[IO.PutFR1["%g", IO.refAny[NEW[Octant Ź octant]]], 1]];
};

RopeFromFace: PUBLIC PROC [face: Face] RETURNS [ROPE] ~ {
RETURN[Rope.Substr[IO.PutFR1["%g", IO.refAny[NEW[Face Ź face]]], 1]];
};
<<>>
RopeFromDirection: PUBLIC PROC [direction: Direction] RETURNS [ROPE] ~ {
RETURN[Rope.Substr[IO.PutFR1["%g", IO.refAny[NEW[Direction Ź direction]]], 1]];
};

RopeFromDirectionSelect: PUBLIC PROC [directionSelect: DirectionSelect] RETURNS [ROPE] ~ {
RETURN[SELECT directionSelect FROM
x=>"x", y=>"y", z=>"z", xy=>"xy", xz=>"xz", yz=>"yz", xyz=>"xyz", ENDCASE=>"none"];
};

OctantFromRope: PUBLIC PROC [rope: ROPE] RETURNS [Octant] ~ {
FOR o: Octant IN Octant DO IF Equal[RopeFromOctant[o], rope] THEN RETURN[o]; ENDLOOP;
ERROR NoSuchOctant[rope];
};

FaceFromRope: PUBLIC PROC [rope: ROPE] RETURNS [Face] ~ {
FOR face: Face IN Face DO
IF Equal[RopeFromFace[face], rope] THEN RETURN[face];
ENDLOOP;
ERROR NoSuchFace[rope];
};

DirectionFromRope: PUBLIC PROC [rope: ROPE] RETURNS [Direction] ~ {
FOR d: Direction IN Direction DO
IF Equal[RopeFromDirection[d], rope] THEN RETURN[d];
ENDLOOP;
RETURN[none];
};

RopeFromRefAny: PUBLIC PROC [ref: REF ANY] RETURNS [rope: ROPE] ~ {
ioRope: ROPE ~ IO.PutFR1["%g", IO.refAny[ref]];
text: REF TEXT ~ RefText.AppendRope[RefText.ObtainScratch[90], ioRope, 1];
IF text[0] IN ['a..'z] THEN text[0] Ź text[0]-('a-'A);
rope Ź Rope.FromRefText[text];
RefText.ReleaseScratch[text];
};

Equal: PROC [r1, r2: ROPE] RETURNS [BOOL] ~ {RETURN[Rope.Equal[r1, r2, FALSE]]};

END.