<> <> <> DIRECTORY Contours, Cubic2, Matrix3d, Real, RealFns, Spline3d, TubeContour, Vector3d, TubeDefs, TubeGeometry, TubeMisc; TubeGeometryImpl: CEDAR PROGRAM IMPORTS Contours, Cubic2, Matrix3d, Real, RealFns, Spline3d, TubeContour, TubeMisc, Vector3d EXPORTS TubeGeometry ~ BEGIN OPEN TubeDefs; <> SetCoeffs: PUBLIC PROC [tube: Tube] ~ { tubeProc: TubeProc ~ {SetSectionCoeffs[tube]}; TubeMisc.ApplyToTube[tube, tubeProc]; }; SetSectionCoeffs: PUBLIC PROC [tube: Tube] ~ { tube.coeffs _ Spline3d.CoeffsFromHermite[ [tube.p0, tube.p1, Vector3d.Mul[tube.v0, tube.tens0], Vector3d.Mul[tube.v1, tube.tens1]], tube.coeffs]; }; <> SetLengths: PUBLIC PROC [tube: Tube] ~ { tubeProc: TubeProc ~ {SetSectionLengths[tube]}; TubeMisc.ApplyToTube[tube, tubeProc]; }; SetSectionLengths: PUBLIC PROC [tube: Tube] ~ { IF tube = NIL THEN RETURN; IF tube.prev # NIL THEN tube.length0 _ tube.prev.length1; tube.length1 _ tube.length0+Spline3d.Length[tube.coeffs]; }; <> Make: PUBLIC PROC [ tube: Tube, epsilon: REAL _ 0.03, scale: REAL _ 1.0, taper: REAL _ 0.0, skin: BOOL _ FALSE, roundTip: BOOL _ FALSE, view: Matrix _ NIL, frameProc: FrameProc _ NIL] ~ { SetCoeffs[tube]; SetLengths[tube]; MakeFrames[tube, epsilon, scale, taper, skin, roundTip, view, frameProc]; }; MakeFrames: PUBLIC PROC [ tube: Tube, epsilon: REAL _ 0.03, scale: REAL _ 1.0, taper: REAL _ 0.0, skin: BOOL _ FALSE, roundTip: BOOL _ FALSE, view: Matrix _ NIL, frameProc: FrameProc _ NIL] ~ { tubeProc: TubeProc ~ { scale0, scale1: REAL _ scale; IF tube = NIL THEN RETURN; IF taper # 0.0 THEN { scale0 _ scale*(1.0-taper*tube.length0); scale1 _ scale*(1.0-taper*tube.length1); }; MakeSectionFrames[tube, epsilon, scale0, scale1, skin, roundTip, view, frameProc]; }; TubeMisc.ApplyToTube[tube, tubeProc]; }; MakeSectionFrames: PUBLIC PROC [ tube: Tube, epsilon: REAL _ 0.03, scale0, scale1: REAL _ 1.0, skin: BOOL _ FALSE, roundTip: BOOL _ FALSE, view: Matrix _ NIL, frameProc: FrameProc _ NIL] ~ { bez: Bezier _ Spline3d.BezierFromCoeffs[tube.coeffs]; maxLen: REAL _ IF view = NIL THEN 0.01*Spline3d.ConvexHullLength[bez] ELSE 0.0; vv: Triple _ IF Vector3d.Equal[bez.b1, bez.b0, 0.0001] THEN Vector3d.Sub[bez.b2, bez.b1] ELSE Vector3d.Sub[bez.b1, bez.b0]; circle: Contour _ Contours.Circle[tube.circleRes]; c0: Contour _ TubeContour.TContour[tube.contours, circle, 0.0]; c1: Contour _ TubeContour.TContour[tube.contours, circle, 1.0]; MakeFrame: PROC [position, velocity: Triple, t: REAL, c: Contour] ~ { tube.frames _ TestFrameSequence[tube.frames]; IF frameProc # NIL THEN tube.frames[tube.frames.length] _ frameProc[position, velocity, t] ELSE { f: Frame _ tube.frames[tube.frames.length]; f.t _ t; f.position _ position; f.triad _ Basis[velocity, vv, tube.refVec]; <> f.scale _ MAX[.000001, (scale0+t*(scale1-scale0))*TubeMisc.Radius[tube, t, roundTip]]; f.twist _ tube.tw0+t*(tube.tw1-tube.tw0); f.matrix _ RefMatrix[f.position, f.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 => Contours.Normals[c]; }; }; tube.frames.length _ tube.frames.length+1; vv _ velocity; tube.refVec _ tube.frames[tube.frames.length-1].triad.n; }; Walker: PROC [bez: Bezier, t0, t1: REAL, c0, c1: Contour] ~ { IF DividedEnough[bez, t0, t1, c0, c1] THEN MakeFrame[ bez.b0, IF Vector3d.Equal[bez.b1, bez.b0, 0.0001] -- degenerate end tangent check THEN Vector3d.Sub[bez.b2, bez.b1] ELSE Vector3d.Sub[bez.b1, bez.b0], t0, c0] ELSE { left, rite: Bezier; t: REAL _ 0.5*(t0+t1); c: Contour _ IF skin THEN TubeContour.TContour[tube.contours, circle, t] ELSE NIL; [left, rite] _ Spline3d.SplitBezier[bez]; Walker[left, t0, t, c0, c]; Walker[rite, 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 Spline3d.ConvexHullLength[bez] < maxLen OR Spline3d.FlatBezier[bez, epsilon] ELSE Cubic2.Flat[ [Matrix3d.TransformD[bez.b0, view], Matrix3d.TransformD[bez.b1, view], Matrix3d.TransformD[bez.b2, view], Matrix3d.TransformD[bez.b3, view]], 100.0*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 { r0: REAL ~ TubeMisc.Radius[tube, t0, TRUE]; r1: REAL ~ TubeMisc.Radius[tube, t1, TRUE]; IF r0 = 0.0 OR r1 = 0.0 THEN RETURN[FALSE]; RETURN[1.0-MIN[r0, r1]/MAX[r0, r1] < 100.0*epsilon]; } ELSE RETURN[TRUE]; }; RETURN[Contours.Similar[c0, c1] > 1.0-10.0*epsilon]; }; IF t1-t0 < 0.005 THEN RETURN[TRUE]; RETURN[ BezSmallEnough[bez] AND TwistSmallEnough[t0, t1] AND ContoursCloseEnough[c0, c1] ]; }; IF tube # NIL AND tube.frames # NIL THEN tube.frames.length _ 0; tube.refVec _ IF tube.prev # NIL THEN tube.prev.refVec ELSE Spline3d.RefVec[tube.coeffs, 0.0]; Walker[bez, 0.0, 1.0, c0, c1]; MakeFrame[ bez.b3, IF Vector3d.Equal[bez.b3, bez.b2, 0.0001] -- test for degenerate end tangent THEN Vector3d.Sub[bez.b2, bez.b1] ELSE Vector3d.Sub[bez.b3, bez.b2], 1.0, c1]; IF view # NIL THEN tube.circleRes _ Real.RoundI[TubeContour.ScreenCircleRes[tube, view]/MAX[0.01, 100*epsilon]]; }; <> 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; IF tube = NIL THEN RETURN; IF taper # 0.0 THEN { scale0 _ scale0*(1.0-taper*tube.length0); scale1 _ scale1*(1.0-taper*tube.length1); }; FOR n: NAT IN [0..NFrames[tube]) DO f: Frame _ tube.frames[n]; f.scale _ MAX[0.000001, scale0+f.t*(scale1-scale0)]; f.matrix _ RefMatrix[f.position, f.triad, f.scale, f.twist, f.matrix]; ENDLOOP; }; TubeMisc.ApplyToTube[tube, tubeProc]; }; SetCircles: PUBLIC PROC [tube: Tube, circleRes: NAT] ~ { tubeProc: TubeProc ~ { tube.circleRes _ circleRes; FOR n: NAT IN [0..NFrames[tube]) DO f: Frame _ tube.frames[n]; f.contour _ circle; f.normals _ circle.normals; ENDLOOP; }; circle: Contour _ Contours.Circle[circleRes]; TubeMisc.ApplyToTube[tube, tubeProc]; }; <> NFrames: PUBLIC PROC [tube: Tube] RETURNS [INTEGER] ~ { RETURN[IF tube.frames # NIL THEN tube.frames.length ELSE 0]; }; TestFrameSequence: PROC [frames: FrameSequence] RETURNS [FrameSequence] ~ { IF frames = NIL OR frames.length >= frames.maxLength THEN frames _ LengthenFrameSequence[frames]; RETURN[frames]; }; LengthenFrameSequence: PROC [frames: FrameSequence] RETURNS [FrameSequence] ~ { nFrames: INTEGER _ IF frames = NIL THEN 0 ELSE frames.length; temp: FrameSequence _ NEW[FrameSequenceRep[MAX[5, Real.RoundI[1.3*nFrames]]]]; FOR i: NAT IN [0..nFrames) DO temp[i] _ frames[i]; ENDLOOP; FOR i: NAT IN [nFrames..temp.maxLength) DO temp[i] _ NEW[FrameRep]; ENDLOOP; temp.length _ nFrames; RETURN[temp]; }; GetFrame: PUBLIC PROC [tube: Tube, t: REAL] RETURNS [Frame] ~ { nFrames: INTEGER _ NFrames[tube]; frames: FrameSequence _ tube.frames; SELECT TRUE FROM nFrames = 0 => RETURN[NIL]; t < frames[0].t => RETURN[frames[0]]; t > tube.frames[nFrames-1].t => RETURN[frames[nFrames-1]]; ENDCASE => { n: NAT _ 0; WHILE frames[n].t < t DO n _ n+1; ENDLOOP; IF frames[n].t = t THEN RETURN[frames[n]] ELSE { frame: Frame _ NEW[FrameRep]; frame0: Frame _ frames[n-1]; frame1: Frame _ frames[n]; alpha: REAL _ (frame.t-frame0.t)/(frame1.t-frame0.t); frame.t _ t; frame.position _ Spline3d.Position[tube.coeffs, t]; frame.triad _ Basis[Spline3d.Velocity[tube.coeffs, t], frame0.triad.v, frame0.triad.n]; frame.scale _ frame0.scale+alpha*(frame1.scale-frame0.scale); frame.twist _ frame0.twist+alpha*(frame1.twist-frame0.twist); frame.matrix _ RefMatrix[frame.position, frame.triad, frame.scale, frame.twist]; frame.contour _ Contours.Interpolate[frame0.contour, frame1.contour, alpha]; frame.normals _ Contours.Normals[frame.contour]; RETURN[frame]; }; }; }; PutFrame: PUBLIC PROC [tube: Tube, frame: Frame] ~ { InsertFrame: PROC [frame: Frame, frames: FrameSequence, n: NAT] ~ { FOR i: NAT DECREASING IN [n..frames.length) DO frames[n+1] _ frames[n]; ENDLOOP; frames[n] _ frame; }; tube.frames _ TestFrameSequence[tube.frames]; SELECT TRUE FROM tube.frames.length = 0 => tube.frames[0] _ frame; frame.t < tube.frames[0].t => InsertFrame[frame, tube.frames, 0]; frame.t > tube.frames[tube.frames.length-1].t => tube.frames[tube.frames.length] _ frame; ENDCASE => { n: NAT _ 0; WHILE tube.frames[n].t < frame.t DO n _ n+1; ENDLOOP; IF tube.frames[n].t = frame.t THEN { tube.frames[n] _ frame; RETURN; }; InsertFrame[frame, tube.frames, n]; }; tube.frames.length _ tube.frames.length+1; }; CopyFrame: PUBLIC PROC [frame: Frame] RETURNS [Frame] ~ { copy: Frame _ NIL; IF frame # NIL THEN { copy _ NEW[FrameRep]; copy.t _ frame.t; copy.position _ frame.position; copy.triad _ frame.triad; copy.scale _ frame.scale; copy.matrix _ Matrix3d.CopyMatrix[frame.matrix]; copy.contour _ Contours.Copy[frame.contour]; }; RETURN[copy]; }; CopyFrames: PUBLIC PROC [frames: FrameSequence] RETURNS [FrameSequence] ~ { copy: FrameSequence _ NIL; IF frames # NIL THEN { copy _ NEW[FrameSequenceRep[frames.length]]; copy.length _ frames.length; FOR n: NAT IN [0..copy.length) DO copy[n] _ CopyFrame[frames[n]]; ENDLOOP; }; RETURN[copy]; }; <> Basis: PUBLIC PROC [v, vv, rv: Triple] RETURNS [Triad] ~ { dot: REAL; b, n: Triple; vv _ Vector3d.Normalize[vv]; v _ Vector3d.Normalize[v]; rv _ Vector3d.Normalize[rv]; dot _ Vector3d.Dot[v, vv]; 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 _ Vector3d.V90[v, n]; -- ensure orthogonality } ELSE { a: REAL _ RealFns.ArcTanDeg[RealFns.SqRt[1.0-dot*dot], dot]; -- == Acos[] axis: Triple _ Vector3d.Cross[v, vv]; n _ Vector3d.Normalize[Vector3d.RotateAbout[rv, axis, a]]; -- incremental kludge }; b _ Vector3d.Normalize[Vector3d.Cross[v, n]]; RETURN[[v, n, b]]; }; RefMatrix: PUBLIC PROC [position: Triple, triad: Triad, scale, twist: REAL, out: Matrix _ NIL] RETURNS [Matrix] ~ { IF ABS[twist] > 0.001 THEN { out _ Matrix3d.MakePureRotate[triad.v, twist, , out]; triad.n _ Matrix3d.TransformVec[triad.n, out]; triad.b _ Matrix3d.TransformVec[triad.b, out]; }; RETURN[ Matrix3d.LocalScale[Matrix3d.MakeFromTriad[triad.n, triad.b, triad.v, position], scale, out]]; }; XYFromFrame: 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] ~ { OXYFromFrame: PROC [f: Frame] RETURNS [tt: ThreeTriples] ~ { m: Matrix ~ f.matrix; tt.o _ f.position; tt.x _ Vector3d.Normalize[[m[0][0], m[0][1], m[0][2]]]; tt.y _ Vector3d.Normalize[[m[1][0], m[1][1], m[1][2]]]; }; nFrames: INTEGER ~ NFrames[tube]; SELECT TRUE FROM nFrames = 0 => RETURN[[origin, origin, origin]]; t < tube.frames[0].t => RETURN[OXYFromFrame[tube.frames[0]]]; t > tube.frames[nFrames-1].t => RETURN[OXYFromFrame[tube.frames[nFrames-1]]]; ENDCASE => { n: NAT _ 0; WHILE tube.frames[n].t < t DO n _ n+1; ENDLOOP; IF tube.frames[n].t = t THEN RETURN[OXYFromFrame[tube.frames[n]]] ELSE { f0: Frame _ tube.frames[n-1]; f1: Frame _ tube.frames[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 _ Spline3d.Position[tube.coeffs, t]; tt.x _ Vector3d.Normalize[Vector3d.Combine[x0, 1.0-alpha, x1, alpha]]; tt.y _ Vector3d.Normalize[Vector3d.Combine[y0, 1.0-alpha, y1, alpha]]; }; }; }; tt: ThreeTriples _ GetOXY[tube, t]; q: Triple _ Vector3d.Sub[point, tt.o]; -- q now in tt space RETURN[[Vector3d.Dot[q, tt.x], Vector3d.Dot[q, tt.y]]]; }; END. .. <> <> <> <> <> <> <> <> <> <> <> <> <> <> <<[n, b] _ Misc3d.Basis[v, vv, s.refVec];>> <> <> <> <> <<};>> <<>> <> <> <> <<>> <> <>> <> <>> <> < tube.frames[numFrames-1].t =>>> <> < {>> <> <> <> <> <> <> <> <> <> <> <> <> <<[m0[0][0], m0[0][1], m0[0][2]], alpha, [m1[0][0], m1[0][1], m1[0][2]], beta];>> <> <<[m0[1][0], m0[1][1], m0[1][2]], alpha, [m1[1][0], m1[1][1], m1[1][2]], beta];>> <> <<[m0[2][0], m0[2][1], m0[2][2]], alpha, [m1[2][0], m1[2][1], m1[2][2]], beta];>> <> <> <> <> <> <> <> <> <> <<};>> <<};>> <<};>> <<>>