Types
Triple: TYPE ~ G3dBasic.Triple;
Patch: TYPE ~ G3dPatch.Patch;
Plane:
TYPE ~ G3dPlane.Plane;
IntTriple: TYPE ~ G3dTriangle.IntTriple;
IntTriples: TYPE ~ G3dTriangle.IntTriples;
Segment: TYPE ~ G3dTriangle.Segment;
Segments: TYPE ~ G3dTriangle.Segments;
Triangle: TYPE ~ G3dTriangle.Triangle;
Triangles: TYPE ~ G3dTriangle.Triangles;
TriangleProc: TYPE ~ G3dTriangle.TriangleProc;
TriIntersection: TYPE ~ G3dTriangle.TriIntersection;
Triangle Operations
TrianglesFromPatch:
PUBLIC PROC [patch: Patch, sRes, tRes:
INT]
RETURNS [triangles: REF Triangles]
~ {
AddTriangle:
PROC [p1, p2, p3: Triple] ~ {
triangles[triangles.length] ¬ MakeTriangle[p1, p2, p3];
triangles.length ¬ triangles.length+1;
};
triangles ¬ NEW[Triangles[2*sRes*tRes]];
FOR s:
INT
IN [0..sRes)
DO
s0: REAL ¬ REAL[s]/REAL[sRes];
s1: REAL ¬ REAL[s+1]/REAL[sRes];
FOR t:
INT
IN [0..tRes)
DO
t0: REAL ¬ REAL[t]/REAL[tRes];
t1: REAL ¬ REAL[t+1]/REAL[tRes];
p00: Triple ¬ G3dPatch.Position[patch, s0, t0];
p01: Triple ¬ G3dPatch.Position[patch, s0, t1];
p10: Triple ¬ G3dPatch.Position[patch, s1, t0];
p11: Triple ¬ G3dPatch.Position[patch, s1, t1];
AddTriangle[p00, p01, p10];
AddTriangle[p10, p01, p11];
ENDLOOP;
ENDLOOP;
};
BoundaryFromPatch:
PUBLIC
PROC
[patch:
Patch,
sRes,
tRes:
INT,
cullColinearPoints:
BOOL
¬
FALSE]
RETURNS [segments: REF Segments]
~ {
oriented to agree with TrianglesFromPatch
AddSeg:
PROC [s, t:
REAL] ~ {
pp: Triple ¬ G3dPatch.Position[patch, s, t];
segments ¬ AddToSegments[segments, [p, pp, G3dVector.Distance[p, pp]]];
p ¬ pp;
};
p: Triple ¬ G3dPatch.Position[patch, 0.0, 0.0];
segments ¬ NEW[Segments[2*(sRes+tRes)]];
FOR t:
INT
IN [1..tRes]
DO
-- s = 0, t from 0 to 1
AddSeg[0.0, REAL[t]/REAL[tRes]];
ENDLOOP;
FOR s:
INT
IN [1..sRes]
DO
-- s from 0 to 1, t = 1
AddSeg[REAL[s]/REAL[sRes], 1.0];
ENDLOOP;
FOR t:
INT DECREASING
IN [0..tRes-1]
DO
-- t s = 1, from 1 to 0
AddSeg[1.0, REAL[t]/REAL[tRes]];
ENDLOOP;
FOR s:
INT
DECREASING
IN [0..sRes-1]
DO
-- s from 1 to 0, t = 0
AddSeg[REAL[s]/REAL[sRes], 0.0];
ENDLOOP;
IF cullColinearPoints THEN CullColinearPoints[segments];
};
CullColinearPoints:
PROC [segments:
REF Segments] ~ {
FOR i1:
INT ¬ 0, i1+1
WHILE i1 <
INT[segments.length]
DO
i2: INT ¬ (i1+1) MOD segments.length;
s1: Segment ¬ segments[i1];
s2: Segment ¬ segments[i2];
IF G3dVector.Collinear[s1.p1, s1.p2, s2.p2]
THEN {
segments[i1].p2 ¬ s2.p2;
segments.length ¬ segments.length-1;
FOR i: INT IN [i2..segments.length) DO segments[i] ¬ segments[i+1]; ENDLOOP;
};
ENDLOOP;
};
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];
};
DistanceToTriangles:
PUBLIC
PROC
[triangles:
REF
Triangles,
p:
Triple]
RETURNS [minD: REAL]
~ {
minAbsD, minSqD: REAL ¬ 100000.0;
FOR i:
INT
IN [0..triangles.length)
DO
t: REF Triangle ¬ triangles[i];
d: REAL ¬ G3dPlane.DistanceToPoint[p, t.plane, FALSE];
absD: REAL ¬ ABS[d];
IF absD < minAbsD
THEN {
q: Triple ~ G3dPlane.ProjectPointToPlane[p, t.plane];
IF InsideTriangle[q, t]
THEN {minAbsD ¬ absD; minD ¬ d; minSqD ¬ d*d}
ELSE {
Test:
PROC [p1, p2: Triple] ~ {
near: Triple ¬ G3dVector.NearestToSegment[p1, p2, p].point;
sqD: REAL ¬ G3dVector.SquareDistance[near, p];
IF sqD < minSqD
THEN {
minSqD ¬ sqD;
minAbsD ¬ RealFns.SqRt[sqD];
minD ¬ IF d < 0.0 THEN -minAbsD ELSE minAbsD;
};
};
Test[t.p1, t.p2]; Test[t.p2, t.p3]; Test[t.p3, t.p1];
};
};
ENDLOOP;
};
TriangleCenter:
PUBLIC
PROC
[t:
REF Triangle]
RETURNS
[c:
Triple]
~
{
c ¬ G3dVector.Midpoint[G3dVector.Midpoint[t.p1, t.p2], G3dVector.Midpoint[t.p2, t.p3]];
};
AverageTriangleNormals:
PUBLIC
PROC [triangles:
REF Triangles]
RETURNS [n: Triple ¬ []]
~ {
FOR i:
INT
IN [0..triangles.length)
DO
n ¬ G3dVector.Add[n, G3dPlane.NormalFromPlane[triangles[i].plane]];
ENDLOOP;
n ¬ G3dVector.Unit[G3dVector.Div[n, triangles.length]];
};
MaxDeviationFromNormal:
PUBLIC
PROC
[normal:
Triple,
triangles:
REF Triangles]
RETURNS [deg: REAL ¬ 1.0]
~ {
FOR i:
INT
IN [0..triangles.length)
DO
deg ¬ MIN[deg, G3dVector.Dot[normal, G3dPlane.NormalFromPlane[triangles[i].plane]]];
ENDLOOP;
deg ¬ G2dBasic.ArcCosDeg[deg];
};
ApplyToTrianglesSides:
PUBLIC
PROC [triangles:
REF Triangles, action: TriangleProc] ~ {
FOR i:
INT
IN [0..triangles.length)
DO
t: REF Triangle ¬ triangles[i];
IF
NOT action[t.p1, t.p2, t]
OR
NOT action[t.p2, t.p3, t]
OR
NOT action[t.p3, t.p1, t]
THEN EXIT;
ENDLOOP;
};
MakeTriangle:
PUBLIC
PROC [p1, p2, p3: Triple, scratch:
REF Triangle ¬
NIL]
RETURNS [t: REF Triangle]
~ {
t ¬ IF scratch = NIL THEN NEW[Triangle] ELSE scratch;
t.p1 ¬ p1; t.p2 ¬ p2; t.p3 ¬ p3;
SetTriangle[t];
};
SetTriangle:
PUBLIC
PROC [t:
REF Triangle] ~ {
q1, q2, q3: G3dBasic.Pair;
t.plane ¬ G3dPlane.FromThreePoints[t.p1, t.p2, t.p3, TRUE];
t.majorPlane ¬ G3dPlane.GetMajorPlane[t.plane];
q1 ¬ G3dPlane.ProjectPointToMajorPlane[t.p1, t.majorPlane];
q2 ¬ G3dPlane.ProjectPointToMajorPlane[t.p2, t.majorPlane];
q3 ¬ G3dPlane.ProjectPointToMajorPlane[t.p3, t.majorPlane];
IF (
SELECT t.majorPlane
FROM xy => t.plane.z, xz => t.plane.y,
ENDCASE => t.plane.x) > 0.0
THEN {
t.l1 ¬ G2dVector.Line[q1, q2];
t.l2 ¬ G2dVector.Line[q2, q3];
t.l3 ¬ G2dVector.Line[q3, q1];
}
ELSE {
-- invert from normal sense
t.l1 ¬ G2dVector.Line[q2, q1];
t.l2 ¬ G2dVector.Line[q3, q2];
t.l3 ¬ G2dVector.Line[q1, q3];
};
};
SetTriangles:
PUBLIC
PROC [triangles:
REF Triangles] ~ {
FOR i: INT IN [0..triangles.length) DO SetTriangle[triangles[i]]; ENDLOOP;
};
InsideTriangle:
PUBLIC
PROC [p: Triple, t:
REF Triangle]
RETURNS [b:
BOOL] ~ {
pair: G3dBasic.Pair ¬ G3dPlane.ProjectPointToMajorPlane[p, t.majorPlane];
Test for inclusion within three bounding lines;
provide some fudge so lines aimed at junctions of triangles will generate an intersection.
To see the problem, intersect the following segment with the following tetrahedron:
segment: [[0.5, -0.324, -0.324], [0.5, -0.418, -0.418]]
tetrahedral points: [0.45,-0.35,-0.35], [0.45,-0.35,-0.25], [0.45,-0.25,-0.35], [0.55,-0.35,-0.35]
b ¬ G2dVector.DistanceToLine[pair, t.l1] > -0.0001
AND
G2dVector.DistanceToLine[pair, t.l2] > -0.0001 AND
G2dVector.DistanceToLine[pair, t.l3] > -0.0001;
};
IntersectTriangle:
PUBLIC
PROC
[p1,
p2:
Triple,
t:
REF
Triangle]
RETURNS [i: TriIntersection]
~ {
d1: REAL ¬ p1.x*t.plane.x+p1.y*t.plane.y+p1.z*t.plane.z+t.plane.w;
d2: REAL ¬ p2.x*t.plane.x+p2.y*t.plane.y+p2.z*t.plane.z+t.plane.w;
i.intersect ¬
IF (d1 > 0) = (d2 > 0)
THEN FALSE
ELSE InsideTriangle[i.point ¬ G3dVector.Interp[d1/(d1-d2), p1, p2], t];
};
Normal:
PUBLIC PROC [p1, p2, p3: Triple]
RETURNS [n: Triple] ~ {
n ¬ G3dVector.Cross[G3dVector.Sub[p2, p1], G3dVector.Sub[p3, p1]];
};
scratchTriangles:
ARRAY [0..10)
OF
REF Triangles ¬
ALL[
NIL];
ObtainTriangles:
PUBLIC ENTRY
PROC RETURNS [
REF Triangles] ~ {
FOR i:
INT
IN [0..10)
DO
scratch: REF Triangles ¬ scratchTriangles[i];
IF scratch # NIL THEN {scratchTriangles[i] ¬ NIL; scratch.length ¬ 0; RETURN[scratch]};
ENDLOOP;
RETURN[NEW[Triangles[12]]];
};
ReleaseTriangles:
PUBLIC ENTRY
PROC [scratch:
REF Triangles] ~ {
FOR i:
INT
IN [0..10)
DO
IF scratchTriangles[i] = NIL THEN {scratchTriangles[i] ¬ scratch; EXIT};
ENDLOOP;
};
AddTriangle:
PUBLIC
PROC [t:
REF Triangle, triangles:
REF Triangles]
RETURNS [r: REF Triangles]
~ {
LengthenTriangles:
PROC
[triangles:
REF Triangles]
RETURNS
[new:
REF Triangles] ~
{
newLength: NAT ¬ MAX[Real.Ceiling[1.3*triangles.maxLength], 3];
new ¬ NEW[Triangles[newLength]];
FOR i: NAT IN [0..triangles.length) DO new[i] ¬ triangles[i]; ENDLOOP;
new.length ¬ triangles.length;
};
IF (r ¬ triangles) = NIL THEN r ¬ NEW[Triangles[1]];
IF r.length = r.maxLength THEN r ¬ LengthenTriangles[r];
r[r.length] ¬ t;
r.length ¬ r.length+1;
};
CopyTriangles:
PUBLIC
PROC [triangles:
REF Triangles]
RETURNS [out:
REF Triangles] ~ {
out ¬ NEW[Triangles[triangles.length]];
FOR i: INT IN [0..triangles.length) DO out[i] ¬ triangles[i]; ENDLOOP;
out.length ¬ triangles.length;
};
AddIntTriple:
PUBLIC
PROC [i: IntTriple, intTriples:
REF IntTriples]
RETURNS [r: REF IntTriples]
~ {
LengthenIntTriples:
PROC
[ints:
REF IntTriples]
RETURNS
[new:
REF IntTriples] ~
{
newLength: NAT ¬ MAX[Real.Ceiling[1.3*ints.maxLength], 3];
new ¬ NEW[IntTriples[newLength]];
FOR i: NAT IN [0..ints.length) DO new[i] ¬ ints[i]; ENDLOOP;
new.length ¬ ints.length;
};
IF (r ¬ intTriples) = NIL THEN r ¬ NEW[IntTriples[1]];
IF r.length = r.maxLength THEN r ¬ LengthenIntTriples[r];
r[r.length] ¬ i;
r.length ¬ r.length+1;
};