-- 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.