LSCurveImpl.mesa
Copyright Ó 1985, 1992 by Xerox Corporation. All rights reserved.
Michael Plass and Maureen Stone Oct-81
Last edited by Michael Plass November 29, 1982 9:52 am
Doug Wyatt, September 5, 1985 1:34:22 pm PDT
DIRECTORY
Complex USING [Abs, Add, Arg, Sub, VEC],
Cubic USING [Bezier, BezierToCoeffs, Coeffs],
LinearSystem USING [ColumnN, MatrixN, MatrixSeq, SolveN, VecSeq],
LSCurve USING [Handle, Patch, PatchSequence, PatchSequenceRec, SplineBasis, SplineBasisRec, SplineFunction, StateRec],
PiecewiseCubic USING [Combine, EnumerateCommonPieces, Eval, EvalAll, Handle, Piece, PieceProc, Zero],
Real USING [RealException],
Seq USING [BooleanSequence, BooleanSequenceRec, ComplexSequence, RealSequence, RealSequenceRec],
Vector USING [Div, Dot];
LSCurveImpl: CEDAR PROGRAM
IMPORTS Complex, Cubic, LinearSystem, PiecewiseCubic, Real, Vector
EXPORTS LSCurve
= BEGIN OPEN Seq;
RealProc: TYPE = PROCEDURE[i:NAT] RETURNS [REAL];
PointNumber: TYPE = NAT;
Patch: TYPE = LSCurve.Patch;
PatchSequence: TYPE = LSCurve.PatchSequence;
PatchSequenceRec: TYPE = LSCurve.PatchSequenceRec;
SplineBasis: TYPE = LSCurve.SplineBasis;
SplineBasisRec: TYPE = LSCurve.SplineBasisRec;
SplineFunction: TYPE = LSCurve.SplineFunction;
Handle: TYPE = LSCurve.Handle;
StateRec: TYPE = LSCurve.StateRec;
Create:
PUBLIC
PROCEDURE [sa: Seq.ComplexSequence]
RETURNS [h: Handle] =
BEGIN
h ¬ SamplesToHandle[sa];
h.weight ¬ NEW[RealSequenceRec[h.z.length]];
FOR i:NAT IN [0..h.z.length) DO h.weight[i] ¬ 1.0 ENDLOOP;
h.t ¬ NEW[RealSequenceRec[h.z.length]];
h.oldt ¬ NEW[RealSequenceRec[h.z.length]];
h.oldoldt ¬ NEW[RealSequenceRec[h.z.length]];
END;
XYat:
PUBLIC
PROCEDURE [h: Handle, t:
REAL]
RETURNS [f: Complex.
VEC] =
BEGIN
patchNumber: NAT ¬ 0;
ValidateCache[h];
f.x ¬ PiecewiseCubic.Eval[h.patchCache.x,t];
f.y ¬ PiecewiseCubic.Eval[h.patchCache.y,t];
END;
SamplesToHandle:
PROCEDURE [s: Seq.ComplexSequence]
RETURNS [h: Handle] =
{h ¬ NEW[StateRec]; BEGIN OPEN h;
z ¬ s;
l ¬ 0;
n ¬ s.length-1;
END};
CountTrues:
PROCEDURE [b: BooleanSequence]
RETURNS [count:
NAT] =
BEGIN
count ¬ 0;
FOR i:NAT IN [0..b.length) DO IF b[i] THEN count ¬ count + 1 ENDLOOP;
END;
EvalBasis:
PUBLIC
PROCEDURE [h: Handle, i:
NAT, t:
REAL]
RETURNS [z: Complex.VEC] = INLINE
{OPEN h.splineBasis[i]; RETURN [[PiecewiseCubic.Eval[x,t],PiecewiseCubic.Eval[y,t]]]};
SolveCurve:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h,LinearSystem;
nFree: NAT ¬ CountTrues[free];
b: ColumnN ¬ NEW[VecSeq[nFree]];
c: ColumnN ¬ NIL;
A: MatrixN ¬ NEW[MatrixSeq[nFree]];
FOR i:
NAT
IN [0..nFree)
DO
A[i] ¬ NEW[VecSeq[nFree]];
FOR j: NAT IN [0..nFree) DO A[i][j] ¬ 0 ENDLOOP;
b[i] ¬ 0;
ENDLOOP;
FOR k:
NAT
IN [l..n]
DO
tk:REAL = t[k];
row: NAT ¬ 0;
FOR i:
NAT
IN [0..a.length)
DO
IF free[i]
THEN
BEGIN
gi:Complex.VEC ¬ EvalBasis[h,i,tk];
IF gi#[0,0]
THEN
BEGIN
column: NAT ¬ 0;
FOR j:
NAT
IN [0..a.length)
DO
gj: Complex.VEC ¬ EvalBasis[h,j,tk];
IF free[j]
THEN
A[row][column] ¬ A[row][column] + weight[k]*Vector.Dot[gi,gj]
ELSE
b[row] ¬ b[row] - weight[k]*a[j]*Vector.Dot[gi,gj];
IF free[j] THEN column ¬ column + 1;
ENDLOOP;
END;
b[row] ¬ b[row] + weight[k]*Vector.Dot[z[k],gi];
row ¬ row + 1;
END ENDLOOP;
ENDLOOP;
c ¬ SolveN[A,b,nFree!
Real.RealException => CHECKED {CONTINUE}];
IF c#
NIL
THEN
BEGIN
row: NAT ¬ 0;
FOR i:
NAT
IN [0..a.length)
DO
IF free[i] THEN a[i] ¬ c[row];
IF free[i] THEN row ¬ row + 1;
ENDLOOP;
patchCacheValid ¬ FALSE;
END;
END;
ArcLengthInitialTValues:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
arcLen: REAL ¬ 0;
t[0] ¬ 0;
FOR i:
NAT
IN [1..z.length)
DO
arcLen ¬ arcLen + Complex.Abs[Complex.Sub[z[i],z[i-1]]];
t[i] ¬ arcLen;
ENDLOOP;
FOR i:
NAT
IN [0..z.length)
DO
t[i] ¬ t[i]/arcLen
ENDLOOP;
END;
UnitInitialTValues:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
IF z.length < 2 THEN RETURN;
FOR i:
NAT
IN [0..z.length)
DO
t[i] ¬ i*1.0/(z.length-1);
ENDLOOP;
END;
AngleInitialTValues:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
twopi: REAL = 2*3.1415926;
m: Complex.VEC ¬ [0,0];
FOR i:
NAT
IN [0..z.length)
DO
m ¬ Complex.Add[m,z[i]];
ENDLOOP;
m ¬ Vector.Div[m,z.length];
FOR i:
NAT
IN [0..z.length)
DO
t[i] ¬ Complex.Arg[Complex.Sub[z[i],m]]/twopi;
WHILE t[i] < 0 DO t[i] ¬ t[i] + 1 ENDLOOP;
WHILE t[i] > 1 DO t[i] ¬ t[i] - 1 ENDLOOP;
ENDLOOP;
END;
AdjustTValues:
PUBLIC
PROCEDURE [h: Handle]
RETURNS [maxChange:
REAL] =
BEGIN OPEN h;
delta: REAL;
{newt: RealSequence ¬ oldoldt; oldoldt ¬ oldt; oldt ¬ t; t ¬ newt};
ValidateCache[h];
maxChange ¬ 0;
FOR i:
NAT
IN [l..n]
DO
delta ¬ Adjustment[z[i],patchCache,oldt[i]];
IF ABS[delta] > maxChange THEN maxChange ¬ ABS[delta];
t[i] ¬ oldt[i]+delta;
IF closedCurve
THEN
BEGIN
WHILE t[i] < 0 DO t[i] ¬ t[i] + 1 ENDLOOP;
WHILE t[i] > 1 DO t[i] ¬ t[i] - 1 ENDLOOP;
END
ENDLOOP;
IF initialEndFree
THEN
BEGIN
t0: REAL ¬ t[l];
FOR i:
NAT
IN [l..n]
DO
t[i] ¬ t[i] - t0;
ENDLOOP;
END;
IF finalEndFree
THEN
BEGIN
interval: REAL ¬ t[n];
FOR i:
NAT
IN [l..n]
DO
t[i] ¬ t[i]/interval;
ENDLOOP;
END;
END;
Adjustment:
PROCEDURE [z: Complex.
VEC, f: SplineFunction, t:
REAL]
RETURNS [delta:REAL] =
BEGIN
epsilon: REAL ← 1.0/8388608;
x,dx,ddx,dddx: REAL;
y,dy,ddy,dddy: REAL;
lo,hi: REAL;
derivSqrDist: REAL;
derivDerivSqrDist: REAL;
[x,dx,ddx,dddx,lo,hi] ¬ PiecewiseCubic.EvalAll[f.x, t];
[y,dy,ddy,dddy,lo,hi] ¬ PiecewiseCubic.EvalAll[f.y, t];
derivSqrDist ¬ (x-z.x)*dx + (y-z.y)*dy;
derivDerivSqrDist ¬ dx*dx + dy*dy + (x-z.x)*ddx + (y-z.y)*ddy;
delta ¬ IF derivDerivSqrDist<=0 THEN 0 ELSE -derivSqrDist/derivDerivSqrDist;
IF t+delta>hi THEN delta ← hi+epsilon-t;
IF t+delta<lo THEN delta ← lo-epsilon-t;
END;
Sqr: PROC [x:REAL] RETURNS [REAL] = INLINE{RETURN[x*x]};
AccelTValues:
PUBLIC
PROCEDURE [h: Handle]
RETURNS [r, maxChange:
REAL] =
BEGIN OPEN h;
sqrOldDelta, deltaDotOldDelta: REAL ¬ 0;
r ¬ 0;
maxChange ¬ AdjustTValues[h];
SolveCurve[h]; maxChange ¬ AdjustTValues[h];
FOR i:
NAT
IN [l..n]
DO
sqrOldDelta ¬ sqrOldDelta + Sqr[oldt[i]-oldoldt[i]];
deltaDotOldDelta ¬ deltaDotOldDelta + (t[i]-oldt[i])*(oldt[i]-oldoldt[i]);
ENDLOOP;
IF sqrOldDelta=0 THEN RETURN;
r ¬ deltaDotOldDelta/sqrOldDelta;
IF ABS[r]>=1 THEN RETURN;
maxChange ¬ 0;
FOR i:
NAT
IN [l..n]
DO
delta: REAL ¬ (t[i] - oldt[i])/(1-r);
t[i] ¬ oldt[i] + delta;
IF ABS[delta] > maxChange THEN maxChange ¬ ABS[delta];
ENDLOOP;
END;
PatchesOf:
PUBLIC
PROCEDURE [h: Handle]
RETURNS [x, y: PatchSequence, knots: RealSequence] =
BEGIN
i,nPatches: NAT ¬ 0;
CountPatches: PiecewiseCubic.PieceProc = {nPatches ¬ nPatches+1};
StorePatch: PiecewiseCubic.PieceProc =
BEGIN
interval: REAL ¬ p.domainEnd-p.domainStart;
b: Cubic.Bezier ¬
[b0: [p.initValue, q.initValue],
b1: [p.initValue + p.initSlope*interval/3, q.initValue + q.initSlope*interval/3],
b2: [p.finalValue - p.finalSlope*interval/3, q.finalValue - q.finalSlope*interval/3],
b3: [p.finalValue, q.finalValue]];
c: Cubic.Coeffs ¬ Cubic.BezierToCoeffs[b];
x[i] ¬ [c0:c.c0.x, c1:c.c1.x, c2:c.c2.x, c3:c.c3.x];
y[i] ¬ [c0:c.c0.y, c1:c.c1.y, c2:c.c2.y, c3:c.c3.y];
knots[i] ¬ p.domainStart;
knots[i+1] ¬ p.domainEnd;
i ¬ i + 1;
END;
ValidateCache[h];
PiecewiseCubic.EnumerateCommonPieces[h.patchCache.x,h.patchCache.y,CountPatches];
x ¬ NEW[PatchSequenceRec[nPatches]];
y ¬ NEW[PatchSequenceRec[nPatches]];
knots ¬ NEW[RealSequenceRec[nPatches+1]];
PiecewiseCubic.EnumerateCommonPieces[h.patchCache.x,h.patchCache.y,StorePatch];
END;
ValidateCache:
PROCEDURE [h: Handle] =
BEGIN OPEN h;
IF patchCacheValid THEN RETURN;
patchCache ¬ [PiecewiseCubic.Zero[],PiecewiseCubic.Zero[]];
FOR i:
NAT
DECREASING
IN [0..a.length)
DO
patchCache ¬ [PiecewiseCubic.Combine[1,patchCache.x,a[i],splineBasis[i].x],
PiecewiseCubic.Combine[1,patchCache.y,a[i],splineBasis[i].y]]
ENDLOOP;
patchCacheValid¬TRUE;
END;
PointSlopeBasis:
PUBLIC
PROCEDURE [h: Handle, b: Cubic.Bezier] =
BEGIN OPEN h,PiecewiseCubic;
vx,vy: REAL;
closedCurve ¬ FALSE;
initialEndFree ¬ FALSE;
finalEndFree ¬ FALSE;
patchCacheValid ¬ FALSE;
splineBasis ¬ NEW[SplineBasisRec[8]];
a ¬ NEW[RealSequenceRec[8]];
free ¬ NEW[BooleanSequenceRec[8]];
(a0,a1) and (a6,a7) are the positions of the endpoints
a2 is the velocity component parallel to b1-b0
a3 is the velocity component perpendicular to b1-b0
a4 is the velocity component parallel to b2-b3
a5 is the velocity component perpendicular to b2-b3
splineBasis[0] ¬ [Piece[0,1,1,0,0,0],Zero[]]; a[0] ¬ b.b0.x;
splineBasis[1] ¬ [Zero[],Piece[0,1,1,0,0,0]]; a[1] ¬ b.b0.y;
[[vx,vy]] ¬ Complex.Sub[b.b1,b.b0];
splineBasis[2] ¬ [Combine[3*vx,Piece[0,1,0,1,0,0],0,Zero[]],
Combine[3*vy,Piece[0,1,0,1,0,0],0,Zero[]]]; a[2] ¬ 1;
splineBasis[3] ¬ [Combine[-3*vy,Piece[0,1,0,1,0,0],0,Zero[]],
Combine[3*vx,Piece[0,1,0,1,0,0],0,Zero[]]]; a[3] ¬ 0;
[[vx,vy]] ¬ Complex.Sub[b.b3,b.b2];
splineBasis[4] ¬ [Combine[3*vx,Piece[0,1,0,0,1,0],0,Zero[]],
Combine[3*vy,Piece[0,1,0,0,1,0],0,Zero[]]]; a[4] ¬ 1;
splineBasis[5] ¬ [Combine[-3*vy,Piece[0,1,0,0,1,0],0,Zero[]],
Combine[3*vx,Piece[0,1,0,0,1,0],0,Zero[]]]; a[5] ¬ 0;
splineBasis[6] ¬ [Piece[0,1,0,0,0,1],Zero[]]; a[6] ¬ b.b3.x;
splineBasis[7] ¬ [Zero[],Piece[0,1,0,0,0,1]]; a[7] ¬ b.b3.y;
FOR i:
NAT
IN [0..8)
DO
free[i] ¬ FALSE;
ENDLOOP;
END;
SmoothBasis:
PUBLIC
PROCEDURE [h: Handle, nKnots:
NAT] =
BEGIN OPEN h;
closedCurve ¬ TRUE;
initialEndFree ¬ FALSE;
finalEndFree ¬ FALSE;
patchCacheValid ¬ FALSE;
splineBasis ¬ NEW[SplineBasisRec[4*nKnots]];
a ¬ NEW[RealSequenceRec[4*nKnots]];
free ¬ NEW[BooleanSequenceRec[4*nKnots]];
FOR i:
NAT
IN [0..nKnots)
DO
even: PiecewiseCubic.Handle ¬ EvenBasis[i*1.0/nKnots,1.0/nKnots];
odd: PiecewiseCubic.Handle ¬ OddBasis[i*1.0/nKnots,1.0/nKnots];
splineBasis[4*i] ¬ [even,PiecewiseCubic.Zero[]];
splineBasis[4*i+1] ¬ [PiecewiseCubic.Zero[],even];
splineBasis[4*i+2] ¬ [odd,PiecewiseCubic.Zero[]];
splineBasis[4*i+3] ¬ [PiecewiseCubic.Zero[],odd];
ENDLOOP;
FOR i:
NAT
IN [0..4*nKnots)
DO
a[i] ¬ 0;
free[i] ¬ TRUE;
ENDLOOP;
END;
EvenBasis:
PROCEDURE [center, radius:
REAL]
RETURNS [PiecewiseCubic.Handle] =
BEGIN OPEN PiecewiseCubic;
l: REAL ¬ center-radius;
IF l<0 THEN l¬l+1;
RETURN [Combine[1,Piece[l,l+radius,0,0,0,1],
1,Piece[center,center+radius,1,0,0,0]]]
END;
OddBasis:
PROCEDURE [center, radius:
REAL]
RETURNS [PiecewiseCubic.Handle] =
BEGIN OPEN PiecewiseCubic;
l: REAL ¬ center-radius;
IF l<0 THEN l¬l+1;
RETURN [Combine[1,Piece[l,l+radius,0,0,1,0],
1,Piece[center,center+radius,0,1,0,0]]]
END;
BSplineBasis:
PUBLIC
PROCEDURE [h: Handle, nKnots:
NAT] =
BEGIN OPEN h;
closedCurve ¬ TRUE;
initialEndFree ¬ FALSE;
finalEndFree ¬ FALSE;
patchCacheValid ¬ FALSE;
splineBasis ¬ NEW[SplineBasisRec[2*nKnots]];
a ¬ NEW[RealSequenceRec[2*nKnots]];
free ¬ NEW[BooleanSequenceRec[2*nKnots]];
FOR i:
NAT
IN [0..nKnots)
DO
elem: PiecewiseCubic.Handle ¬ BSplineBasisElement[i,nKnots];
splineBasis[2*i] ¬ [elem,PiecewiseCubic.Zero[]];
splineBasis[2*i+1] ¬ [PiecewiseCubic.Zero[],elem];
ENDLOOP;
FOR i:
NAT
IN [0..2*nKnots)
DO
a[i] ¬ 0;
free[i] ¬ TRUE;
ENDLOOP;
END;
BSplineBasisElement:
PROCEDURE [i, n:
NAT]
RETURNS [f: PiecewiseCubic.Handle] =
BEGIN
k: REAL ¬ i-2;
IF k<0 THEN k¬k+n;
f ¬ PiecewiseCubic.Piece[k/n,(k+1)/n,0,0,0.75*n,0.25];
k¬k+1; IF k>n-1 THEN k¬k-n;
f ¬ PiecewiseCubic.Combine[1,f,1,PiecewiseCubic.Piece[k/n,(k+1)/n,0.25,0.75*n,0,1]];
k¬k+1; IF k>n-1 THEN k¬k-n;
f ¬ PiecewiseCubic.Combine[1,f,1,PiecewiseCubic.Piece[k/n,(k+1)/n,1,0,-0.75*n,0.25]];
k¬k+1; IF k>n-1 THEN k¬k-n;
f ¬ PiecewiseCubic.Combine[1,f,1,PiecewiseCubic.Piece[k/n,(k+1)/n,0.25,-0.75*n,0,0]];
END;
CubicBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
zero: PiecewiseCubic.Handle ¬ PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[-10,10,1,0,0,1];
iden: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[-10,10,-10,1,1,10];
sqr: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[-10,10,100,-20,20,100];
cube: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[-10,10,-1000,300,300,1000];
closedCurve ¬ FALSE;
initialEndFree ¬ TRUE;
finalEndFree ¬ TRUE;
patchCacheValid ¬ FALSE;
splineBasis ¬ NEW[SplineBasisRec[8]];
a ¬ NEW[RealSequenceRec[8]];
free ¬ NEW[BooleanSequenceRec[8]];
splineBasis[0] ¬ [one,zero];
splineBasis[1] ¬ [zero,one];
splineBasis[2] ¬ [iden,zero];
splineBasis[3] ¬ [zero,iden];
splineBasis[4] ¬ [sqr,zero];
splineBasis[5] ¬ [zero,sqr];
splineBasis[6] ¬ [cube,zero];
splineBasis[7] ¬ [zero,cube];
FOR i:
NAT
IN [0..8)
DO
a[i] ¬ 0;
free[i] ¬ TRUE;
ENDLOOP;
END;
CubicPieceBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
zero: PiecewiseCubic.Handle ¬ PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[0,1,1,0,0,1];
iden: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[0,1,0,1,1,1];
sqr: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[0,1,0,0,2,1];
cube: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[0,1,0,0,3,1];
closedCurve ¬ FALSE;
initialEndFree ¬ TRUE;
finalEndFree ¬ TRUE;
patchCacheValid ¬ FALSE;
splineBasis ¬ NEW[SplineBasisRec[8]];
a ¬ NEW[RealSequenceRec[8]];
free ¬ NEW[BooleanSequenceRec[8]];
splineBasis[0] ¬ [one,zero];
splineBasis[1] ¬ [zero,one];
splineBasis[2] ¬ [iden,zero];
splineBasis[3] ¬ [zero,iden];
splineBasis[4] ¬ [sqr,zero];
splineBasis[5] ¬ [zero,sqr];
splineBasis[6] ¬ [cube,zero];
splineBasis[7] ¬ [zero,cube];
FOR i:
NAT
IN [0..8)
DO
a[i] ¬ 0;
free[i] ¬ TRUE;
ENDLOOP;
END;
ParabolaBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
zero: PiecewiseCubic.Handle ¬ PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[-10,10,1,0,0,1];
iden: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[-10,10,-10,1,1,10];
sqr: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[-10,10,100,-20,20,100];
closedCurve ¬ FALSE;
initialEndFree ¬ TRUE;
finalEndFree ¬ TRUE;
patchCacheValid ¬ FALSE;
splineBasis ¬ NEW[SplineBasisRec[6]];
a ¬ NEW[RealSequenceRec[6]];
free ¬ NEW[BooleanSequenceRec[6]];
splineBasis[0] ¬ [one,zero];
splineBasis[1] ¬ [zero,one];
splineBasis[2] ¬ [iden,zero];
splineBasis[3] ¬ [zero,iden];
splineBasis[4] ¬ [sqr,zero];
splineBasis[5] ¬ [zero,sqr];
FOR i:
NAT
IN [0..6)
DO
a[i] ¬ 0;
free[i] ¬ TRUE;
ENDLOOP;
END;
CircleBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
twopi: REAL ¬ 2*3.1415926;
zero: PiecewiseCubic.Handle ¬ PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[0,1,1,0,0,1];
cos: PiecewiseCubic.Handle ¬ CosBasis[];
sin: PiecewiseCubic.Handle ¬ SinBasis[];
closedCurve ¬ TRUE;
initialEndFree ¬ FALSE;
finalEndFree ¬ FALSE;
patchCacheValid ¬ FALSE;
splineBasis ¬ NEW[SplineBasisRec[3]];
a ¬ NEW[RealSequenceRec[3]];
free ¬ NEW[BooleanSequenceRec[3]];
splineBasis[0] ¬ [one,zero];
splineBasis[1] ¬ [zero,one];
splineBasis[2] ¬ [cos,sin];
FOR i:
NAT
IN [0..3)
DO
a[i] ¬ 0;
free[i] ¬ TRUE;
ENDLOOP;
END;
EllipseBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
twopi: REAL ¬ 2*3.1415926;
zero: PiecewiseCubic.Handle ¬ PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ¬ PiecewiseCubic.Piece[0,1,1,0,0,1];
cos: PiecewiseCubic.Handle ¬ CosBasis[];
sin: PiecewiseCubic.Handle ¬ SinBasis[];
closedCurve ¬ TRUE;
initialEndFree ¬ FALSE;
finalEndFree ¬ FALSE;
patchCacheValid ¬ FALSE;
splineBasis ¬ NEW[SplineBasisRec[5]];
a ¬ NEW[RealSequenceRec[5]];
free ¬ NEW[BooleanSequenceRec[5]];
splineBasis[0] ¬ [one,zero];
splineBasis[1] ¬ [zero,one];
splineBasis[2] ¬ [cos,sin];
splineBasis[3] ¬ [sin,zero];
splineBasis[4] ¬ [zero,cos];
FOR i:
NAT
IN [0..5)
DO
a[i] ¬ 0;
free[i] ¬ TRUE;
ENDLOOP;
END;
CosBasis:
PROC
RETURNS[PiecewiseCubic.Handle] =
{RETURN[LIST[[domainStart: 0, domainEnd: 0.25, initValue: 99.98286, initSlope: 1.302177e-2, finalSlope: -662.7453, finalValue: -1.397745e-4], [domainStart: 0.25, domainEnd: 0.5, initValue: -1.397745e-4, initSlope: -662.7453, finalSlope: 2.099349e-3, finalValue: -99.9828], [domainStart: 0.5, domainEnd: 0.75, initValue: -99.9828, initSlope: 2.099349e-3, finalSlope: 662.7528, finalValue: 6.351456e-6], [domainStart: 0.75, domainEnd: 1, initValue: 6.351456e-6, initSlope: 662.7528, finalSlope: 1.302177e-2, finalValue: 99.98286]]]};
SinBasis:
PROC
RETURNS[PiecewiseCubic.Handle] =
{RETURN[LIST[[domainStart: 0, domainEnd: 0.25, initValue: 1.728265e-4, initSlope: 662.7498, finalSlope: -4.030657e-4, finalValue: 99.98302], [domainStart: 0.25, domainEnd: 0.5, initValue: 99.98302, initSlope: -4.030657e-4, finalSlope: -662.7458, finalValue: 1.005444e-4], [domainStart: 0.5, domainEnd: 0.75, initValue: 1.005444e-4, initSlope: -662.7458, finalSlope: 1.04512e-3, finalValue: -99.98254], [domainStart: 0.75, domainEnd: 1, initValue: -99.98254, initSlope: 1.04512e-3, finalSlope: 662.7498, finalValue: 1.728265e-4]]]};
END.
Quick approximations
cos ¬ PiecewiseCubic.Combine[1,PiecewiseCubic.Combine[1,PiecewiseCubic.Combine[1,
PiecewiseCubic.Piece[0,0.25,1,0,-twopi,0],1,
PiecewiseCubic.Piece[0.25,0.5,0,-twopi,0,-1]],1,
PiecewiseCubic.Piece[0.5,0.75,-1,0,twopi,0]],1,
PiecewiseCubic.Piece[0.75,1,0,twopi,0,1]];
sin ¬ PiecewiseCubic.Combine[1,PiecewiseCubic.Combine[1,PiecewiseCubic.Combine[1,
PiecewiseCubic.Piece[0,0.25,0,twopi,0,1],1,
PiecewiseCubic.Piece[0.25,0.5,1,0,-twopi,0]],1,
PiecewiseCubic.Piece[0.5,0.75,0,-twopi,0,-1]],1,
PiecewiseCubic.Piece[0.75,1,-1,0,twopi,0]];