G3dTubeGeometryImpl.mesa
Copyright Ó 1985, 1992 by Xerox Corporation. All rights reserved.
Bloomenthal, July 15, 1992 4:11 pm PDT
Glassner, June 29, 1989 3:55:56 pm PDT
DIRECTORY Cubic2, G2dBasic, G2dContour, G2dSpline, G3dBasic, G3dCurve, G3dMatrix, G3dPlane, G3dSpline, G3dVector, MessageWindow, Real, RealFns, G3dTube, Vector2;
G3dTubeGeometryImpl: CEDAR PROGRAM
IMPORTS Cubic2, G2dBasic, G2dContour, G3dCurve, G3dMatrix, G3dPlane, G3dSpline, G3dVector, MessageWindow, Real, RealFns, G3dTube, Vector2
EXPORTS G3dTube
~ BEGIN
Types
Contour:     TYPE ~ G2dContour.Contour;
ContourRep:    TYPE ~ G2dContour.ContourRep;
ContourSequence:  TYPE ~ G2dContour.ContourSequence;
ContourSequenceRep: TYPE ~ G2dContour.ContourSequenceRep;
NearSpline2d:   TYPE ~ G2dSpline.NearSpline;
Pair:      TYPE ~ G3dBasic.Pair;
PairSequence:   TYPE ~ G3dBasic.PairSequence;
Ray:      TYPE ~ G3dBasic.Ray;
Triad:      TYPE ~ G3dBasic.Triad;
Triple:     TYPE ~ G3dBasic.Triple;
TripleSequence:   TYPE ~ G3dBasic.TripleSequence;
Curve:     TYPE ~ G3dCurve.Curve;
Matrix:     TYPE ~ G3dMatrix.Matrix;
SplineSequence:   TYPE ~ G3dSpline.SplineSequence;
Spline:     TYPE ~ G3dSpline.Spline;
Bezier:     TYPE ~ G3dSpline.Bezier;
XSection:     TYPE ~ G3dTube.XSection;
XSectionRep:    TYPE ~ G3dTube.XSectionRep;
XSectionProc:    TYPE ~ G3dTube.XSectionProc;
XSectionSequence:  TYPE ~ G3dTube.XSectionSequence;
XSectionSequenceRep: TYPE ~ G3dTube.XSectionSequenceRep;
RadiusMode:    TYPE ~ G3dTube.RadiusMode;
Tube:      TYPE ~ G3dTube.Tube;
TubePlace:    TYPE ~ G3dTube.TubePlace;
TubeProc:     TYPE ~ G3dTube.TubeProc;
TubeRep:     TYPE ~ G3dTube.TubeRep;
origin:     Triple ¬ [];
Creation
NewTube: PUBLIC PROC [
knots: TripleSequence,
tension: REAL ¬ 1.0,
tangent0, tangent1: Triple ¬ [],
cyclic: BOOL ¬ FALSE]
RETURNS [head: Tube ¬ NIL]
~ {
IF knots # NIL AND knots.length > 0 THEN {
SetBez: TubeProc ~ {
tube.bezier ¬ G3dSpline.BezierFromHermite[tube.p0, tube.p1, tube.v0, tube.v1];
};
MakeTube: PROC [p0: Triple, spline: Spline, parent: Tube] RETURNS [new: Tube] ~ {
v0: Triple ¬ G3dSpline.Velocity[spline, 0.0];
new ¬ NEW[TubeRep ¬ [
p0: p0,
v0: v0,
spline: spline,
prev: parent,
tens0: tension,
tens1: tension,
r0: 1.0,
r1: 1.0,
splineValid: TRUE
]];
IF parent # NIL THEN {parent.p1 ¬ p0; parent.v1 ¬ v0};
};
splines: SplineSequence ¬ IF cyclic
THEN G3dSpline.InterpolateCyclic[knots, tension]
ELSE G3dSpline.Interpolate[knots, tangent0, tangent1, tension];
tube: Tube ¬ head ¬ MakeTube[knots[0], splines[0], NIL];
FOR n: NAT IN [1..splines.length) DO 
tube.next ¬ MakeTube[knots[n], splines[n], tube]; 
tube ¬ tube.next;
ENDLOOP;
tube.p1 ¬ knots[knots.length-1];
tube.v1 ¬ G3dSpline.Velocity[splines[splines.length-1], 1.0];
G3dTube.ApplyToTube[head, SetBez];
};
};
Evaluation
Spot: TYPE ~ RECORD [tube: Tube, t: REAL];
GetSpot: PROC [tube: Tube, t: REAL] RETURNS [s: Spot] ~ {
nTubes: INTEGER ¬ G3dTube.NTubes[tube];
seg: INT ¬ IF t >= 1.0 THEN nTubes-1 ELSE Real.Floor[t*nTubes];
s.tube ¬ tube;
s.t ¬ IF t >= 1.0 THEN 1.0 ELSE t-REAL[seg]/REAL[nTubes];
FOR n: INT IN [0..seg) DO s.tube ¬ s.tube.next; ENDLOOP;
};
EvalPosition: PUBLIC PROC [tube: Tube, t: REAL] RETURNS [Triple] ~ {
spot: Spot ¬ GetSpot[tube, t];
RETURN[G3dSpline.Position[spot.tube.spline, spot.t]];
};
EvalXSection: PUBLIC PROC [tube: Tube, t: REAL] RETURNS [XSection] ~ {
spot: Spot ¬ GetSpot[tube, t];
IF NOT tube.xSectionsValid THEN MakeXSections[tube];
RETURN[GetXSection[spot.tube, spot.t]];
};
Splines
SetSpline: PUBLIC PROC [tube: Tube] ~ {
tubeProc: TubeProc ~ {SetSectionSpline[tube]};
G3dTube.ApplyToTube[tube, tubeProc];
};
SetSectionSpline: PUBLIC PROC [tube: Tube] ~ {
IF tube.splineValid THEN RETURN;
IF tube.v0 = [] THEN tube.v0 ¬ G3dVector.Sub[tube.p1, tube.p0];
IF tube.v1 = [] THEN tube.v1 ¬ tube.v0;
tube.spline ¬ G3dSpline.SplineFromHermite[
[tube.p0, tube.p1, G3dVector.Mul[tube.v0, tube.tens0], G3dVector.Mul[tube.v1, tube.tens1]],
tube.spline];
tube.bezier ¬ G3dSpline.BezierFromHermite[tube.p0, tube.p1, tube.v0, tube.v1];
tube.splineValid ¬ TRUE;
tube.xSectionsValid ¬ tube.lengthsValid ¬ FALSE;
};
ExtractPartTube: PUBLIC PROC [in: Tube, t0, t1: REAL, out: Tube ¬ NIL] RETURNS [Tube] ~ {
IF in # out THEN out ¬ G3dTube.CopyTubeSection[in, out];
out.spline ¬ G3dSpline.Reparameterize[in.spline, t0, t1, in.spline];
out.p0 ¬ G3dSpline.Position[in.spline, 0.0];
out.p1 ¬ G3dSpline.Position[in.spline, 1.0];
out.v0 ¬ G3dSpline.Velocity[in.spline, 0.0];
out.v1 ¬ G3dSpline.Velocity[in.spline, 1.0];
out.r0 ¬ Radius[in, t0];
out.r1 ¬ Radius[in, t1];
out.tens0 ¬ Tension[in, t0];
out.tens1 ¬ Tension[in, t1];
out.tw0 ¬ Twist[in, t0];
out.tw1 ¬ Twist[in, t1];
RETURN[out];
};
CurveFromTube: PUBLIC PROC [tube: Tube] RETURNS [Curve] ~ {
nSections: NAT ¬ 0;
splines: G3dSpline.SplineSequence;
FOR t: Tube ¬ tube, t.next WHILE t # NIL DO nSections ¬ nSections+1; ENDLOOP;
splines ¬ NEW[G3dSpline.SplineSequenceRep[nSections]];
FOR n: NAT IN [0..splines.length ¬ nSections) DO
splines[n] ¬ tube.spline;
tube ¬ tube.next;
ENDLOOP;
RETURN[G3dCurve.CurveFromSplines[splines]];
};
Epsilon, Scale, Tension, Twist
SetTubeProperties: PUBLIC PROC [
tube: Tube, scale, taper, epsilon: REAL, recursively: BOOL ¬ TRUE]
~ {
Set: TubeProc ~ {
tube.scale ¬ scale;
tube.taper ¬ taper;
tube.epsilon ¬ epsilon;
};
IF recursively THEN G3dTube.ApplyToTube[tube, Set] ELSE [] ¬ Set[tube];
};
PropagateEpsilon: PUBLIC PROC [tube: Tube, epsilon: REAL] ~ {
tubeProc: TubeProc ~ {tube.epsilon ¬ epsilon};
G3dTube.ApplyToTube[tube, tubeProc];
};
PropagateScale: PUBLIC PROC [tube: Tube, scale: REAL] ~ {
tubeProc: TubeProc ~ {tube.scale ¬ scale};
G3dTube.ApplyToTube[tube, tubeProc];
};
Tension: PUBLIC PROC [tube: Tube, t: REAL] RETURNS [REAL] ~ {
RETURN[(1.0-t)*tube.tens0+t*tube.tens1];
};
Twist: PUBLIC PROC [tube: Tube, t: REAL] RETURNS [REAL] ~ {
RETURN[(1.0-t)*tube.tw0+t*tube.tw1];
};
Contour Procedures
NContours: PUBLIC PROC [tube: Tube] RETURNS [INTEGER] ~ {
RETURN[IF tube.contours # NIL THEN tube.contours.length ELSE 0];
};
TestContourSequence: PROC [contours: ContourSequence] RETURNS [ContourSequence] ~ {
IF contours = NIL OR contours.length >= contours.maxLength
THEN contours ¬ LengthenContourSequence[contours];
RETURN[contours];
};
LengthenContourSequence: PROC [contours: ContourSequence]
RETURNS [ContourSequence]
~ {
nContours: CARDINAL ¬ IF contours = NIL THEN 0 ELSE contours.length;
new: ContourSequence ¬ NEW[ContourSequenceRep[MAX[5, Real.Round[1.3*nContours]]]];
FOR i: CARDINAL IN[0..nContours) DO new[i] ¬ contours[i]; ENDLOOP;
FOR i: CARDINAL IN[nContours..new.maxLength) DO new[i] ¬ NEW[ContourRep]; ENDLOOP;
new.length ¬ nContours;
RETURN[new];
};
AddContour: PUBLIC PROC [tube: Tube, contour: Contour] ~ {
InsertContour: PROC [contour: Contour, contours: ContourSequence, n: NAT] ~ {
FOR i: INT DECREASING IN [n..contours.length) DO contours[n+1] ¬ contours[n]; ENDLOOP;
contours[n] ¬ contour;
};
new: Contour ¬ NEW[ContourRep ¬ contour­];
tube.contours ¬ TestContourSequence[tube.contours];
SELECT TRUE FROM
tube.contours.length = 0 =>
tube.contours[0] ¬ new;
new.t < tube.contours[0].t =>
InsertContour[new, tube.contours, 0];
new.t > tube.contours[tube.contours.length-1].t =>
tube.contours[tube.contours.length] ¬ new;
ENDCASE => {
n: NAT ¬ 0;
WHILE tube.contours[n].t < new.t DO n ¬ n+1; ENDLOOP;
IF tube.contours[n].t = new.t THEN {
tube.contours[n] ¬ new;
RETURN;
};
InsertContour[new, tube.contours, n];
};
tube.contours.length ¬ tube.contours.length+1;
};
DivideContours: PUBLIC PROC [tube, tube0, tube1: Tube, t: REAL] ~ {
circle: Contour ¬ G2dContour.Circle[tube.circleRes];
NewTs: PROC [tube: Tube, t0, t1: REAL] ~ {
div: REAL ¬ IF t0 = t1 THEN 1.0 ELSE 1.0/(t1-t0);
FOR n: NAT IN [0..tube.contours.length) DO
tube.contours[n].t ¬ (tube.contours[n].t-t0)*div;
ENDLOOP;
};
IF tube.contours # NIL THEN {
contours: ContourSequence ¬ tube.contours;    -- tube may = tube0 or tube1
divideContour: Contour ¬ TContour[contours, circle, t];
tube0.contours ¬ tube1.contours ¬ NIL;
AddContour[tube1, divideContour];
FOR n: NAT IN [0..contours.length) DO
SELECT contours[n].t FROM
<= t => AddContour[tube0, contours[n]];
>= t => AddContour[tube1, contours[n]];
ENDCASE;
ENDLOOP;
AddContour[tube0, divideContour];
NewTs[tube0, 0.0, t];
NewTs[tube1, t, 1.0];
};
};
Alpha: PROC [t, t0, t1: REAL] RETURNS [REAL] ~ {
RETURN[IF t0 = t1 THEN t0 ELSE (t-t0)/(t1-t0)];
};
TContour: PUBLIC PROC [contours: ContourSequence, circle: Contour, t: REAL]
RETURNS [Contour] ~ {
IF contours = NIL
THEN RETURN[circle]
ELSE {
c0, c1: Contour ¬ NIL;
FOR n: NAT IN [0..contours.length) WHILE contours[n].t <= t DO
c0 ¬ contours[n];
ENDLOOP;
FOR n: NAT DECREASING IN [0..contours.length) WHILE contours[n].t >= t DO
c1 ¬ contours[n];
ENDLOOP;
RETURN[SELECT TRUE FROM
c0 = NIL AND c1 = NIL => circle,
c0 = NIL => G2dContour.Interpolate[circle, c1, Alpha[t, 0.0, c1.t]],
c1 = NIL => G2dContour.Interpolate[c0, circle, Alpha[t, c0.t, 1.0]],
ENDCASE => G2dContour.Interpolate[c0, c1, Alpha[t, c0.t, c1.t]]];
};
};
Skin: PUBLIC PROC [tube: Tube] ~ {
InnerSkin: TubeProc ~ {
circle: Contour ¬ G2dContour.Circle[tube.circleRes];
FOR n: INTEGER IN [0..G3dTube.NXSections[tube]) DO
f: XSection ¬ tube.xSections[n];
f.contour ¬ TContour[tube.contours, circle, f.t];
f.normals ¬ SELECT TRUE FROM
f.contour = circle => f.contour.pairs,
ENDCASE => G2dContour.Normals[f.contour];
ENDLOOP;
};
G3dTube.ApplyToTube[tube, InnerSkin];
};
PropagateCircleRes: PUBLIC PROC [tube: Tube, circleRes: INTEGER] ~ {
tubeProc: TubeProc ~ {
IF tube = NIL THEN RETURN;
tube.circleRes ¬ circleRes;
FOR n: NAT IN [0..G3dTube.NXSections[tube]) DO
f: XSection ¬ tube.xSections[n];
IF f.contour # NIL AND f.contour.circle THEN f.contour ¬ circle;
ENDLOOP;
};
circle: Contour ¬ G2dContour.Circle[circleRes];
G3dTube.ApplyToTube[tube, tubeProc];
};
ScreenCircleRes: PUBLIC PROC [tube: Tube, view: Matrix] RETURNS [INTEGER] ~ {
CircleRes: PROC [p: Triple, r: REAL] RETURNS [NAT] ~ {
q: Vector2.VEC ¬ G3dMatrix.TransformD[p, view];
q1: Vector2.VEC ¬ G3dMatrix.TransformD[[p.x+r, p.y, p.z], view];
q2: Vector2.VEC ¬ G3dMatrix.TransformD[[p.x, p.y+r, p.z], view];
l: REAL ¬ MAX[Vector2.Length[Vector2.Sub[q, q1]], Vector2.Length[Vector2.Sub[q, q2]]];
RETURN[3+Real.Round[0.4*l]];
};
res: INTEGER ¬ IF tube.xSections = NIL
THEN 0
ELSE MAX[
CircleRes[tube.p0, tube.xSections[0].scale],
CircleRes[tube.p1, tube.xSections[tube.xSections.length-1].scale]];
IF res MOD 2 = 1 THEN res ¬ res+1;
RETURN[res];
};
Radius Procedures
RoundTipInfo: TYPE ~ RECORD [needTo: BOOL, end: Triple];
GetRoundTipInfo: PROC [tube: Tube, t: REAL] RETURNS [i: RoundTipInfo] ~ {
i ¬ IF t < 0.5
THEN IF tube.prev # NIL
THEN [FALSE, origin]
ELSE [TRUE, tube.p0]
ELSE IF tube.next # NIL OR G3dTube.NBranches[tube] # 0
THEN [FALSE, origin]
ELSE [TRUE, tube.p1];
};
Radius: PUBLIC PROC [tube: Tube, t: REAL, roundTip: BOOL ¬ FALSE] RETURNS [r: REAL] ~ {
r ¬ (1.0-t)*tube.r0+t*tube.r1;
IF tube.scale # 1.0 THEN r ¬ r*tube.scale;
IF tube.taper # 0.0 THEN r ¬ r*(1.0-tube.taper*Length[tube, t]);
IF roundTip THEN {
i: RoundTipInfo ¬ GetRoundTipInfo[tube, t];
IF i.needTo THEN {
p: Triple ¬ G3dSpline.Position[tube.spline, t];
d, dSq, rSq: REAL;
IF (dSq ¬ G3dVector.SquareDistance[i.end, p]) > (rSq ¬ r*r) THEN RETURN;
IF (d ¬ r-RealFns.SqRt[dSq]) > 0.0 THEN r ¬ RealFns.SqRt[rSq-d*d];
};
};
};
SetRadii: PUBLIC PROC [tube: Tube, radius: REAL] ~ {
tubeProc: TubeProc ~ {
tube.r0 ¬ tube.r1 ¬ radius;
IF tube.xSectionsValid THEN ReScale[tube, tube.scale, tube.taper];
};
G3dTube.ApplyToTube[tube, tubeProc];
};
SetDescendantRadii: PUBLIC PROC [tube: Tube, mode: RadiusMode] ~ {
tubeProc: TubeProc ~ {
summer: TubeProc ~ {
radius0: REAL ¬ Radius[tube, 0.0];
sum ¬ sum+(IF mode = square THEN radius0*radius0 ELSE radius0)};
sum: REAL ¬ 0.0;
G3dTube.ApplyToBranches[tube, summer];
IF sum # 0.0 THEN {
radius1: REAL ¬ Radius[tube, 1.0];
factor: REAL~IF mode=square THEN RealFns.SqRt[radius1*radius1/sum] ELSE radius1/sum;
setter: TubeProc ~ {
taper: REAL ~ IF tube.r0 # 0.0 THEN tube.r1/tube.r0 ELSE 1.0;
tube.r0 ¬ factor*tube.r0;
tube.r1 ¬ tube.r0*taper;
IF tube.xSectionsValid THEN ReScale[tube, tube.scale, tube.taper];
};
G3dTube.ApplyToBranches[tube, setter];
};
};
G3dTube.ApplyToTube[tube, tubeProc];
};
Length Procedures
Length: PUBLIC PROC [tube: Tube, t: REAL] RETURNS [REAL] ~ {
IF NOT tube.lengthsValid THEN SetSectionLengths[tube];
RETURN[tube.length0+t*tube.length];
};
SetLengths: PUBLIC PROC [tube: Tube] ~ {
tubeProc: TubeProc ~ {SetSectionLengths[tube]};
G3dTube.ApplyToTube[tube, tubeProc];
};
SetSectionLengths: PUBLIC PROC [tube: Tube] ~ {
IF tube = NIL OR tube.lengthsValid THEN RETURN;
tube.length0 ¬ IF tube.prev = NIL THEN 0.0 ELSE tube.prev.length1;
tube.length ¬ G3dSpline.Length[tube.spline];
tube.length1 ¬ tube.length0+tube.length;
tube.lengthsValid ¬ TRUE;
AdjustXSectionLengths[tube];
};
AdjustXSectionLengths: PUBLIC PROC [tube: Tube] ~ {
IF tube.length # 0.0 AND tube.xSections # NIL AND tube.xSections.length > 0 THEN {
factor: REAL ¬ tube.length/tube.xSections[tube.xSections.length-1].length;
FOR n: NAT IN [0..tube.xSections.length) DO
tube.xSections[n].length ¬ factor*tube.xSections[n].length;
ENDLOOP;
};
};
Plane Procedures
SetPlanes: PUBLIC PROC [tube: Tube] ~ {
tubeProc: TubeProc ~ {SetPlane[tube]};
G3dTube.ApplyToTube[tube, tubeProc];
};
SetPlane: PUBLIC PROC [tube: Tube] ~ {
tube.plane ¬ G3dPlane.FromPointAndNormal[tube.p1, tube.v1, TRUE];
tube.planeValid ¬ TRUE;
};
XSections
Make: PUBLIC PROC [
tube: Tube,
skin: BOOL ¬ FALSE,
roundTip: BOOL ¬ FALSE,
view: Matrix ¬ NIL,
epsilon: REAL ¬ 0.03,
xSectionProc: XSectionProc ¬ NIL]
~ {
tube.epsilon ¬ epsilon;
SetSpline[tube];
SetLengths[tube];
SetPlanes[tube];
MakeXSections[tube, skin, roundTip, view, xSectionProc];
};
MakeXSections: PUBLIC PROC [
tube: Tube,
skin: BOOL ¬ FALSE,
roundTip: BOOL ¬ FALSE,
view: Matrix ¬ NIL,
xSectionProc: XSectionProc ¬ NIL]
~ {
tubeProc: TubeProc ~ {MakeSectionXSections[tube, skin, roundTip, view, xSectionProc]};
G3dTube.ApplyToTube[tube, tubeProc];
};
MakeSectionXSections: PUBLIC PROC [
tube: Tube,
skin: BOOL ¬ FALSE,
roundTip: BOOL ¬ FALSE,
view: Matrix ¬ NIL,
xSectionProc: XSectionProc ¬ NIL]
~ {
tEpsilon: REAL ~ 0.0005;
bez: Bezier ← G3dSpline.BezierFromSpline[tube.spline];
bez: Bezier ← BezierFromHermite[tube.p0, tube.p1, tube.v0, tube.v1];
bez: Bezier ¬ tube.bezier;
maxLen: REAL ¬ IF view = NIL THEN 0.01*G3dSpline.ConvexHullLength[bez] ELSE 0.0;
vv: Triple ¬ G3dVector.Unit[IF G3dVector.Equal[bez.b1, bez.b0, 0.0001]
THEN G3dVector.Sub[bez.b2, bez.b1]
ELSE G3dVector.Sub[bez.b1, bez.b0]];
circle: Contour ¬ G2dContour.Circle[tube.circleRes];
c0: Contour ¬ TContour[tube.contours, circle, 0.0];
c1: Contour ¬ TContour[tube.contours, circle, 1.0];
MakeXSection: PROC [position, velocity: Triple, t: REAL, c: Contour] ~ {
velocity ¬ G3dVector.Unit[velocity];
tube.xSections ¬ TestXSectionSequence[tube.xSections];
IF xSectionProc # NIL
THEN tube.xSections[tube.xSections.length] ¬ xSectionProc[position, velocity, t]
ELSE {
fPrevious: XSection ¬ IF tube.xSections.length = 0
THEN NIL ELSE tube.xSections[tube.xSections.length-1];
f: XSection ¬ tube.xSections[tube.xSections.length];
f.t ¬ t;
f.frame.position ¬ position;
f.length ¬ IF fPrevious = NIL
THEN 0.0
ELSE fPrevious.length+
G3dVector.Distance[f.frame.position, fPrevious.frame.position];
f.frame.triad ¬ Basis[velocity, vv, tube.refVec];
Keep scale non-zero (to avoid subsequent surface normal computation problems):
f.scale ¬ MAX[.000001, Radius[tube, t, roundTip]];
f.twist ¬ tube.tw0+t*(tube.tw1-tube.tw0);
f.matrix ¬ RefMatrix[f.frame.position, f.frame.triad, f.scale, f.twist, f.matrix];
IF skin THEN {
f.contour ¬ c;
f.normals ¬ SELECT TRUE FROM
c.circle => c.pairs,
c.normals # NIL => c.normals,
ENDCASE => G2dContour.Normals[c];
};
};
tube.xSections.length ¬ tube.xSections.length+1;
vv ¬ velocity;
tube.refVec ¬ tube.xSections[tube.xSections.length-1].frame.triad.n;
};
Walker: PROC [bez: Bezier, t0, t1: REAL, c0, c1: Contour] ~ {
IF DividedEnough[bez, t0, t1, c0, c1]
THEN MakeXSection[
bez.b0,
IF G3dVector.Equal[bez.b1, bez.b0, 0.0001] -- degenerate end tangent check
THEN G3dVector.Sub[bez.b2, bez.b1]
ELSE G3dVector.Sub[bez.b1, bez.b0],
t0,
c0]
ELSE {
bez0, bez1: Bezier;
t: REAL ¬ 0.5*(t0+t1);
c: Contour ¬ IF skin THEN TContour[tube.contours, circle, t] ELSE NIL;
[bez0, bez1] ¬ G3dSpline.SplitBezier[bez];
Walker[bez0, t0, t, c0, c];
Walker[bez1, t, t1, c, c1];
};
};
DividedEnough: PROC [bez: Bezier, t0, t1: REAL, c0, c1: Contour] RETURNS [BOOL] ~ {
BezSmallEnough: PROC [bez: Bezier] RETURNS[BOOL] ~ {
RETURN[IF view = NIL
THEN G3dSpline.ConvexHullLength[bez] < maxLen OR G3dSpline.FlatBezier[bez, tube.epsilon]
ELSE Cubic2.Flat[
[G3dMatrix.TransformD[bez.b0, view], G3dMatrix.TransformD[bez.b1, view],
G3dMatrix.TransformD[bez.b2, view], G3dMatrix.TransformD[bez.b3, view]],
100.0*tube.epsilon]];
};
TwistSmallEnough: PROC [t0, t1: REAL] RETURNS[BOOL] ~ {
RETURN[ABS[(t1-t0)*(tube.tw1-tube.tw0)] < 20.0];
};
ContoursCloseEnough: PROC [c0, c1: Contour] RETURNS [ret: BOOL] ~ {
AllCircles: PROC [contours: ContourSequence] RETURNS [BOOL] ~ {
IF contours # NIL THEN FOR n: NAT IN [0..contours.length) DO
IF NOT contours[n].circle THEN RETURN [FALSE];
ENDLOOP;
RETURN[TRUE];
};
IF AllCircles[tube.contours]
THEN {
IF roundTip
THEN {
i0: RoundTipInfo ¬ GetRoundTipInfo[tube, t0];
i1: RoundTipInfo ¬ GetRoundTipInfo[tube, t1];
IF NOT i0.needTo AND NOT i1.needTo
THEN RETURN[TRUE]
ELSE {
Angle: PROC [t: REAL] RETURNS [REAL] ~ {
p: Triple ¬ G3dSpline.Position[tube.spline, t];
d: REAL ¬ r-G3dVector.Distance[p, i0.end];
RETURN[IF d < 0.0 THEN 90.0 ELSE G2dBasic.ArcCosDeg[d/r]];
};
r, ang0, ang1: REAL ¬ IF t1 > 0.5 THEN tube.r1 ELSE tube.r0;
RETURN[r = 0.0 OR ABS[Angle[t0]-Angle[t1]] < 15.0+1500.0*tube.epsilon];
}
}
ELSE RETURN[TRUE];
};
RETURN[G2dContour.Similar[c0, c1] > 1.0-10.0*tube.epsilon];
};
Note the tEpsilon may not be small enough for long, thin cylinders:
IF t1-t0 < tEpsilon THEN RETURN[TRUE];
RETURN[
BezSmallEnough[bez] AND
TwistSmallEnough[t0, t1] AND
ContoursCloseEnough[c0, c1]
];
};
IF tube.xSections # NIL THEN tube.xSections.length ¬ 0;
tube.refVec ¬ SELECT TRUE FROM
tube.prev # NIL => tube.prev.refVec,     -- maintain continuity with parent
tube.refVec # [] => tube.refVec,      -- client-defined xSection orientation
ENDCASE => G3dVector.Unit[
G3dSpline.RefVec[tube.spline, 0.0]];    -- a vector othogonal to base direction
Walker[bez, 0.0, 1.0, c0, c1];
MakeXSection[
bez.b3,
IF G3dVector.Equal[bez.b3, bez.b2, 0.0001] -- test for degenerate end tangent
THEN G3dVector.Sub[bez.b2, bez.b1]
ELSE G3dVector.Sub[bez.b3, bez.b2],
1.0,
c1];
IF view # NIL THEN tube.circleRes ¬
Real.Round[ScreenCircleRes[tube, view]/MAX[0.01, 100*tube.epsilon]];
AdjustXSectionLengths[tube];
tube.xSectionsValid ¬ TRUE;
};
Altering Tube XSections
ReScale: PUBLIC PROC [tube: Tube, scale: REAL ¬ 1.0, taper: REAL ¬ 0.0] ~ {
tubeProc: TubeProc ~ {
scale0: REAL ¬ scale*tube.r0;
scale1: REAL ¬ scale*tube.r1;
tube.scale ¬ scale;
tube.taper ¬ taper;
IF taper # 0.0 THEN {
scale0 ¬ scale0*(1.0-taper*tube.length0);
scale1 ¬ scale1*(1.0-taper*tube.length1);
};
FOR n: NAT IN [0..G3dTube.NXSections[tube]) DO
f: XSection ¬ tube.xSections[n];
f.scale ¬ MAX[0.000001, scale0+f.t*(scale1-scale0)];
f.matrix ¬ RefMatrix[f.frame.position, f.frame.triad, f.scale, f.twist, f.matrix];
ENDLOOP;
};
G3dTube.ApplyToTube[tube, tubeProc];
};
SetCircles: PUBLIC PROC [tube: Tube, circleRes: NAT] ~ {
tubeProc: TubeProc ~ {
tube.circleRes ¬ circleRes;
FOR n: NAT IN [0..G3dTube.NXSections[tube]) DO
f: XSection ¬ tube.xSections[n];
f.contour ¬ circle;
f.normals ¬ circle.normals;
ENDLOOP;
};
circle: Contour ¬ G2dContour.Circle[circleRes];
G3dTube.ApplyToTube[tube, tubeProc];
};
Nearness Operations
NearestTube: PUBLIC PROC [q: Triple, tube: Tube, nearValid: BOOL ¬ TRUE]
RETURNS [nearest: TubePlace]
~ {
tubeProc: TubeProc ~ {
IF NOT nearValid
THEN tube.nearSpline ¬ G3dSpline.NearestPoint[q, tube.spline, 0.0, 1.0, 0.0];
IF tube.nearSpline.distance > minDist THEN RETURN;
minDist ¬ tube.nearSpline.distance;
nearest.tube ¬ tube;
nearest.t ¬ tube.nearSpline.t;
};
minDist: REAL ¬ 10000.0;
G3dTube.ApplyToTube[tube, tubeProc];
};
NearestTube: PUBLIC PROC [q: Triple, tube: Tube, nearValid: BOOL ← TRUE]
RETURNS [nearest: TubePlace]
~ {
minDist: REAL ← Real.LargestNumber;
IF nearValid
THEN {
tubeProc: TubeProc ~ {
IF tube.nearSpline.distance > minDist THEN RETURN;
minDist ← tube.nearSpline.distance;
nearest.tube ← tube;
nearest.t ← tube.nearSpline.t;
};
G3dTube.ApplyToTube[tube, tubeProc];
}
ELSE {
HullDistance: TYPE ~ RECORD [sqDistance: REAL, tube: Tube];
HullDistances: TYPE ~ RECORD [element: SEQUENCE maxLength: NAT OF HullDistance];
tubeProc: TubeProc ~ {
SqDistanceToHull: PROC RETURNS [REAL] ~ {
Need to treat bezier points as a tetrahedron; if q is outside tetrahedron, then return distance from q to nearest face; but if q is inside, would need to return actual distance to curve. So bag it for now.
};
i: NAT ← 0;
sqDistance: REAL ← SqDistanceToHull[tube.bezier, q];
WHILE i < n DO
IF sqDistance < hullDistances[i].sqDistance THEN {
FOR ii: NAT DECREASING IN (i..n] DO
hullDistances[ii] ← hullDistances[ii-1]; -- make space
ENDLOOP;
EXIT;
};
i ← i+1;
ENDLOOP;
hullDistances[i] ← [sqDistance, tube];  -- insert
n ← n+1;
};
n: NAT ← 0;
nTubes: NAT ← G3dTube.NTubes[tube];
hullDistances: REF HullDistances ← NEW[HullDistances[nTubes]];
G3dTube.ApplyToTube[tube, tubeProc];
FOR n: NAT IN [0..nTubes) DO
d: REAL;
IF hullDistances[n].sqDistance > minDist THEN EXIT;
IF (d ← G3dSpline.NearestPoint[q, tube.spline, 0.0, 1.0, 0.0]) < minDist THEN {
minDist ← d;
nearest ← [tube, t];
};
ENDLOOP;
};
};
NearestTube2d: PUBLIC PROC [tube: Tube, screen: Pair, mouseDown: BOOL, view: Matrix]
RETURNS [near: TubePlace]
~ {
nilXSpline: BOOL ¬ FALSE;
t, dist, minDist: REAL ¬ 1000.0;
hasPersp: BOOL ¬ G3dMatrix.HasPerspective[view];
tubeProc: TubeProc ~ {
near2d: NearSpline2d;
IF tube.xSpline = NIL THEN {
tube.xSpline ¬ G3dMatrix.Mul[tube.spline, view, tube.xSpline];
nilXSpline ¬ TRUE;
};
near2d ¬ G3dSpline.NearestPair[screen, tube.xSpline, hasPersp];
IF near2d.distance < minDist THEN {
minDist ¬ near2d.distance;
near.t ¬ SELECT near2d.t FROM < 0.05 => 0.0, > 0.95 => 1.0, ENDCASE => near2d.t;
near.tube ¬ tube;
};
};
IF mouseDown THEN {
tubeProc: TubeProc ~ {tube.xSpline ¬ G3dMatrix.Mul[tube.spline, view, tube.xSpline]};
G3dTube.ApplyToTube[tube, tubeProc];
};
G3dTube.ApplyToTube[tube, tubeProc];
IF nilXSpline THEN {MessageWindow.Append["NilXSpline", TRUE]; MessageWindow.Blink[]};
};
Rotation Operations
Rotate: PUBLIC PROC [tube: Tube, line: Ray, angle: REAL, degrees, tubeAlso: BOOL ¬ TRUE]
~ {
IF tube # NIL THEN {
tubeProc: TubeProc ~ {
tube.p0 ¬ G3dMatrix.Transform[tube.p0, m];
tube.p1 ¬ G3dMatrix.Transform[tube.p1, m];
tube.v0 ¬ G3dMatrix.TransformVec[tube.v0, m];
tube.v1 ¬ G3dMatrix.TransformVec[tube.v1, m];
G3dTube.ApplyToBranches[tube, tubeProc];
};
m: Matrix ¬ G3dMatrix.ObtainMatrix[];
m ¬ G3dMatrix.MakeRotateAbout[line.axis, angle, degrees, line.base, m];
IF tubeAlso THEN [] ¬ tubeProc[tube];
G3dTube.ApplyToBranches[tube, tubeProc];
G3dMatrix.ReleaseMatrix[m];
};
};
XSection Operations
TotalNXSections: PUBLIC PROC [tube: Tube] RETURNS [nXSections: NAT ¬ 0] ~ {
CountNXSections: TubeProc ~ {nXSections ¬ nXSections+G3dTube.NXSections[tube]};
G3dTube.ApplyToTube[tube, CountNXSections];
};
TestXSectionSequence: PROC [xSections: XSectionSequence] RETURNS [XSectionSequence] ~ {
IF xSections = NIL OR xSections.length >= xSections.maxLength
THEN xSections ¬ LengthenXSectionSequence[xSections];
RETURN[xSections];
};
LengthenXSectionSequence: PROC [xSections: XSectionSequence] RETURNS [XSectionSequence] ~ {
nXSections: INTEGER ¬ IF xSections = NIL THEN 0 ELSE xSections.length;
temp: XSectionSequence ¬ NEW[XSectionSequenceRep[MAX[5, Real.Round[1.3*nXSections]]]];
FOR i: NAT IN [0..nXSections) DO temp[i] ¬ xSections[i]; ENDLOOP;
FOR i: NAT IN [nXSections..temp.maxLength) DO temp[i] ¬ NEW[XSectionRep]; ENDLOOP;
temp.length ¬ nXSections;
RETURN[temp];
};
GetXSection: PUBLIC PROC [tube: Tube, t: REAL] RETURNS [XSection] ~ {
nXSections: INTEGER ¬ G3dTube.NXSections[tube];
xSections: XSectionSequence ¬ tube.xSections;
SELECT TRUE FROM
nXSections = 0 =>
RETURN[NIL];
t < xSections[0].t =>
RETURN[xSections[0]];
t > tube.xSections[nXSections-1].t =>
RETURN[xSections[nXSections-1]];
ENDCASE => {
n: NAT ¬ 0;
WHILE xSections[n].t < t DO n ¬ n+1; ENDLOOP;
IF xSections[n].t = t
THEN RETURN[xSections[n]]
ELSE {
xSection: XSection ¬ NEW[XSectionRep];
xSection0: XSection ¬ xSections[n-1];
xSection1: XSection ¬ xSections[n];
alpha: REAL ¬ (t-xSection0.t)/(xSection1.t-xSection0.t);
unitTangent: Triple ¬ G3dVector.Unit[G3dSpline.Velocity[tube.spline, t]];
xSection.t ¬ t;
xSection.frame.position ¬ G3dSpline.Position[tube.spline, t];
xSection.frame.triad ¬
Basis[unitTangent, xSection0.frame.triad.v, xSection0.frame.triad.n];
xSection.scale ¬ xSection0.scale+alpha*(xSection1.scale-xSection0.scale);
xSection.twist ¬ xSection0.twist+alpha*(xSection1.twist-xSection0.twist);
xSection.length ¬ xSection0.length+alpha*(xSection1.length-xSection0.length);
xSection.matrix ¬ RefMatrix[
xSection.frame.position, xSection.frame.triad, xSection.scale, xSection.twist];
xSection.contour ¬ G2dContour.Interpolate[xSection0.contour, xSection1.contour, alpha];
xSection.normals ¬ G2dContour.Normals[xSection.contour];
RETURN[xSection];
};
};
};
PutXSection: PUBLIC PROC [tube: Tube, xSection: XSection] ~ {
InsertXSection: PROC [xSection: XSection, xSections: XSectionSequence, n: NAT] ~ {
FOR i: NAT DECREASING IN [n..xSections.length) DO xSections[n+1] ¬ xSections[n]; ENDLOOP;
xSections[n] ¬ xSection;
};
tube.xSections ¬ TestXSectionSequence[tube.xSections];
SELECT TRUE FROM
tube.xSections.length = 0 =>
tube.xSections[0] ¬ xSection;
xSection.t < tube.xSections[0].t =>
InsertXSection[xSection, tube.xSections, 0];
xSection.t > tube.xSections[tube.xSections.length-1].t =>
tube.xSections[tube.xSections.length] ¬ xSection;
ENDCASE => {
n: NAT ¬ 0;
WHILE tube.xSections[n].t < xSection.t DO n ¬ n+1; ENDLOOP;
IF tube.xSections[n].t = xSection.t THEN {
tube.xSections[n] ¬ xSection;
RETURN;
};
InsertXSection[xSection, tube.xSections, n];
};
tube.xSections.length ¬ tube.xSections.length+1;
};
CopyXSection: PUBLIC PROC [xSection: XSection] RETURNS [XSection] ~ {
copy: XSection ¬ NIL;
IF xSection # NIL THEN {
copy ¬ NEW[XSectionRep];
copy.t ¬ xSection.t;
copy.frame ¬ xSection.frame;
copy.length ¬ xSection.length;
copy.scale ¬ xSection.scale;
copy.twist ¬ xSection.twist;
copy.matrix ¬ G3dMatrix.CopyMatrix[xSection.matrix];
copy.contour ¬ G2dContour.Copy[xSection.contour];
copy.normals ¬ G2dBasic.CopyPairSequence[xSection.normals];
};
RETURN[copy];
};
CopyXSections: PUBLIC PROC [xSections: XSectionSequence] RETURNS [XSectionSequence] ~ {
copy: XSectionSequence ¬ NIL;
IF xSections # NIL THEN {
copy ¬ NEW[XSectionSequenceRep[xSections.length]];
copy.length ¬ xSections.length;
FOR n: NAT IN [0..copy.length) DO copy[n] ¬ CopyXSection[xSections[n]]; ENDLOOP;
};
RETURN[copy];
};
Basis: PUBLIC PROC [v, vv, rv: Triple] RETURNS [Triad] ~ {
UnitTest: PROC [v: Triple] ~ {IF ABS[G3dVector.Square[v]-1.0] > 0.001 THEN ERROR};
v, vv, rv are presumed unit length!
dot: REAL ¬ G3dVector.Dot[v, vv];
b, n: Triple;
vv ← G3dVector.Unit[vv];
v ← G3dVector.Unit[v];
rv ← G3dVector.Unit[rv];
UnitTest[v]; UnitTest[vv]; UnitTest[rv];
IF ABS[dot] > 0.9999 THEN {        -- tangents v, vv colinear
n ¬ IF dot > 0.0 THEN rv ELSE [-rv.x, -rv.y, -rv.z];
n ¬ G3dVector.V90[v, n];        -- ensure orthogonality
}
ELSE {
a: REAL ¬ G2dBasic.ArcCos[dot];
axis: Triple ¬ G3dVector.Cross[vv, v];
n ¬ G3dVector.Unit[G3dVector.RotateAbout[rv, axis, a, FALSE]]; -- the kludge
};
b ¬ G3dVector.Unit[G3dVector.Cross[v, n]];
RETURN[[v, n, b]];
};
RefMatrix: PUBLIC PROC [position: Triple, triad: Triad, scale, twist: REAL, out: Matrix ¬ NIL]
RETURNS [m: Matrix]
~ {
x: Triple ¬ G3dVector.Unit[triad.n];
y: Triple ¬ G3dVector.Unit[triad.b];
z: Triple ¬ G3dVector.Unit[triad.v];
IF ABS[twist] > 0.001 THEN {
out ¬ G3dMatrix.MakeRotate[z, twist,, out];
x ¬ G3dMatrix.TransformVec[x, out];
y ¬ G3dMatrix.TransformVec[y, out];
};
m ¬ G3dMatrix.LocalScale[G3dMatrix.MakeFromTriad[x, y, z, position], scale, out];
};
XYFromXSection: PUBLIC PROC [point: Triple, tube: Tube, t: REAL] RETURNS [Pair] ~ {
ThreeTriples: TYPE ~ RECORD [o, x, y: Triple]; -- origin, x and y vectors
GetOXY: PROC [tube: Tube, t: REAL] RETURNS [tt: ThreeTriples] ~ {
OXYFromXSection: PROC [f: XSection] RETURNS [tt: ThreeTriples] ~ {
m: Matrix ~ f.matrix;
tt.o ¬ f.frame.position;
tt.x ¬ G3dVector.Unit[[m[0][0], m[0][1], m[0][2]]];
tt.y ¬ G3dVector.Unit[[m[1][0], m[1][1], m[1][2]]];
};
nXSections: INTEGER ~ G3dTube.NXSections[tube];
SELECT TRUE FROM
nXSections = 0 =>      RETURN[[origin, origin, origin]];
t < tube.xSections[0].t =>    RETURN[OXYFromXSection[tube.xSections[0]]];
t > tube.xSections[nXSections-1].t => RETURN[OXYFromXSection[tube.xSections[nXSections-1]]];
ENDCASE => {
n: NAT ¬ 0;
WHILE tube.xSections[n].t < t DO n ¬ n+1; ENDLOOP;
IF tube.xSections[n].t = t THEN RETURN[OXYFromXSection[tube.xSections[n]]]
ELSE {
f0: XSection ¬ tube.xSections[n-1];
f1: XSection ¬ tube.xSections[n];
alpha: REAL ~ (t-f0.t)/(f1.t-f0.t);
x0: Triple ~ [f0.matrix[0][0], f0.matrix[0][1], f0.matrix[0][2]];
x1: Triple ~ [f1.matrix[0][0], f1.matrix[0][1], f1.matrix[0][2]];
y0: Triple ~ [f0.matrix[1][0], f0.matrix[1][1], f0.matrix[1][2]];
y1: Triple ~ [f1.matrix[1][0], f1.matrix[1][1], f1.matrix[1][2]];
tt.o ¬ G3dSpline.Position[tube.spline, t];
tt.x ¬ G3dVector.Unit[G3dVector.Combine[x0, 1.0-alpha, x1, alpha]];
tt.y ¬ G3dVector.Unit[G3dVector.Combine[y0, 1.0-alpha, y1, alpha]];
};
};
};
tt: ThreeTriples ¬ GetOXY[tube, t];
q: Triple ¬ G3dVector.Sub[point, tt.o];   -- q now in tt space
RETURN[[G3dVector.Dot[q, tt.x], G3dVector.Dot[q, tt.y]]];
};
END.
..
MakeMatrices: PUBLIC PROC [s: SplineData, scale: REAL] ~ {
sc0: REAL ← s.rad0*scale;
sc1: REAL ← s.rad1*scale;
tw0: REAL ← s.twist0;
tw1: REAL ← s.twist1;
vv: Triple ← s.v0;
mats: MatrixSeq ← s.mats;
s.refVec ← IF s.prev # NIL THEN s.prev.refVec ELSE G3dSpline.RefVec[s.c, 0.0];
mats[0] ←Misc3d.RefMatrix[s.p0, s.refVec, G3dVector.Cross[vv, s.refVec], vv, sc0, tw0, mats[0]];
FOR i: NAT IN[1..s.axialRes) DO
n, b: Triple;
t: REAL ← i/REAL[s.axialRes-1.0];
p: Triple ← G3dSpline.Position[s.c, t];
v: Triple ← G3dSpline.Tangent[s.c, t];
[n, b] ← Misc3d.Basis[v, vv, s.refVec];
mats[i] ← Misc3d.RefMatrix[p, n, b, v, (1.0-t)*sc0+t*sc1, (1.0-t)*tw0+t*tw1, mats[i]];
vv ← v;
s.refVec ← n;
ENDLOOP;
};
GetXSection: PUBLIC PROC [tube: Tube, t: REAL] RETURNS [XSection] ~ {
xSections: XSectionSequence ← tube.xSections;
numXSections: INT ← IF xSections # NIL THEN xSections.maxLength ELSE 0;
SELECT TRUE FROM
numXSections = 0 =>
RETURN[NIL];
xSection.t < xSections[0].t =>
RETURN[tube.xSections[0]];
xSection.t > tube.xSections[numXSections-1].t =>
RETURN[tube.xSections[numXSections-1];
ENDCASE => {
n: NAT ← 0;
WHILE xSections[n].t < xSection.t DO tube.xSections[n] ← xSections[n]; n ← n+1; ENDLOOP;
IF tube.xSections[n].t = t
THEN RETURN[tube.xSections[n]]
ELSE {    -- the following is sloppy:
xSection0: XSection ← tube.xSections[n-1];
xSection1: XSection ← tube.xSections[n];
m0: Matrix ← xSection0.m;
m1: Matrix ← xSection1.m;
alpha: REAL ← (xSection.t-xSection0.t)/(xSection1.t-xSection0.t);
beta: REAL ← 1.0-alpha;
x: Triple ← G3dVector.Combine[
[m0[0][0], m0[0][1], m0[0][2]], alpha, [m1[0][0], m1[0][1], m1[0][2]], beta];
y: Triple ← G3dVector.Combine[
[m0[1][0], m0[1][1], m0[1][2]], alpha, [m1[1][0], m1[1][1], m1[1][2]], beta];
z: Triple ← G3dVector.Combine[
[m0[2][0], m0[2][1], m0[2][2]], alpha, [m1[2][0], m1[2][1], m1[2][2]], beta];
v: Triple ← G3dVector.SameMag[z, G3dSpline.Velocity[tube.spline, t]];
n: Triple ← G3dVector.SameMag[x, G3dVector.Cross[v, y]];
b: Triple ← G3dVector.SameMag[y, G3dVector.Cross[v, n]];
p: Triple ← G3dSpline.Position[tube.spline, t];
xSection: XSection ← NEW[XSectionRep];
xSection.t ← t;
xSection.c ← Contour.Interpolate[xSection0.c, xSection1.c, alpha];
xSection.m ← RefMatrix[p, n, b, v, 1.0, 0.0];
RETURN[xSection];
};
};
};