--M. Stone September 26, 1980 5:42 PM
--subdivision, patricks test for change in tangent
-- Last Edited by: Stone, February 9, 1983 4:52 pm
DIRECTORY

 RealFns: FROM "RealFns",
 SplineDefs: FROM "SplineDefs",
 PressLineDefs: FROM "PressLineDefs",
 GriffinMemoryDefs USING [CZone],
 PointDefs: FROM "PointDefs";

PressLines: PROGRAM
IMPORTS RealFns,SplineDefs, GriffinMemoryDefs
EXPORTS PressLineDefs=
BEGIN OPEN PointDefs,PressLineDefs, GriffinMemoryDefs;

--calling proc must START Float and Init Splines

diff,eps,magtstart,W,deltaT: REAL;
prec: INTEGER ← 20;
t,p,start,tstart,next: ObjPt;
firstpt: ObjPt;
end,End0: End;
PutCurve: PutProc;
lastx,lasty,Tol: REAL;
ONE: REAL ← 1; E2: REAL ← 100;
half: REALONE/2;
Vals: TYPE=
RECORD[Fx,Cx,Tx: REAL, Fy,Cy,Ty: REAL];

RABS: PROCEDURE[r: REAL] RETURNS [REAL] = INLINE BEGIN
IF r<0 THEN RETURN[-r] ELSE RETURN[r];
END;

ContourCubic: PUBLIC PROCEDURE[coeffs: SplineDefs.Coeffs, w: REAL, end0,end1: End, putCurve: PutProc]=
BEGIN
p0,p1,t0,t1: ObjPt;
i: CARDINAL;
val0,val1: Vals;
dtanx,dtany: REAL;
zeropt: ObjPt ← [0,0];
FOR i IN [X..Y] DO
 p0[i] ← coeffs.t0[i];
 p1[i] ← coeffs.t0[i]+coeffs.t1[i]+coeffs.t2[i]+coeffs.t3[i];
 t0[i] ← coeffs.t1[i];
 t1[i] ← coeffs.t1[i]+2*coeffs.t2[i]+3*coeffs.t3[i];
ENDLOOP;
--do lines directly
IF coeffs.t3=coeffs.t2 AND coeffs.t2=zeropt THEN BEGIN
 ContourLine[p0,p1,w,end0,end1,putCurve];
RETURN;
END;
firstpt ← p0;
start ← p0;
tstart ← t0;
deltaT ← 0;
magtstart ← tstart[X]*tstart[X]+tstart[Y]*tstart[Y];
--can't run the algorithm with magtstart=0. Is probably a badly parameterized
--Bezier curve. Punt
IF magtstart=0 THEN BEGIN
 ContourLine[p0,p1,w,end0,end1,putCurve];
RETURN;
END;
End0← end0;
W ← w;
eps ← prec/E2; --from Draw
PutCurve ← putCurve;

--set initial values. assume t0=0, t1=1
--F[0], F[1], T[0],T[1], C ← (6at+2b)/2
--d is 1/2, but things will be div by 4 first time thru, so don't divide here
 val0.Fx ← p0[X];
 val0.Tx ← t0[X];
 val0.Cx ← coeffs.t2[X];
 val0.Fy ← p0[Y];
 val0.Ty ← t0[Y];
 val0.Cy ← coeffs.t2[Y];

 lastx ← val1.Fx ← p1[X];
 val1.Tx ← t1[X];
 val1.Cx ← (3*coeffs.t3[X]+coeffs.t2[X]);
 lasty ← val1.Fy ← p1[Y];
 val1.Ty ← t1[Y];
 val1.Cy ← (3*coeffs.t3[Y]+coeffs.t2[Y]);

 dtanx ← (3*coeffs.t3[X]); 
 dtany ← (3*coeffs.t3[Y]);
 Tol ← 2*w;
DrawCurve[val0, val1,dtanx,dtany,1]
END;

--recursive subdivision
DrawCurve: PROCEDURE[val1,val2: Vals, dtanX,dtanY: REAL, thisT: REAL ]=
BEGIN
DO
SELECT TRUE FROM
--is small or flat. check tangent
(RABS[val1.Cx] < half AND RABS[val1.Cy] < half AND RABS[val2.Cx] < half AND RABS[val2.Cy] < half) OR (RABS[val1.Fx-val2.Fx] <=Tol AND RABS[val1.Fy-val2.Fy] <=Tol)=> BEGIN
 deltaT ← deltaT+thisT;
 t[X] ← val2.Tx;
 t[Y] ← val2.Ty;
 diff ← ABS[(t[Y]*tstart[X]-t[X]*tstart[Y])/magtstart];
IF diff< eps AND (RABS[lastx-val2.Fx]>=1 OR RABS[lasty-val2.Fy]>=1)
  THEN RETURN;
 p[X] ← val2.Fx;
 p[Y] ← val2.Fy;
 end ← (IF p=firstpt THEN End0 ELSE round);
--scale the tangents
 next ← t;
 t[X] ← t[X]*deltaT; t[Y] ←t[Y]*deltaT;
 tstart[X] ← tstart[X]*deltaT; tstart[Y] ←tstart[Y]*deltaT;
 ContourCurve[start,p,tstart,t,W,end,round,PutCurve];
 start ← p;
 tstart ← next;
 deltaT ← 0;
 magtstart ← tstart[X]*tstart[X]+tstart[Y]*tstart[Y];
RETURN;
END;
ENDCASE => BEGIN
 valc: Vals;
 thisT ← thisT/2;
 val1.Cx ← val1.Cx/4; val2.Cx ← val2.Cx/4;
 valc.Cx ← (val1.Cx+val2.Cx)/2;
 valc.Fx ← ((val1.Fx+val2.Fx)/2) - valc.Cx;

 val1.Cy ← val1.Cy/4; val2.Cy ← val2.Cy/4;
 valc.Cy ← (val1.Cy+val2.Cy)/2;
 valc.Fy ← ((val1.Fy+val2.Fy)/2) - valc.Cy;

 dtanX ← dtanX/4; dtanY ← dtanY/4;
 valc.Tx ← ((val1.Tx+val2.Tx)/2) - dtanX;
 valc.Ty ← ((val1.Ty+val2.Ty)/2) - dtanY;
 DrawCurve[val1,valc,dtanX,dtanY,thisT];
 val1 ←valc; --instead of a second call to DrawCurve
END;
ENDLOOP;
END;

ContourCurve: PROCEDURE[p0,p1,t0,t1: ObjPt, w: REAL, end0,end1: End, putCurve: PutProc]=
BEGIN
length: REAL;
nx0,ny0,nx1,ny1: REAL;
coeffs: SplineDefs.Coeffs;
ncoeffs: SplineDefs.CoeffsSequence; --used by makespline
pt0,pt1: ObjPt;
knots: SplineDefs.KnotSequence ← CZone.NEW[SplineDefs.KnotSequenceRec[4]];
length ← RealFns.SqRt[t0[X]*t0[X]+t0[Y]*t0[Y]];
IF length = 0 THEN BEGIN t0[X] ← 1; t0[Y] ← 0; length ← 1; END;
ny0 ← w*t0[X]/length;
nx0 ← -w*t0[Y]/length;
--[nx0,ny0] is a vector normal to t0, w long
length ← RealFns.SqRt[t1[X]*t1[X]+t1[Y]*t1[Y]];
IF length = 0 THEN BEGIN t1[X] ← 1; t1[Y] ← 0; length ← 1; END;
ny1 ← w*t1[X]/length;
nx1 ← -w*t1[Y]/length;
--[nx1,ny1] is a vector normal to t1, w long
pt0 ← [p0[X]+nx0,p0[Y]+ny0];
IF end0=square THEN pt0 ← [pt0[X]-ny0,pt0[Y]+nx0]; --move it back a unit
pt1 ← [p1[X]+nx1,p1[Y]+ny1];
IF end1=square THEN pt1 ← [pt1[X]+ny1,pt1[Y]-nx1]; --move it up a unit
coeffs ← MakeCoeffs[pt0,pt1,t0,t1];
putCurve[coeffs,TRUE]; --start new object

--first end
SELECT end1 FROM
 flat,square => BEGIN
  coeffs.t0 ← pt1;
  coeffs.t1 ← [-2*nx1,-2*ny1];
  coeffs.t3 ← coeffs.t2 ← [0,0];
  pt1 ← [pt1[X]+coeffs.t1[X],pt1[Y]+coeffs.t1[Y]];
  END;
 round => BEGIN
  knots[0] ← pt1;
  knots[3] ← [p1[X]-nx1,p1[Y]-ny1];
  knots[1] ← [knots[0][X]+ny1*4/3,knots[0][Y]-nx1*4/3];
  knots[2] ← [knots[3][X]+ny1*4/3,knots[3][Y]-nx1*4/3];
  ncoeffs ← SplineDefs.MakeSpline[knots,bezier];
  coeffs ← ncoeffs[0];
  pt1 ← knots[3];
  END;
ENDCASE;
putCurve[coeffs,FALSE];

--second side
pt0 ← [p0[X]-nx0,p0[Y]-ny0];
IF end0=square THEN pt0 ← [pt0[X]-ny0,pt0[Y]+nx0]; --move it back a unit
t0 ← [-t0[X],-t0[Y]];
t1 ← [-t1[X],-t1[Y]];
coeffs ← MakeCoeffs[pt1,pt0,t1,t0];
putCurve[coeffs,FALSE];

--second end
SELECT end0 FROM
 flat,square => BEGIN
  coeffs.t0 ← pt0;
  coeffs.t1 ← [2*nx0,2*ny0];
  coeffs.t3 ← coeffs.t2 ← [0,0];
  END;
 round => BEGIN
  knots[0] ← pt0;
  knots[3] ← [p0[X]+nx0,p0[Y]+ny0];
  knots[1] ← [knots[0][X]-ny0*4/3,knots[0][Y]+nx0*4/3];
  knots[2] ← [knots[3][X]-ny0*4/3,knots[3][Y]+nx0*4/3];
  ncoeffs ← SplineDefs.MakeSpline[knots,bezier];
  coeffs ← ncoeffs[0];
  END;
ENDCASE;
putCurve[coeffs,FALSE];
END;

--2 points and 2 tangents to make a curve (Hermite)
MakeCoeffs: PROCEDURE[p0,p1,t0,t1: ObjPt] RETURNS[coeffs: SplineDefs.Coeffs] =
BEGIN
i: CARDINAL;
FOR i IN [X..Y] DO
 coeffs.t3[i] ← 2*p0[i]-2*p1[i]+t0[i]+t1[i];
 coeffs.t2[i] ← -3*p0[i]+3*p1[i]-2*t0[i]-t1[i];
 coeffs.t1[i] ← t0[i];
 coeffs.t0[i] ← p0[i];
ENDLOOP;
END;

ContourLine: PUBLIC PROCEDURE[p0,p1: ObjPt, w: REAL, end0,end1: End, putCurve: PutProc]=
BEGIN
length: REAL;
dx,dy,nx,ny: REAL;
coeffs: SplineDefs.Coeffs;
ncoeffs: SplineDefs.CoeffsSequence; --used by makespline
pt: ObjPt;
knots: SplineDefs.KnotSequence ← CZone.NEW[SplineDefs.KnotSequenceRec[4]];
dx ← p1[X]-p0[X];
dy ← p1[Y]-p0[Y];
length ← RealFns.SqRt[dx*dx+dy*dy];
IF length = 0 THEN BEGIN dx ← 1; dy ← 0; length ← 1; END;
ny ← w*dx/length;
nx ← -w*dy/length;
--[nx,ny] is a vector normal to the line, w long
coeffs ← [[0,0],[0,0],[0,0],[0,0]];
coeffs.t0 ← [p0[X]+nx,p0[Y]+ny];
IF end0=square THEN coeffs.t0 ← [coeffs.t0[X]-ny,coeffs.t0[Y]+nx]; --move it back a unit
pt ← [p1[X]+nx,p1[Y]+ny];
IF end1=square THEN pt ← [pt[X]+ny,pt[Y]-nx]; --move it up a unit
coeffs.t1 ← [pt[X]-coeffs.t0[X], pt[Y]-coeffs.t0[Y]];
putCurve[coeffs,TRUE]; --start new object

--first end
SELECT end1 FROM
 flat,square => BEGIN
  coeffs.t0 ← pt;
  coeffs.t1 ← [-2*nx,-2*ny];
  pt ← [pt[X]+coeffs.t1[X],pt[Y]+coeffs.t1[Y]];
  END;
 round => BEGIN
  knots[0] ← pt;
  knots[3] ← [p1[X]-nx,p1[Y]-ny];
  knots[1] ← [knots[0][X]+ny*4/3,knots[0][Y]-nx*4/3];
  knots[2] ← [knots[3][X]+ny*4/3,knots[3][Y]-nx*4/3];
  ncoeffs ← SplineDefs.MakeSpline[knots,bezier];
  coeffs ← ncoeffs[0];
  pt ← knots[3];
  END;
ENDCASE;
putCurve[coeffs,FALSE];

--second side
coeffs ← [[0,0],[0,0],[0,0],[0,0]];
coeffs.t0 ← pt;
pt ← [p0[X]-nx,p0[Y]-ny];
IF end0=square THEN pt ← [pt[X]-ny,pt[Y]+nx]; --move it back a unit
coeffs.t1 ← [pt[X]-coeffs.t0[X], pt[Y]-coeffs.t0[Y]];
putCurve[coeffs,FALSE];

--second end
SELECT end0 FROM
 flat,square => BEGIN
  coeffs.t0 ← pt;
  coeffs.t1 ← [2*nx,2*ny];
  END;
 round => BEGIN
  knots[0] ← pt;
  knots[3] ← [p0[X]+nx,p0[Y]+ny];
  knots[1] ← [knots[0][X]-ny*4/3,knots[0][Y]+nx*4/3];
  knots[2] ← [knots[3][X]-ny*4/3,knots[3][Y]+nx*4/3];
  ncoeffs ← SplineDefs.MakeSpline[knots,bezier];
  coeffs ← ncoeffs[0];
  END;
ENDCASE;
putCurve[coeffs,FALSE];
END;

END.