Types and Constants
Triple: TYPE ~ G3dBasic.Triple;
Plane: TYPE ~ G3dPlane.Plane;
IntTriple: TYPE ~ G3dTriangle.IntTriple;
Segment: TYPE ~ G3dTriangle.Segment;
Segments: TYPE ~ G3dTriangle.Segments;
Triangle: TYPE ~ G3dTriangle.Triangle;
Triangles: TYPE ~ G3dTriangle.Triangles;
Cube: TYPE ~ ImplicitCubeTet.Cube;
Corner: TYPE ~ ImplicitCubeTet.Corner;
EventType: TYPE ~ ImplicitCubeTet.EventType;
Event: TYPE ~ ImplicitCubeTet.Event;
Events: TYPE ~ ImplicitCubeTet.Events;
TetEdge: TYPE ~ ImplicitCubeTet.TetEdge;
TetFace: TYPE ~ ImplicitCubeTet.TetFace;
TetIntersection: TYPE ~ ImplicitCubeTet.TetIntersection;
TetIntersections: TYPE ~ ImplicitCubeTet.TetIntersections;
TetIntersectProc: TYPE ~ ImplicitCubeTet.TetIntersectProc;
TetProc: TYPE ~ ImplicitCubeTet.TetProc;
Tetrahedron: TYPE ~ ImplicitCubeTet.Tetrahedron;
ThreeEdges: TYPE ~ ImplicitCubeTet.ThreeEdges;
ThreeCorners: TYPE ~ ImplicitCubeTet.ThreeCorners;
TwoCorners:
TYPE ~ ImplicitCubeTet.TwoCorners;
sqRt3: REAL ¬ RealFns.SqRt[3.0];
Tetrahedra
IntersectWithSurface:
PUBLIC
PROC [tet:
REF Tetrahedron, action: TetIntersectProc] ~ {
Pos: PROC [c: REF Corner, i: INT] RETURNS [j: INT] ~ {j ¬ IF c.state = positive THEN i ELSE 0};
[Artwork node; type 'Artwork on' to command tool]
SELECT (Pos[tet.a, 8]+Pos[tet.b, 4]+Pos[tet.c, 2]+Pos[tet.d, 1])
FROM
1 => action[3, bd, cd, ad, ab];
2 => action[3, ac, cd, bc, ab];
3 => action[4, ad, bd, bc, ac];
4 => action[3, ab, bc, bd, ab];
5 => action[4, ad, ab, bc, cd];
6 => action[4, ab, ac, cd, bd];
7 => action[3, ab, ac, ad, ab];
8 => action[3, ab, ad, ac, ab];
9 => action[4, ab, bd, cd, ac];
10 => action[4, ab, ad, cd, bc];
11 => action[3, ab, bd, bc, ab];
12 => action[4, ad, ac, bc, bd];
13 => action[3, cd, ac, bc, ab];
14 => action[3, bd, ad, cd, ab];
ENDCASE; -- 0, 15
};
FaceName:
PUBLIC
PROC [f: TetFace]
RETURNS [r:
ROPE] ~ {
r ¬ SELECT f FROM abc=>"abc", abd=>"abd", acd=>"acd", ENDCASE=>"bcd";
};
EdgeName:
PUBLIC
PROC [e: TetEdge]
RETURNS [r:
ROPE] ~ {
r ¬ SELECT e FROM ab=>"ab",ac=>"ac",ad=>"ad",bc=>"bc",bd=>"bd",ENDCASE=>"cd";
};
CornerName:
PUBLIC
PROC [c:
REF Corner, tet:
REF Tetrahedron]
RETURNS [r:
ROPE] ~ {
r ¬ SELECT c FROM tet.a => "a", tet.b => "b", tet.c => "c", tet.d => "d", ENDCASE => "?";
};
MakeTetrahedron:
PUBLIC
PROC
[
a, b, c, d: REF Corner,
which: INT,
scratch: REF Tetrahedron ¬ NIL]
RETURNS [tet: REF Tetrahedron]
~ {
IF (tet ¬ scratch) = NIL THEN tet ¬ NEW[Tetrahedron];
tet.which ¬ which;
tet.a ¬ a; tet.b ¬ b; tet.c ¬ c; tet.d ¬ d;
tet.triangles[abc] ¬ G3dTriangle.MakeTriangle[a.pos, c.pos, b.pos, tet.triangles[abc]];
tet.triangles[abd] ¬ G3dTriangle.MakeTriangle[a.pos, b.pos, d.pos, tet.triangles[abd]];
tet.triangles[acd] ¬ G3dTriangle.MakeTriangle[a.pos, d.pos, c.pos, tet.triangles[acd]];
tet.triangles[bcd] ¬ G3dTriangle.MakeTriangle[b.pos, c.pos, d.pos, tet.triangles[bcd]];
};
TrianglesInTetrahedron:
PUBLIC
PROC [
candidates: REF Triangles,
tet: REF Tetrahedron,
scratch: REF Triangles]
RETURNS [ret: REF Triangles]
~ {
IF (ret ¬ scratch) # NIL THEN ret.length ¬ 0;
FOR i:
INT
IN [0..candidates.length)
DO
t: REF Triangle ¬ candidates[i];
IF TriangleInTetrahedron[t, tet] THEN ret ¬ G3dTriangle.AddTriangle[t, ret];
ENDLOOP;
};
TriangleInTetrahedron:
PUBLIC
PROC [t:
REF Triangle, tet:
REF Tetrahedron]
RETURNS [BOOL]
~ {
This should be as fast as possible.
PlaneInTetrahedron:
PROC [p: Plane]
RETURNS [b:
BOOL] ~ {
Pos:
PROC [c:
REF Corner, type: CornerType]
RETURNS [b:
BOOL] ~ {
d: REAL ¬ c.pos.x*p.x+c.pos.y*p.y+c.pos.z*p.z+p.w;
cInfo[type] ¬ [TRUE, d, b ¬ d > 0.0, c];
};
s: BOOL ¬ Pos[tet.a, a];
b ¬ Pos[tet.b, b] # s OR Pos[tet.c, c] # s OR Pos[tet.d, d] # s;
};
Test if plane passes through cube:
CornerInfo: TYPE ~ RECORD [set: BOOL, distance: REAL, pos: BOOL, c: REF Corner];
CornerType: TYPE ~ {a, b, c, d};
cInfo: ARRAY CornerType OF CornerInfo ¬ ALL[[FALSE, 0.0, FALSE, NIL]];
IF PlaneInTetrahedron[t.plane]
THEN {
Code: TYPE ~ ARRAY TetFace OF BOOL;
SetCodes:
PROC
RETURNS [
ARRAY
[0..3)
OF Code]
~
{
Set:
PROC
[t:
Triple]
RETURNS
[r:
Code]
~
{
FOR face: TetFace
IN TetFace
DO
plane: Plane ¬ tet.triangles[face].plane; -- planes point inwards
r[face] ¬ t.x*plane.x+t.y*plane.y+t.z*plane.z+plane.w < 0.0;
ENDLOOP;
};
RETURN[[Set[t.p1], Set[t.p2], Set[t.p3]]];
};
NoneSet:
PROC [c: Code]
RETURNS [b:
BOOL ¬
TRUE] ~ {
FOR face: TetFace IN TetFace DO IF c[face] THEN RETURN[FALSE]; ENDLOOP;
};
Cross:
PROC
[p1, p2: Triple]
RETURNS
[b:
BOOL ¬
FALSE] ~
{
FOR face: TetFace
IN TetFace
DO
IF G3dTriangle.IntersectTriangle[p1, p2, tet.triangles[face]].intersect
THEN RETURN[TRUE];
ENDLOOP;
};
SetInfo:
PROC [c:
REF Corner]
RETURNS [i: CornerInfo] ~ {
d: REAL ¬ c.pos.x*t.plane.x+c.pos.y*t.plane.y+c.pos.z*t.plane.z+t.plane.w;
RETURN[[TRUE, d, d > 0.0, c]];
};
Hit:
PROC [t1, t2: CornerType]
RETURNS [b:
BOOL] ~ {
c1: CornerInfo ¬ cInfo[t1];
c2: CornerInfo ¬ cInfo[t2];
a: REAL ¬ c1.distance/(c1.distance-c2.distance);
b ¬ G3dTriangle.InsideTriangle[InterpCorners[a, c1.c, c2.c], t];
};
Test if entire triangle beyond either of tetrahedron's four faces:
out: ARRAY [0..3) OF Code ¬ SetCodes[];
FOR face: TetFace
IN TetFace
DO
IF out[0][face] AND out[1][face] AND out[2][face] THEN RETURN[FALSE];
ENDLOOP;
Test if any point inside tetrahedron:
FOR i: INT IN [0..3) DO IF NoneSet[out[i]] THEN RETURN[TRUE]; ENDLOOP;
Test if any triangle edge penetrates a tetrahedron face:
IF Cross[t.p1, t.p2] OR Cross[t.p2, t.p3] OR Cross[t.p3, t.p1] THEN RETURN[TRUE];
Test if any tetrahedron edge penetrates the triangle (note, at least a and b already set):
-- ab -- IF cInfo[a].pos # cInfo[b].pos AND Hit[a, b] THEN RETURN[TRUE];
-- ac -- IF NOT cInfo[c].set THEN cInfo[c] ¬ SetInfo[tet.c];
IF cInfo[a].pos # cInfo[c].pos AND Hit[a, c] THEN RETURN[TRUE];
-- ad -- IF NOT cInfo[d].set THEN cInfo[d] ¬ SetInfo[tet.d];
IF cInfo[a].pos # cInfo[d].pos AND Hit[a, d] THEN RETURN[TRUE];
-- bc -- IF cInfo[b].pos # cInfo[c].pos AND Hit[b, c] THEN RETURN[TRUE];
-- bd -- IF cInfo[b].pos # cInfo[d].pos AND Hit[b, d] THEN RETURN[TRUE];
-- cd -- IF cInfo[c].pos # cInfo[d].pos AND Hit[c, d] THEN RETURN[TRUE];
};
RETURN[FALSE];
};
InsideTetrahedron:
PUBLIC
PROC [t:
REF Tetrahedron, p: Triple]
RETURNS [
BOOL] ~ {
Test:
PROC
[face:
TetFace]
RETURNS
[
BOOL]
~
{
RETURN[G3dPlane.Side[p, t.triangles[face].plane] = positive];
};
RETURN[Test[abc] AND Test[abd] AND Test[acd] AND Test[bcd]];
};
EdgesFromFace:
PUBLIC
PROC [f: TetFace]
RETURNS [r: ThreeEdges] ~ {
r ¬
SELECT f
FROM
abc => [ab, ac, bc], abd => [ab, ad, bd], acd => [ac, ad, cd], ENDCASE => [bc, bd, cd];
};
CornersFromFace:
PUBLIC
PROC
[f:
TetFace,
t:
REF
Tetrahedron]
RETURNS
[r:
ThreeCorners]
~
{
r ¬
SELECT f
FROM
abc => [t.a, t.b, t.c], abd => [t.a, t.b, t.d], acd => [t.a, t.c, t.d], ENDCASE => [t.b, t.c, t.d];
};
CornersFromEdge:
PUBLIC
PROC [e:
TetEdge, t:
REF Tetrahedron]
RETURNS
[r:
TwoCorners]
~
{
r ¬
SELECT e
FROM
ab => [t.a, t.b], ac => [t.a, t.c], ad => [t.a, t.d],
bc => [t.b, t.c], bd => [t.b, t.d], ENDCASE => [t.c, t.d];
};
MidEdge:
PUBLIC
PROC [tet:
REF Tetrahedron, e: TetEdge]
RETURNS [Triple] ~ {
tc: TwoCorners ¬ CornersFromEdge[e, tet];
RETURN[MidCorner[tc.a, tc.b]];
};
LRFaces:
PUBLIC
PROC [e: TetEdge]
RETURNS [l, r: TetFace] ~ {
RETURN[LFace[e], RFace[e]]};
LFace:
PUBLIC
PROC [e: TetEdge]
RETURNS [f: TetFace] ~ {
f ¬ SELECT e FROM ab=>abd, ac=>abc, ad=>acd, bc=>bcd, bd=>abd, ENDCASE=>bcd;
};
RFace:
PUBLIC
PROC [e: TetEdge]
RETURNS [f: TetFace] ~ {
f ¬ SELECT e FROM ab=>abc, ac=>acd, ad=>abd, bc=>abc, bd=>bcd, ENDCASE=>acd;
};
Cubes
Decompose:
PUBLIC
PROC [c:
REF Cube, action: TetProc, scratch:
REF Tetrahedron ¬
NIL]
~
{
odd: BOOL ¬ (c.id.i+c.id.j+c.id.k) MOD 2 # 0;
FOR i: INT IN [0..5) DO IF NOT action[InnerGetTetN[c, odd, i, scratch]] THEN EXIT; ENDLOOP;
};
For those desiring an aid to the imagination, here are cut-outs for the tetrahedra that result
from the decomposition; the first is the single equilateral tetrahedra in the center of the cube,
the second represents the four corner tetrahedra. Try scaling by five.
[Artwork node; type 'Artwork on' to command tool]
GetTetN:
PUBLIC
PROC [c:
REF Cube, which:
INT, scratch:
REF Tetrahedron ¬
NIL]
RETURNS [t: REF Tetrahedron]
~ {
t ¬ InnerGetTetN[c, (c.id.i+c.id.j+c.id.k) MOD 2 # 0, which, scratch];
};
InnerGetTetN:
PROC
[cube:
REF
Cube,
odd:
BOOL,
which:
INT,
scratch:
REF
Tetrahedron
¬
NIL]
RETURNS [tet: REF Tetrahedron]
~ {
Inner:
PROC [a, b, c, d:
REF Corner] ~ {
Note: b, c, d to appear cw viewed from a.
tet ¬ MakeTetrahedron[a, b, c, d, which, scratch];
tet.cube ¬ cube;
};
IF odd
THEN SELECT which FROM
ELSE
SELECT which
FROM
[Artwork node; type 'Artwork on' to command tool]
0 => Inner[cube.lbf, cube.rbf, cube.ltf, cube.lbn];
1 => Inner[cube.ltn, cube.ltf, cube.rtn, cube.lbn];
2 => Inner[cube.rbn, cube.rbf, cube.lbn, cube.rtn];
3 => Inner[cube.rtf, cube.rtn, cube.ltf, cube.rbf];
4 => Inner[cube.lbn, cube.ltf, cube.rtn, cube.rbf];
ENDCASE => ERROR;
};
IDFromPoint:
PUBLIC
PROC [p: Triple, size:
REAL, offset: Triple ¬ []]
RETURNS [id: IntTriple]
~ {
id ¬ [Real.Floor[p.x+0.5], Real.Floor[p.y+0.5], Real.Floor[p.z+0.5]];
id ¬ [Real.Round[p.x/size], Real.Round[p.y/size], Real.Round[p.z/size]];
id ¬ [Real.Round[(p.x-offset.x)/size+0.5],
Real.Round[(p.y-offset.y)/size+0.5],
Real.Round[(p.z-offset.z)/size+0.5]];
};
PointFromID:
PUBLIC
PROC [id: IntTriple, size:
REAL]
RETURNS [p: Triple] ~ {
p ¬ [id.i*size, id.j*size, id.k*size];
};
NewCube:
PROC
[id:
IntTriple,
size:
REAL,
offset: Triple ¬ [], scratch:
REF
Cube
¬
NIL]
RETURNS [c: REF Cube]
~ {
c ¬ IF scratch # NIL THEN scratch ELSE NEW[Cube];
c.id ¬ id;
c.diag ¬ size*sqRt3;
c.l ¬ offset.x+(REAL[id.i]-0.5)*size; c.r ¬ c.l+size;
c.b ¬ offset.y+(REAL[id.j]-0.5)*size; c.t ¬ c.b+size;
c.n ¬ offset.z+(REAL[id.k]-0.5)*size; c.f ¬ c.n+size;
};
SetCube:
PUBLIC
PROC [
id: IntTriple,
lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf: REF Corner,
offset: Triple ¬ [],
scratch: REF Cube ¬ NIL]
RETURNS [c: REF Cube]
~ {
c ¬ NewCube[id, rbn.pos.x-lbn.pos.x, offset, scratch];
c.lbn ¬ lbn; c.lbf ¬ lbf; c.ltn ¬ ltn; c.ltf ¬ ltf;
c.rbn ¬ rbn; c.rbf ¬ rbf; c.rtn ¬ rtn; c.rtf ¬ rtf;
};
MakeCube:
PUBLIC
PROC
[id:
IntTriple,
size:
REAL,
offset: Triple ¬ [], scratch:
REF Cube
¬
NIL]
RETURNS [c: REF Cube]
~ {
Set:
PROC [x, y, z:
REAL, i, j, k:
INT]
RETURNS [c:
REF Corner] ~ {
c ¬ NEW[Corner ¬ [[i, j, k], [x, y, z], 0, unset, FALSE]];
};
c ¬ NewCube[id, size, offset, scratch];
c.lbn ¬ Set[c.l, c.b, c.n, id.i, id.j, id.k];
c.lbf ¬ Set[c.l, c.b, c.f, id.i, id.j, id.k+1];
c.ltn ¬ Set[c.l, c.t, c.n, id.i, id.j+1, id.k];
c.ltf ¬ Set[c.l, c.t, c.f, id.i, id.j+1, id.k+1];
c.rbn ¬ Set[c.r, c.b, c.n, id.i+1, id.j, id.k];
c.rbf ¬ Set[c.r, c.b, c.f, id.i+1, id.j, id.k+1];
c.rtn ¬ Set[c.r, c.t, c.n, id.i+1, id.j+1, id.k];
c.rtf ¬ Set[c.r, c.t, c.f, id.i+1, id.j+1, id.k+1];
};
TriangleInCube:
PUBLIC PROC [t:
REF Triangle, c:
REF Cube]
RETURNS [
BOOL] ~ {
This should be as fast as possible.
PlaneInCube:
PROC [p: Plane]
RETURNS [in:
BOOL] ~ {
Pos:
PROC [c:
REF Corner, type: CornerType]
RETURNS [b:
BOOL] ~ {
d: REAL ¬ c.pos.x*p.x+c.pos.y*p.y+c.pos.z*p.z+p.w;
cInfo[type] ¬ [TRUE, d, b ¬ d > 0.0, c];
};
s: BOOL ¬ Pos[c.lbn, lbn];
in ¬
IF
ABS[cInfo[lbn].distance] > c.diag
THEN FALSE
ELSE Pos[c.lbf, lbf] # s OR Pos[c.ltn, ltn] # s OR Pos[c.ltf, ltf] # s OR
Pos[c.rbn, rbn] # s OR Pos[c.rbf, rbf] # s OR Pos[c.rtn, rtn] # s OR
Pos[c.rtf, rtf] # s;
};
Test if plane passes through cube:
CornerInfo: TYPE ~ RECORD [set: BOOL, distance: REAL, pos: BOOL, c: REF Corner];
CornerType: TYPE ~ {lbn, lbf, ltn, ltf, rbn, rbf, rtn, rtf};
cInfo: ARRAY CornerType OF CornerInfo ¬ ALL[[FALSE, 0.0, FALSE, NIL]];
IF PlaneInCube[t.plane]
THEN {
[Artwork node; type 'Artwork on' to command tool]
Face: TYPE ~ {l, r, b, t, n, f};
Code: TYPE ~ ARRAY Face OF BOOL;
SetCodes:
PROC
RETURNS [
ARRAY
[0..3)
OF Code]
~
{
Set: PROC [t: Triple] RETURNS [r: Code] ~ {r¬[t.x<c.l,t.x>c.r,t.y<c.b,t.y>c.t,t.z<c.n,t.z>c.f]};
RETURN[[Set[t.p1], Set[t.p2], Set[t.p3]]];
};
NoneSet:
PROC [c: Code]
RETURNS [b:
BOOL ¬
TRUE] ~ {
FOR face: Face IN Face DO IF c[face] THEN RETURN[FALSE]; ENDLOOP;
};
SetInfo:
PROC [c:
REF Corner]
RETURNS [i: CornerInfo] ~ {
d: REAL ¬ c.pos.x*t.plane.x+c.pos.y*t.plane.y+c.pos.z*t.plane.z+t.plane.w;
RETURN[[TRUE, d, d > 0.0, c]];
};
Hit:
PROC [t1, t2: CornerType]
RETURNS [b:
BOOL] ~ {
c1: CornerInfo ¬ cInfo[t1];
c2: CornerInfo ¬ cInfo[t2];
a: REAL ¬ c1.distance/(c1.distance-c2.distance);
b ¬ G3dTriangle.InsideTriangle[InterpCorners[a, c1.c, c2.c], t];
};
out: ARRAY [0..3) OF Code ¬ SetCodes[];
Test if entire triangle beyond either of cube's six faces:
FOR face: Face
IN Face
DO
IF out[0][face] AND out[1][face] AND out[2][face] THEN RETURN[FALSE];
ENDLOOP;
Test if any point inside cube:
FOR i: INT IN [0..3) DO IF NoneSet[out[i]] THEN RETURN[TRUE]; ENDLOOP;
Test if any triangle edge penetrates a cube face:
FOR i:
INT
IN [0..3)
DO
Cross:
PROC
[face: Face, e,
e1,
e2,
a1,
a2,
amn,
amx,
b1,
b2,
bmn,
bmx:
REAL]
RETURNS [b: BOOL ¬ FALSE] ~ {
IF out[i][face] # out[j][face]
THEN {
a: REAL ¬ (e-e1)/(e2-e1);
IF amn > amx OR bmn > bmx THEN ERROR;
b ¬ (a1+a*(a2-a1)) IN [amn..amx] AND (b1+a*(b2-b1)) IN [bmn..bmx];
};
};
p, q: Triple;
j: INT ¬ (i+1) MOD 3;
SELECT i
FROM
0 => {p ¬ t.p1; q ¬ t.p2};
1 => {p ¬ t.p2; q ¬ t.p3};
ENDCASE => {p ¬ t.p3; q ¬ t.p1};
Note: l < r, b < t, n < f
IF Cross[
l, c.
l, p.
x, q.
x, p.
y, q.
y, c.
b, c.
t, p.
z, q.
z, c.
n, c.
f]
OR
Cross[r, c.r, p.x, q.x, p.y, q.y, c.b, c.t, p.z, q.z, c.n, c.f] OR
Cross[b, c.b, p.y, q.y, p.x, q.x, c.l, c.r, p.z, q.z, c.n, c.f] OR
Cross[t, c.t, p.y, q.y, p.x, q.x, c.l, c.r, p.z, q.z, c.n, c.f] OR
Cross[n, c.n, p.z, q.z, p.x, q.x, c.l, c.r, p.y, q.y, c.b, c.t] OR
Cross[
f, c.
f, p.
z, q.
z, p.
x, q.
x, c.
l, c.
r, p.
y, q.
y, c.
b, c.
t]
THEN RETURN[TRUE];
ENDLOOP;
Test if any cube edge penetrates the triangle (note, at least lbn and lbf already set):
-- lb -- IF cInfo[lbn].pos # cInfo[lbf].pos AND Hit[lbn, lbf] THEN RETURN[TRUE];
-- ln -- IF NOT cInfo[ltn].set THEN cInfo[ltn] ¬ SetInfo[c.ltn];
IF cInfo[lbn].pos # cInfo[ltn].pos AND Hit[lbn, ltn] THEN RETURN[TRUE];
-- lt -- IF NOT cInfo[ltf].set THEN cInfo[ltf] ¬ SetInfo[c.ltf];
IF cInfo[ltn].pos # cInfo[ltf].pos AND Hit[ltn, ltf] THEN RETURN[TRUE];
-- lf -- IF cInfo[ltf].pos # cInfo[lbf].pos AND Hit[ltf, lbf] THEN RETURN[TRUE];
-- rb -- IF NOT cInfo[rbn].set THEN cInfo[rbn] ¬ SetInfo[c.rbn];
IF NOT cInfo[rbf].set THEN cInfo[rbf] ¬ SetInfo[c.rbf];
IF cInfo[rbn].pos # cInfo[rbf].pos AND Hit[rbn, rbf] THEN RETURN[TRUE];
-- rf -- IF NOT cInfo[rtf].set THEN cInfo[rtf] ¬ SetInfo[c.rtf];
IF cInfo[rtf].pos # cInfo[rbf].pos AND Hit[rtf, rbf] THEN RETURN[TRUE];
-- rt -- IF NOT cInfo[rtn].set THEN cInfo[rtn] ¬ SetInfo[c.rtn];
IF cInfo[rtn].pos # cInfo[rtf].pos AND Hit[rtn, rtf] THEN RETURN[TRUE];
-- rn -- IF cInfo[rtn].pos # cInfo[rbn].pos AND Hit[rtn, rbn] THEN RETURN[TRUE];
-- bn -- IF cInfo[lbn].pos # cInfo[rbn].pos AND Hit[lbn, rbn] THEN RETURN[TRUE];
-- bf -- IF cInfo[lbf].pos # cInfo[rbf].pos AND Hit[lbf, rbf] THEN RETURN[TRUE];
-- tn -- IF cInfo[ltn].pos # cInfo[rtn].pos AND Hit[ltn, rtn] THEN RETURN[TRUE];
-- tf -- IF cInfo[ltf].pos # cInfo[rtf].pos AND Hit[ltf, rtf] THEN RETURN[TRUE];
};
RETURN[FALSE];
};
TrianglesInCube:
PUBLIC
PROC [
candidates: REF Triangles,
cube: REF Cube,
scratch: REF Triangles ¬ NIL]
RETURNS [ret: REF Triangles]
~ {
IF (ret ¬ scratch) # NIL THEN ret.length ¬ 0;
FOR i:
INT
IN [0..candidates.length)
DO
t: REF Triangle ¬ candidates[i];
IF TriangleInCube[t, cube] THEN ret ¬ G3dTriangle.AddTriangle[t, ret];
ENDLOOP;
};
Boundary Procedures
Side: TYPE ~ {in, on, out};
EndState:
TYPE ~
RECORD [side: Side, face: TetFace];
EventName:
PUBLIC
PROC [e: Event]
RETURNS [r: Rope.
ROPE] ~ {
r ¬ SELECT e.type FROM entry=> "entry", exit=> "exit", inside=> "inside", ENDCASE=> "graze";
};
scratchEvents:
ARRAY [0..10)
OF
REF Events ¬
ALL[
NIL];
ObtainEvents:
ENTRY
PROC RETURNS [
REF Events] ~ {
FOR i:
INT
IN [0..10)
DO
s: REF Events ¬ scratchEvents[i];
IF s # NIL THEN {scratchEvents[i] ¬ NIL; s.length ¬ 0; RETURN[s]};
ENDLOOP;
RETURN[NEW[Events[12]]];
};
ReleaseEvents:
ENTRY
PROC [scratch:
REF Events] ~ {
FOR i:
INT
IN [0..10)
DO
IF scratchEvents[i] = NIL THEN {scratchEvents[i] ¬ scratch; EXIT};
ENDLOOP;
};
AddEvent:
PROC [e: Event, events:
REF Events]
RETURNS [r:
REF Events] ~ {
LengthenEvents:
PROC
[events:
REF Events]
RETURNS
[new:
REF Events] ~
{
newLength: NAT ¬ MAX[Real.Ceiling[1.3*events.maxLength], 3];
new ¬ NEW[Events[newLength]];
FOR i: NAT IN [0..events.length) DO new[i] ¬ events[i]; ENDLOOP;
new.length ¬ events.length;
};
IF (r ¬ events) = NIL THEN r ¬ NEW[Events[1]];
IF r.length = r.maxLength THEN r ¬ LengthenEvents[r];
r[r.length] ¬ e;
r.length ¬ r.length+1;
};
GetEvents:
PUBLIC PROC
[
tet: REF Tetrahedron,
segments: REF Segments,
epsilon: REAL,
insideToo: BOOL,
scratch: REF Events ¬ NIL]
RETURNS [events: REF Events]
~ {
Info: TYPE ~ RECORD [nInts: INT, p1, p2: Triple, f1, f2: TetFace];
GetEndState:
PROC [p: Triple]
RETURNS [s: EndState ¬ [in, abc]] ~ {
near: BOOL ¬ FALSE;
nearFace: TetFace;
FOR l:
LIST
OF TetFace ¬
LIST[abc, abd, acd, bcd], l.rest
WHILE l #
NIL
DO
d: REAL ¬ G3dPlane.DistanceToPoint[p, tet.triangles[l.first].plane, FALSE];
IF d < -epsilon THEN RETURN[[out, abc]];
IF ABS[d] < epsilon THEN {near ¬ TRUE; nearFace ¬ l.first};
ENDLOOP;
RETURN[IF near THEN [on, nearFace] ELSE [in, abc]];
};
end1: EndState ¬ GetEndState[segments[0].p1];
IF (events ¬ scratch) # NIL THEN events.length ¬ 0;
FOR s:
INT
IN [0..segments.length)
DO
Eq: PROC [p, q: Triple] RETURNS [b: BOOL] ~ {b ¬ G3dVector.Distance[p, q] < epsilon};
Add:
PROC [type: EventType, p: Triple, face: TetFace ¬ abc] ~ {
e: Event ¬ [p, type, face, s];
IF events #
NIL
AND events.length > 0
AND Eq[p, events[events.length-1].p]
THEN {IF events[events.length-1].type = inside THEN events[events.length-1] ¬ e}
ELSE events ¬ AddEvent[e, events];
};
GetInfo:
PROC
RETURNS [i: Info] ~ {
Possibilities:
no intersections
one intersection (grazing), one intersection (entry), one intersection (exit)
two intersections (entry and exit)
GetAlpha:
PROC [f: TetFace] ~ {
t: REF Triangle ¬ tet.triangles[f];
d1: REAL ¬ G3dPlane.DistanceToPoint[seg.p1, t.plane];
d2: REAL ¬ G3dPlane.DistanceToPoint[seg.p2, t.plane];
IF (d1 > 0.0) # (d2 > 0.0)
THEN {
p: Triple ¬ G3dVector.Interp[a2 ¬ d1/(d1-d2), seg.p1, seg.p2];
IF NOT G3dTriangle.InsideTriangle[p, t] THEN RETURN;
IF (i.nInts ¬ i.nInts+1) = 1
THEN {a1 ¬ a2; i.p1 ¬ p; i.f1 ¬ f}
ELSE
IF a2 < a1
THEN {i.p2 ¬ i.p1; i.f2 ¬ i.f1; i.p1 ¬ p; i.f1 ¬ f}
ELSE {i.p2 ¬ p; i.f2 ¬ f};
IF i.nInts = 2 AND Eq[i.p1, i.p2] THEN i.nInts ¬ 1;
IF i.nInts = 2
AND (end1.side # out
OR end2.side # out)
THEN i.nInts ¬ 1;
(InsideTriangle[] fudge has yielded coincident intersections along tet edge)
};
};
a1, a2: REAL;
FOR f: TetFace IN TetFace DO IF i.nInts = 2 THEN EXIT ELSE GetAlpha[f]; ENDLOOP;
};
seg: Segment ¬ segments[s];
end2: EndState ¬ GetEndState[seg.p2];
i: Info ¬ GetInfo[];
SELECT i.nInts
FROM
0 =>
SELECT
TRUE
FROM
end1.side = out AND end2.side = out => NULL; -- completely out
end1.side = in AND end2.side = out,
end1.side = out AND end2.side = in => ERROR; -- should have nInt > 0
end1.side = in AND end2.side = in,
end1.side = on AND end2.side = in,
end1.side = in AND end2.side = on,
end1.side = on AND end2.side = on => Add[inside, seg.p1, end1.face];
end1.side = on AND end2.side = out => Add[exit, seg.p1, end1.face];
end1.side = out AND end2.side = on => Add[entry, seg.p2, end2.face];
ENDCASE;
1 =>
SELECT
TRUE
FROM
end1.side = in AND end2.side = in => ERROR; -- no ints, completely in
end1.side = out AND end2.side = out => Add[graze, i.p1, i.f1];
end1.side = out AND end2.side = on,
end1.side = out AND end2.side = in => Add[entry, i.p1, i.f1];
end1.side = on AND end2.side = in,
end1.side = in AND end2.side = on,
end1.side = on AND end2.side = on => Add[inside, seg.p1, end1.face];
end1.side = on AND end2.side = out,
end1.side = in
AND end2.side = out => {
Add[inside, seg.p1, end1.face];
Add[exit, i.p1, i.f1];
};
ENDCASE;
2 =>
SELECT
TRUE
FROM
end1.side = on AND end2.side = on => Add[inside, seg.p1, end1.face];
end1.side = in AND end2.side = in,
end1.side = in AND end2.side = out,
end1.side = out AND end2.side = in,
end1.side = on AND end2.side = in,
end1.side = in AND end2.side = on => ERROR; -- couldn't have 2 ints
end1.side = on
AND end2.side = out => {
Add[inside, seg.p1, end1.face];
Add[exit, i.p2, i.f2];
};
end1.side = out AND end2.side = on => Add[entry, i.p1, i.f1];
end1.side = out
AND end2.side = out => {
Add[entry, i.p1, i.f1];
Add[exit, i.p2, i.f2];
};
ENDCASE;
ENDCASE;
end1 ¬ end2;
ENDLOOP;
};
BoundaryIntersectTet:
PUBLIC PROC
[
tet: REF Tetrahedron,
segments: REF Segments,
insideToo: BOOL ¬ TRUE,
scratch: REF TetIntersections ¬ NIL]
RETURNS [REF TetIntersections]
~ {
Search all segments for intersections against all tetrahedral faces. Often this is not needed, but the illustration below, right shows anomolous cases (at patch corners, for example).
[Artwork node; type 'Artwork on' to command tool]
We assume, however, that the arbitrary nature of the boundary is constrained in one way. If there are no edge intersections, we place no constraints on the boundary (below, left), and any boundary-tetrahedron intersection point serves as entry and exit. If, however, there are edge intersections, we assume that the boundary intersects the tetrahedron in at most a single continuous section (below, right); here, the edge intersection must be incorporated into the boundary section according to the entry and exit points; we assume that these points are separated by the longest outside section (in this case, outside4).
[Artwork node; type 'Artwork on' to command tool]
Next: PROC [i: INT] RETURNS [INT] ~ {RETURN[(i+1) MOD events.length]};
maxD: REAL ¬ 0.0;
epsilon: REAL ¬ 0.01*G3dVector.Distance[tet.a.pos, tet.b.pos];
events: REF Events ¬ GetEvents[tet, segments, epsilon, insideToo, ObtainEvents[]];
mark, lastExit, start: INT ¬ -1;
IF scratch = NIL THEN scratch ¬ NEW[TetIntersections];
scratch.length ¬ 0;
IF events = NIL OR events.length = 0 THEN RETURN[scratch];
Find beginning of chain by finding longest gap from exit to entry:
FOR i:
INT
IN [0..events.length)
DO
type: EventType ¬ events[i].type;
IF type = graze OR type = exit THEN {mark ¬ lastExit ¬ i; EXIT};
ENDLOOP;
IF mark # -1
THEN
FOR i:
INT ¬ Next[mark], Next[i]
DO
CheckStart:
PROC ~ {
DistanceBetweenEntries:
PROC [e1, e2: Event]
RETURNS [d:
REAL] ~ {
d ¬ G3dVector.Distance[e1.p, segments[e1.seg].p2]+
G3dVector.Distance[e2.p, segments[e2.seg].p1];
FOR s: INT IN (e1.seg..e2.seg) DO d ¬ d+segments[s].length; ENDLOOP;
};
d: REAL ¬ DistanceBetweenEntries[events[lastExit], events[i]];
IF d > maxD THEN {maxD ¬ d; start ¬ i};
};
SELECT events[i].type
FROM
graze => {CheckStart[]; lastExit ¬ i};
entry => CheckStart[];
exit => lastExit ¬ i;
ENDCASE;
IF i = mark THEN EXIT;
ENDLOOP;
Begin at start and collect points all along chain:
IF start # -1
THEN
FOR i:
INT
IN [0..events.length)
DO
e: Event ¬ events[start];
scratch ¬ AddTetIntersection[scratch, [e.p, e.type = inside, e.face]];
start ¬ Next[start];
ENDLOOP;
ReleaseEvents[events];
RETURN[scratch];
};
AddTetIntersection:
PROC [seq:
REF TetIntersections, int: TetIntersection]
RETURNS [REF TetIntersections]
~ {
IF seq = NIL THEN seq ¬ NEW[TetIntersections[3]];
IF seq.length = seq.maxLength
THEN {
old: REF TetIntersections ¬ seq;
seq ¬ NEW[TetIntersections[MAX[Real.Ceiling[1.3*old.maxLength], 3]]];
FOR i: INT IN [0..old.length) DO seq[i] ¬ old[i]; ENDLOOP;
seq.length ¬ old.length;
};
seq[seq.length] ¬ int;
seq.length ¬ seq.length+1;
RETURN[seq];
};
BoundaryIntersectTet: PUBLIC PROC [
tet: REF Tetrahedron,
segments: REF Segments,
mode: TetIntersectMode ¬ all,
insideToo: BOOL ¬ TRUE,
f1, f2: TetFace ¬ abc,
scratch: REF TetIntersections ¬ NIL]
RETURNS [REF TetIntersections]
~ {
Info: TYPE ~ RECORD [p: Triple, f: TetFace, s: INT];
Infos: TYPE ~ RECORD [nInfos: [0..2] ¬ 0, infos: ARRAY [0..2) OF Info];
Eq: PROC [s, t: Triple] RETURNS [b: BOOL] ~ {b ¬ G3dVector.Distance[s, t] < epsilon};
GetInfos: PROC [start, stop: INT] RETURNS [Infos] ~ {
GetInfo: PROC [s: INT] RETURNS [i: Infos] ~ {
GetAlpha: PROC [f: TetFace] ~ {
t: REF Triangle ¬ tet.triangles[f];
d1: REAL ¬ G3dPlane.DistanceToPoint[seg.p1, t.plane];
d2: REAL ¬ G3dPlane.DistanceToPoint[seg.p2, t.plane];
IF (d1 > 0.0) # (d2 > 0.0) THEN {
p: Triple ¬ G3dVector.Interp[a2 ¬ d1/(d1-d2), seg.p1, seg.p2];
IF NOT Support.InsideTriangle[p, t] THEN RETURN;
IF (i.nInfos ¬ i.nInfos+1) = 1
THEN {
found ¬ TRUE;
save ¬ IF mode = two THEN (IF f = f1 THEN f2 ELSE f1) ELSE f;
a1 ¬ a2;
i.infos[0] ¬ [p, f];
}
ELSE IF a2 < a1
THEN {i.infos[1] ¬ i.infos[0]; i.infos[0] ¬ [p, f]}
ELSE i.infos[1] ¬ [p, f];
IF i.nInfos = 2 AND Eq[i.infos[0].p, i.infos[1].p] THEN {
InsideTriangle[] fudge has yielded coincident intersections along tet edge,
so don't play trick of limiting faces searched for exit point.
i.nInfos ¬ 1;
found ¬ FALSE;
};
};
};
a1, a2: REAL;
seg: Segment ¬ segments[s];
SELECT TRUE FROM
found AND mode # all => GetAlpha[save];
found AND mode = all =>
FOR f: TetFace IN TetFace DO IF f # save AND nInt < 2 THEN GetAlpha[f]; ENDLOOP;
mode = two => {GetAlpha[f1]; GetAlpha[f2]};
ENDCASE => FOR f: TetFace IN TetFace DO IF nInt < 2 THEN GetAlpha[f]; ENDLOOP;
};
start, stop: INT;
save: TetFace ¬ abc;
found: BOOL ¬ FALSE;
epsilon: REAL ¬ 0.01*G3dVector.Distance[tet.a.pos, tet.b.pos];
IF scratch = NIL THEN scratch ¬ ObtainTetIntersections[];
scratch.length ¬ 0;
FOR s: INT IN [0..segments.length) DO
infos: Infos ¬ GetInfo[s];
start ¬ stop ¬ s;
SELECT infos.nInfos FROM
1 => {
infoStart, infoStop: Info ¬ info1;
fwd: BOOL ¬ Support.InsideTetrahedron[tet, segments[s].p2];
FOR i: INT ¬ IF fwd THEN s+1 ELSE segments.length-1, IF fwd THEN i+1 ELSE i-1 DO
IF i NOT IN [0..INT[segments.length]) THEN EXIT;
IF GetInfo[i] = 0 THEN LOOP;
IF fwd THEN {stop ¬ i; infoStop ¬ info1} ELSE {start ¬ i; infoStart ¬ info1};
EXIT;
ENDLOOP;
IF start = stop OR Eq[infoStart.p, infoStop.p] THEN {
[] ¬ action[info1.p, graze, info1.f];
RETURN;
};
IF NOT action[infoStart.p, entry, infoStart.f] THEN RETURN;
IF insideToo THEN FOR i: INT IN [start..stop) DO
p: Triple ¬ segments[i].p2;
IF Eq[p,infoStart.p] OR Eq[p,infoStop.p] OR (i>start AND Eq[p,segments[i-1].p2])
THEN LOOP;
IF NOT action[segments[i].p2, inside, abc] THEN RETURN;
ENDLOOP;
[] ¬ action[infoStop.p, exit, infoStop.f];
};
2 => IF action[info1.p, entry, info1.f] THEN [] ¬ action[info2.p, exit, info2.f];
ENDCASE => LOOP;
EXIT;
ENDLOOP;
RETURN[scratch];
};
TetCenter:
PUBLIC PROC [tet:
REF Tetrahedron]
RETURNS [t: Triple] ~ {
t.x ¬ 0.25*(tet.a.pos.x+tet.b.pos.x+tet.c.pos.x+tet.d.pos.x);
t.y ¬ 0.25*(tet.a.pos.y+tet.b.pos.y+tet.c.pos.y+tet.d.pos.y);
t.z ¬ 0.25*(tet.a.pos.z+tet.b.pos.z+tet.c.pos.z+tet.d.pos.z);
};
AddToSegments:
PROC
[segments:
REF
Segments,
s:
Segment]
RETURNS
[
REF
Segments]
~
{
LengthenSegments:
PROC [segments:
REF Segments]
RETURNS [new:
REF Segments] ~ {
newLength: NAT ¬ MAX[Real.Ceiling[1.3*segments.maxLength], 3];
new ¬ NEW[Segments[newLength]];
FOR i: NAT IN [0..segments.length) DO new[i] ¬ segments[i]; ENDLOOP;
new.length ¬ segments.length;
};
IF segments = NIL THEN segments ¬ NEW[Segments[1]];
IF segments.length = segments.maxLength THEN segments ¬ LengthenSegments[segments];
segments[segments.length] ¬ s;
segments.length ¬ segments.length+1;
RETURN[segments];
};
OrderRandomBoundaryIndices: PROC [] ~ {
size: INT ¬ LengthOfBoundary[];
IF nRandomIndices > size/2
THEN {
do it the hard way: allocate a bit buffer, etc.
}
ELSE {
Action: PROC [i: INT] ~ {
IF i = 0 THEN hitZero ¬ TRUE;
IF i = size-1 THEN hitSize ¬ TRUE;
min1 ¬ MIN[min1, i];
max1 ¬ MAX[max1, i];
min2 ¬ MIN[min2, (i+size/2) MOD size];
max2 ¬ MAX[max2, (i+size/2) MOD size];
};
hitZero, hitSize: BOOL ¬ FALSE;
min1, min2: INT ¬ 1000000;
max1, max2: INT ¬ 0;
GetRandomBoundaryIndices[Action];
IF hitZero AND hitSize
THEN {
max2 ¬ (max2+size/2) MOD size;
min2 ¬ (min2+size/2) MOD size;
Process[min2, max2];
}
ELSE Process[min1, max1];
};
};
scratchTetIntersections:
ARRAY [0..10)
OF
REF TetIntersections ¬
ALL[
NIL];
ObtainTetIntersections:
PUBLIC ENTRY
PROC RETURNS [
REF TetIntersections] ~ {
FOR i:
INT
IN [0..10)
DO
s: REF TetIntersections ¬ scratchTetIntersections[i];
IF s # NIL THEN {scratchTetIntersections[i] ¬ NIL; s.length ¬ 0; RETURN[s]};
ENDLOOP;
RETURN[NEW[TetIntersections[12]]];
};
ReleaseTetIntersections:
PUBLIC ENTRY
PROC [scratch:
REF TetIntersections] ~ {
FOR i:
INT
IN [0..10)
DO
IF scratchTetIntersections[i] = NIL THEN {scratchTetIntersections[i] ¬ scratch; EXIT};
ENDLOOP;
};