-- NodeImpl.mesa
--last edited by Michael Plass August 30, 1982 2:06 pm
DIRECTORY Complex, Cubic, Curve, DynFit, FitDisplayUtilities, JaMFnsDefs, LSPiece, Seq, Vector, RealFns;
NodeImpl: PROGRAM IMPORTS Complex, Cubic, Curve, DynFit, FitDisplayUtilities, JaMFnsDefs, LSPiece, Vector, RealFns =
BEGIN OPEN JaMFnsDefs;
OpenCurve: PROC = {
Curve.AddNode[Curve.defaultHandle, 0];
Curve.AddNode[Curve.defaultHandle, Curve.CurrentSamples[Curve.defaultHandle].length-1];
};
SetTangents: PROC = {-- range err maxit => . Sets tangents at nodes by locally fitting 2*range + 1 samples
maxit: NAT ← PopInteger[];
err: REAL ← GetReal[];
range: NAT ← PopInteger[];
node,cusp: Seq.NatSequence;
z: Seq.ComplexSequence ← Curve.CurrentSamples[Curve.defaultHandle];
n: NAT ← z.length;
t: Seq.RealSequence ← NEW[Seq.RealSequenceRec[n]];
bezier: Cubic.Bezier;
nextCusp,cindex: INT ← 0;
forward: PROC[pt: NAT, range: INTEGER] RETURNS[INTEGER] = {
np,nc: INTEGER;
np ← (pt+range) MOD z.length;
IF cusp.length = 0 THEN RETURN[np]
ELSE {
i: INTEGER ← IF pt=nextCusp THEN (cindex+1) MOD cusp.length ELSE cindex;
nc ← cusp[i];
RETURN[MIN[np,nc]];
};
};
back: PROC[pt: NAT, range: INTEGER] RETURNS[INTEGER] = {
np,nc: INTEGER;
np ← pt-range;
IF np<0 THEN np ← z.length+np;
IF cusp.length=0 THEN RETURN[np]
ELSE {
nc ← IF cindex-1 <0 THEN cusp[cusp.length-1] ELSE cusp[cindex-1];
RETURN[MAX[np,nc]];
};
};
[nodes:node] ← Curve.CurrentNodes[Curve.defaultHandle];
[cusps: cusp] ← Curve.CurrentCusps[Curve.defaultHandle];
IF cusp.length=0 THEN nextCusp ← -1 ELSE nextCusp ← cusp[0];
FOR i: NAT IN [0..node.length) DO
q: NAT ← node[i];
in: Complex.Vec;
IF q=nextCusp THEN {
p: INTEGER ← back[q,2*range];
r: INTEGER ← forward[q,2*range];
[b:bezier] ← LSPiece.FitPiece[z: z, t: t, from: p, thru: q,
eps: .0001, maxd: err, maxit: maxit, initFree: TRUE, finalFree: TRUE];
in ← TangentVector[bezier,t[q]];
[b:bezier] ← LSPiece.FitPiece[z: z, t: t, from: q, thru: r,
eps: .0001, maxd: err, maxit: maxit, initFree: TRUE, finalFree: TRUE];
Curve.AddCusp[Curve.defaultHandle, q,in,TangentVector[bezier,t[q]]];
cindex ← cindex+1;
nextCusp ← IF cindex < cusp.length THEN cusp[cindex] ELSE -1;
}
ELSE {
p: INTEGER ← back[q,range];
r: INTEGER ← forward[q,range];
[b:bezier] ← LSPiece.FitPiece[z: z, t: t, from: p, thru: r,
eps: .0001, maxd: err, maxit: maxit, initFree: TRUE, finalFree: TRUE];
Curve.AddNode[Curve.defaultHandle, q,TangentVector[bezier,t[q]]];
};
IF GetJaMBreak[] THEN EXIT;
ENDLOOP;
};
CubicTangents: PROC = {-- err maxit => . Sets tangents at nodes by fitting between neighboring nodes.
maxit: NAT ← PopInteger[];
err: REAL ← GetReal[];
node: Seq.NatSequence;
z: Seq.ComplexSequence ← Curve.CurrentSamples[Curve.defaultHandle];
n: NAT ← z.length;
t: Seq.RealSequence ← NEW[Seq.RealSequenceRec[n]];
bezier, oldBezier: Cubic.Bezier;
[nodes:node] ← Curve.CurrentNodes[Curve.defaultHandle];
oldBezier ← [[-1,-1],[-1,-1],[-1,-1],[-1,-1]];
FOR i: NAT IN [0..node.length) DO
p: NAT ← node[IF i=0 THEN node.length-2 ELSE i-1];
q: NAT ← node[i];
r: NAT ← node[IF i=node.length-1 THEN 1 ELSE i+1];
e: REAL;
nPoints: NAT ← IF r>p THEN r-p+1 ELSE z.length-p+r+1;
IF nPoints < minPoints THEN
Curve.AddNode[Curve.defaultHandle, q, CheapTangent[z[p], z[q], z[r], maxSin]]
ELSE {
[b:bezier, maxDev: e] ← LSPiece.FitPiece[z: z, t: t, from: p, thru: r,
eps: .0001, maxd: err/4, maxit: maxit, initFree: freeEnds, finalFree: freeEnds];
Curve.AddNode[Curve.defaultHandle, q, IF e>err THEN [0,0] ELSE TangentVector[bezier,t[q]]];
FitDisplayUtilities.ShowBezierInverted[oldBezier];
FitDisplayUtilities.ShowBezierInverted[bezier];
oldBezier ← bezier;
IF GetJaMBreak[] THEN EXIT;
};
ENDLOOP;
FitDisplayUtilities.ShowBezierInverted[oldBezier];
};
freeEnds: BOOLEAN ← TRUE;
minPoints: NAT ← 4;
maxSin: REAL ← 0.5;
TangentVector: PROC [b: Cubic.Bezier, t: REAL] RETURNS [tang: Complex.Vec] =
BEGIN
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;
END;
DynNodes: PROC = {-- penalty => . Finds nodes using DynFit
endpoints: DynFit.Segments;
[segments:endpoints] ← DynFit.FitSegments[Curve.CurrentSamples[Curve.defaultHandle], GetReal[]]; Curve.ResetNodes[Curve.defaultHandle];
Curve.AddNode[Curve.defaultHandle, 0];
FOR i:NAT IN [0..endpoints.length) DO
Curve.AddNode[Curve.defaultHandle, endpoints[i]]
ENDLOOP;
};
QuickTangents: PROC = {
z: Seq.ComplexSequence ← Curve.CurrentSamples[Curve.defaultHandle];
node: Seq.NatSequence ← Curve.CurrentNodes[Curve.defaultHandle].nodes;
maxAngle: REAL ← GetReal[];
maxSin: REAL ← RealFns.SinDeg[maxAngle];
initTan: Complex.Vec ←
IF z[node[0]] = z[node[node.length-1]] THEN CheapTangent[z[node[node.length-2]], z[node[0]], z[node[1]], maxSin] ELSE [0,0];
Curve.AddNode[Curve.defaultHandle, node[0], initTan];
FOR i: NAT IN [1..node.length-1) DO
Curve.AddNode[Curve.defaultHandle, node[i], CheapTangent[z[node[i-1]], z[node[i]], z[node[i+1]], maxSin]];
ENDLOOP;
Curve.AddNode[Curve.defaultHandle, node[node.length-1], initTan];
};
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];
};
SquareTangents: PROC = {
z: Seq.ComplexSequence ← Curve.CurrentSamples[Curve.defaultHandle];
node: Seq.NatSequence ← Curve.CurrentNodes[Curve.defaultHandle].nodes;
maxAngle: REAL ← GetReal[];
maxSin: REAL ← RealFns.SinDeg[maxAngle];
initTan: Complex.Vec ←
IF z[node[0]] = z[node[node.length-1]] THEN SquareTangent[z[node[node.length-2]], z[node[0]], z[node[1]], maxSin] ELSE [0,0];
Curve.AddNode[Curve.defaultHandle, node[0], initTan];
FOR i: NAT IN [1..node.length-1) DO
Curve.AddNode[Curve.defaultHandle, node[i], SquareTangent[z[node[i-1]], z[node[i]], z[node[i+1]], maxSin]];
ENDLOOP;
Curve.AddNode[Curve.defaultHandle, node[node.length-1], initTan];
};
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]]];
};
JaMFnsDefs.Register[".opencurve"L,OpenCurve]; -- => . Sets tangentless nodes at endpoints
JaMFnsDefs.Register[".settangents"L,SetTangents];-- range maxdev maxit => . Sets tangents at nodes by locally fitting 2*range + 1 samples
JaMFnsDefs.Register[".cubictangents"L,CubicTangents]; -- err maxit => . Sets tangents at nodes by fitting between neighboring nodes.
JaMFnsDefs.Register[".dynnodes"L,DynNodes]; -- penalty => . Finds nodes using DynFit
JaMFnsDefs.Register[".quicktangents"L,QuickTangents]; -- maxangle => . computes tangents by differencing neighbors
JaMFnsDefs.Register[".squaretangents"L,SquareTangents]; -- maxangle => . computes tangents, weighting the longer edges more.
END.
Michael Plass, August 30, 1982 2:06 pm. Added CubicTangents.