<> <> <> DIRECTORY Complex, Cubic, DynFit, LSPiece, Seq, Vector, Nodes, RealFns; NodeImpl: CEDAR PROGRAM IMPORTS Complex, Cubic, DynFit, LSPiece, Vector, RealFns = BEGIN DynNodes: PUBLIC PROC [samples: Seq.ComplexSequence, closed: BOOLEAN, penalty: REAL] RETURNS [nodes: Seq.NatSequence] = { <> endpoints: DynFit.Segments; [segments: endpoints] _ DynFit.FitSegments[samples, penalty]; nodes _ NEW[Seq.NatSequenceRec[endpoints.length]]; nodes[0] _ 0; FOR i:NAT IN [0..endpoints.length) DO nodes[i] _ endpoints[i]; ENDLOOP; }; CubicTangents: PUBLIC PROC [samples: Seq.ComplexSequence, closed: BOOLEAN, err: REAL, maxit: INT, nodes: Seq.NatSequence _ NIL] RETURNS [tangents: Seq.ComplexSequence] = { <> n: NAT _ samples.length; t: Seq.RealSequence _ NEW[Seq.RealSequenceRec[n]]; bezier: Cubic.Bezier; prev: PROC [index: NAT] RETURNS [NAT] = { IF closed THEN RETURN[IF index=0 THEN nodes.length-2 ELSE index-1] ELSE RETURN[index]; }; next: PROC [index: NAT] RETURNS [NAT] = { IF closed THEN RETURN[IF index=nodes.length-1 THEN 1 ELSE index+1] ELSE RETURN[index]; }; pointsBetween: PROC [from, to: NAT] RETURNS [NAT] = { IF closed THEN RETURN[IF to>from THEN to-from+1 ELSE samples.length-from+to+1] ELSE RETURN[to-from]; }; FOR i: NAT IN [0..nodes.length) DO p: NAT _ prev[i]; q: NAT _ nodes[i]; r: NAT _ next[i]; e: REAL; nPoints: NAT _ pointsBetween[p,r]; IF nPoints < minPoints THEN tangents[q] _ CheapTangent[samples[p], samples[q], samples[r], maxSin] ELSE { [b:bezier, maxDev: e] _ LSPiece.FitPiece[z: samples, t: t, from: p, thru: r, eps: .0001, maxd: err/4, maxit: maxit, initFree: freeEnds, finalFree: freeEnds]; tangents[q] _IF e>err THEN [0,0] ELSE TangentVector[bezier,t[q]]; }; ENDLOOP; }; ICubicTangents: PUBLIC PROC [return: PROC[tangent: Complex.Vec, cubic: Cubic.Bezier], samples: Seq.ComplexSequence, closed: BOOLEAN, err: REAL, maxit: INT, nodes: Seq.NatSequence _ NIL] = { <> n: NAT _ samples.length; t: Seq.RealSequence _ NEW[Seq.RealSequenceRec[n]]; lastTangent: Complex.Vec; bezier: Cubic.Bezier; prev: PROC [index: NAT] RETURNS [NAT] = { IF closed THEN RETURN[IF index=0 THEN nodes.length-2 ELSE index-1] ELSE RETURN[index]; }; next: PROC [index: NAT] RETURNS [NAT] = { IF closed THEN RETURN[IF index=nodes.length-1 THEN 1 ELSE index+1] ELSE RETURN[index]; }; pointsBetween: PROC [from, to: NAT] RETURNS [NAT] = { IF closed THEN RETURN[IF to>from THEN to-from+1 ELSE samples.length-from+to+1] ELSE RETURN[to-from]; }; FOR i: NAT IN [0..nodes.length) DO p: NAT _ prev[i]; q: NAT _ nodes[i]; r: NAT _ next[i]; nPoints: NAT _ pointsBetween[p,r]; e: REAL; IF nPoints < minPoints THEN lastTangent _ CheapTangent[samples[p], samples[q], samples[r], maxSin] ELSE { [b:bezier, maxDev: e] _ LSPiece.FitPiece[z: samples, t: t, from: p, thru: r, eps: .0001, maxd: err/4, maxit: maxit, initFree: freeEnds, finalFree: freeEnds]; lastTangent _IF e>err THEN [0,0] ELSE TangentVector[bezier,t[q]]; return[lastTangent,bezier]; }; ENDLOOP; }; freeEnds: BOOLEAN _ TRUE; minPoints: NAT _ 4; maxSin: REAL _ 0.5; TangentVector: PROC [b: Cubic.Bezier, t: REAL] RETURNS [tang: Complex.Vec] = { c: Cubic.Coeffs _ Cubic.BezierToCoeffs[b]; tang.x _ (3*c.c3.x*t+2*c.c2.x)*t + c.c1.x; tang.y _ (3*c.c3.y*t+2*c.c2.y)*t + c.c1.y; }; QuickTangents: PUBLIC PROC [samples: Seq.ComplexSequence, closed: BOOLEAN, maxAngle: REAL, nodes: Seq.NatSequence _ NIL] RETURNS [tangents: Seq.ComplexSequence] = { <> tangents _ DiffTangents[samples, closed, maxAngle, nodes, CheapTangent]; }; SquareTangents: PUBLIC PROC [samples: Seq.ComplexSequence, closed: BOOLEAN, maxAngle: REAL, nodes: Seq.NatSequence _ NIL] RETURNS [tangents: Seq.ComplexSequence] = { <> tangents _ DiffTangents[samples, closed, maxAngle, nodes, SquareTangent]; }; TangentProc: TYPE = PROC[a, b, c: Complex.Vec, maxSin: REAL] RETURNS [t: Complex.Vec]; DiffTangents: PROC [samples: Seq.ComplexSequence, closed: BOOLEAN, maxAngle: REAL, nodes: Seq.NatSequence, tangentProc: TangentProc] RETURNS [tangents: Seq.ComplexSequence]= { maxSin: REAL _ RealFns.SinDeg[maxAngle]; initTan: Complex.Vec _ IF samples[nodes[0]] = samples[nodes[nodes.length-1]] THEN CheapTangent[samples[nodes[nodes.length-2]], samples[nodes[0]], samples[nodes[1]], maxSin] ELSE [0,0]; tangents[0] _ initTan; FOR i: NAT IN [1..nodes.length-1) DO tangents[0] _ CheapTangent[samples[nodes[i-1]], samples[nodes[i]], samples[nodes[i+1]], maxSin]; ENDLOOP; }; CheapTangent: PROC [a, b, c: Complex.Vec, maxSin: REAL] RETURNS [t: Complex.Vec] = { d1: Complex.Vec _ Complex.Sub[b,a]; d2: Complex.Vec _ Complex.Sub[c,b]; absd1d2: REAL _ RealFns.SqRt[Complex.SqrAbs[d1]*Complex.SqrAbs[d2]]; IF Vector.Dot[d1, d2] < 0 OR ABS[Vector.Cross[d1, d2]] >= maxSin * absd1d2 THEN t _ [0,0] ELSE t _ Complex.Sub[c,a]; }; SquareTangent: PROC [a, b, c: Complex.Vec, maxSin: REAL] RETURNS [t: Complex.Vec] = { d1: Complex.Vec _ Complex.Sub[b,a]; d2: Complex.Vec _ Complex.Sub[c,b]; absd1d2: REAL _ RealFns.SqRt[Complex.SqrAbs[d1]*Complex.SqrAbs[d2]]; IF Vector.Dot[d1, d2] < 0 OR ABS[Vector.Cross[d1, d2]] >= maxSin * absd1d2 THEN t _ [0,0] ELSE t _ Complex.Add[Vector.Mul[d1, Complex.Abs[d1]], Vector.Mul[d2, Complex.Abs[d2]]]; }; END. Michael Plass, August 30, 1982 2:06 pm. Added CubicTangents. < . Sets tangents at nodes by locally fitting 2*range + 1 samples>> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <<};>> <<};>> <> <> <> <> <> <> <> <> <<};>> <<};>> <<[nodes:node] _ Curve.CurrentNodes[Curve.defaultHandle];>> <<[cusps: cusp] _ Curve.CurrentCusps[Curve.defaultHandle];>> <> <> <> <> <> <> <> <<[b:bezier] _ LSPiece.FitPiece[z: z, t: t, from: p, thru: q,>> <> <> <<[b:bezier] _ LSPiece.FitPiece[z: z, t: t, from: q, thru: r,>> <> <> <> <> <<}>> <> <> <> <<[b:bezier] _ LSPiece.FitPiece[z: z, t: t, from: p, thru: r,>> <> <> <<};>> <> <> <<};>> <<>>