-- Analyze.mesa
-- M. Stone November 3, 1983 10:16 am 
-- M. Plass March 12, 1982 11:21 am
DIRECTORY
	Graphics,
	Real,
	Seq,
	Runtime,
	RealFns,
	LSFit,
	ConvertUnsafe USING [ToRope],
	DynFit,
	Vector USING[Sub,Cross,Dot,Mag,Vec,Mul,Det,Matrix],
	CurveDefs,
	TJaMGraphics USING[Painter],
	AnalyzeDefs,
	Cubic USING[Coeffs,Bezier,CoeffsToBezier],
	JaMFnsDefs USING [PushReal, PushInteger, PopInteger,Register, GetReal,PopBoolean,PopString];

Analyze: PROGRAM
	IMPORTS RealFns, TJaMGraphics, Graphics, Vector, JaMFnsDefs, CurveDefs, LSFit, DynFit, Runtime, ConvertUnsafe, Real, Cubic
	EXPORTS CurveDefs,AnalyzeDefs  = {
	OPEN JaMFnsDefs, Vector,CurveDefs;


--JaM registered procedures

BaseLine: PROC = {
	baseLine ← ~baseLine;
	};

SAOff: PROC = {saOff.y ← GetReal[];saOff.x ← GetReal[]};
VSOff: PROC = {vOff ← GetReal[];sOff ← GetReal[]};

RDrawLocs: PROC = {
	DrawLocs[];
	};
SaScale: PROC = {saScale ← GetReal[];};
VScale: PROC = { vScale ← GetReal[];};
SScale: PROC = { sScale ← GetReal[];};

RLocateVS: PROC = {
	y: REAL ← GetReal[];
	x: REAL ← GetReal[];
	LocateVS[x,y];
	};

RLocateXY: PROC = {
	v: REAL ← GetReal[];
	s: REAL ← GetReal[];
	LocateXY[s,v];
	};

RPlotValues: PROCEDURE = {
	PlotValues[values];
	lastValues ← values;
	};

RPlotDValues: PROCEDURE = {
	PlotValues[dValues];
	lastValues ← dValues;
	};

RFilterValues: PROCEDURE = {
	FilterValues[values];
	};

FilterDValues: PROCEDURE = {
	FilterValues[dValues];
	};

DValues: PROC = {
	i: INTEGER ← 0;
	n: INTEGER ← values.l;
	IF n<2 THEN RETURN;
	maxDValue ← Real.SmallestNormalizedNumber;
	minDValue ← Real.LargestNumber;
	dValues ← NEW[ValueSeq[n]];
	FOR i IN [0..n-1) DO
--dvalue.v is difference between adjacent values.  dvalue.s is half way between
--the s of adjacent values
		dValues[i].v ← values[i+1].v-values[i].v;
		dValues[i].s ← (values[i+1].s+values[i].s)/2;
		IF dValues[i].v >maxDValue THEN maxDValue ← dValues[i].v;
		IF dValues[i].v <minDValue THEN minDValue ← dValues[i].v;
		ENDLOOP;
	dValues[n-1].v ← values[0].v-values[n-1].v;
--don't have the information to do this calculation exactly right
	dValues[n-1].s ← values[n-1].s;
	IF dValues[n-1].v >maxDValue THEN maxDValue ← dValues[n-1].v;
	IF dValues[n-1].v <minDValue THEN minDValue ← dValues[n-1].v;
	PushReal[maxDValue];
	PushReal[minDValue];
	};

DValuesToValues: PROCEDURE = {
	values ← dValues;
	dValues ← NIL;
	maxValue ← maxDValue;
	minValue ← minDValue;
	};

ValuesToDValues: PROCEDURE = {
	dValues ← values;
	values ← NIL;
	maxDValue ← maxValue;
	minDValue ← minValue;
	};

ValuesToSamples: PROC = {
	sa: SampleHandle ← sList.samples;
	FOR i: INTEGER IN [0..values.l) DO
	  IF sa=NIL THEN EXIT;
	  sa.value ← values[i].v;
	  sa ← sa.next;
	  ENDLOOP;
	};

SamplesToValues: PROCEDURE = {
	sa: SampleHandle ← sList.samples;
	i: INTEGER ← 0;
	values ← NEW[ValueSeq[sList.nsamples]];
	maxValue ← Real.SmallestNormalizedNumber;
	minValue ← Real.LargestNumber;
	FOR i IN [0..sList.nsamples) DO
		values[i] ← [s: sa.s, v: sa.value];
		IF values[i].v >maxValue THEN maxValue ← values[i].v;
		IF values[i].v <minValue THEN minValue ← values[i].v;
		sa ← sa.next;
		ENDLOOP;
	};

SamplesToSequence: PROCEDURE [s: CurveDefs.SampleHandle]
    RETURNS [z: Seq.ComplexSequence] =
    BEGIN
    sa: CurveDefs.SampleHandle;
    i: NAT;
    n:NAT ← 0;
    FOR sa:CurveDefs.SampleHandle ← s, sa.next UNTIL sa=NIL DO
        n ← n+1;
	ENDLOOP;
    z ← NEW[Seq.ComplexSequenceRec[n]];
    i ← 0;
    FOR sa←s, sa.next UNTIL sa=NIL DO
		z[i] ← sa.xy;
		i←i+1;
		ENDLOOP;
    END;

--JaM setting of analysis procedures

SetCos: PROCEDURE = {
	sampleProc ← CheckCos;
	};

SetDX: PROCEDURE = {
	sampleProc ← CheckDx;
	};

SetDY: PROCEDURE = {
	sampleProc ← CheckDy;
	};

SetMag: PROCEDURE = {
	sampleProc ← CheckMag;
	};

SetAngle: PROCEDURE = {
	sampleProc ← CheckAngle;
	};

SetDxAngle: PROCEDURE = {
	sampleProc ← CheckDxAngle;
	};

SetDyAngle: PROCEDURE = {
	sampleProc ← CheckDyAngle;
	};

SetArcLength: PROCEDURE = {
	sampleProc ← ArcLength;
	};

SetCurvature: PROCEDURE = {
	sampleProc ← CheckCurvature;
	};

SetCCurve: PROCEDURE = {
	sampleProc ← CCurve;
	};

SetCCurveAngle: PROCEDURE = {
	sampleProc ← CCurveAngle;
	};

--analysis procedures

SampleProc: TYPE = PROC[sa: SampleHandle];
sampleProc: SampleProc ← CheckAngle;

Prev: PROC [sa: SampleHandle] RETURNS[SampleHandle] = {
	IF sa.prev#NIL THEN RETURN[sa.prev];
	IF sList.closed THEN RETURN[sList.lastSample] ELSE RETURN[sa];
	};

Next: PROC [sa: SampleHandle] RETURNS[SampleHandle] = {
	IF sa.next#NIL THEN RETURN[sa.next];
	IF sList.closed THEN RETURN[sList.samples] ELSE RETURN[sa];
	};
CheckCos: SampleProc = {
	v0: Vec ← Sub[sa.xy,Prev[sa].xy];
	v1: Vec ← Sub[Next[sa].xy,sa.xy];
	value: REAL ← Dot[v0,v1];
	sa.value ← value/(Mag[v0]*Mag[v1]);
	};

GetSin: SampleProc = {
	v0: Vec ← Sub[sa.xy,Prev[sa].xy];
	v1: Vec ← Sub[Next[sa].xy,sa.xy];
	value: REAL  ← Cross[v0,v1];
	sa.value ← value/(Mag[v0]*Mag[v1]);
	};

CheckMag: SampleProc = {
	sa.value ←sa.s;
	};

CheckAngle: SampleProc = {
	x,y: REAL;
	CheckCos[sa];
	x ← sa.value;
	GetSin[sa];
	y ← sa.value;
	sa.value ← RealFns.ArcTanDeg[y,x];
	};

CheckDxAngle: SampleProc = {
	dx,ds: REAL;
	dx ← sa.xy.x-Prev[sa].xy.x;
	ds ← sa.s-Prev[sa].s;
	sa.value ← RealFns.ArcTanDeg[dx,ds];
	};

CheckDyAngle: SampleProc = {
	dy,ds: REAL;
	dy ← sa.xy.y-Prev[sa].xy.y;
	ds ← sa.s-Prev[sa].s;
	sa.value ← RealFns.ArcTanDeg[dy,ds];
	};

--this is angle/(m0+m1)/2  in degrees
CheckCurvature: SampleProc = {
	v0: Vec ← Sub[sa.xy,Prev[sa].xy];
	v1: Vec ← Sub[Next[sa].xy,sa.xy];
	m0: REAL ← ABS[sa.s-Prev[sa].s];
	m1: REAL ← ABS[Next[sa].s-sa.s];
	x: REAL  ← Dot[v0,v1];
	y: REAL  ← Cross[v0,v1];
	sa.value ← RealFns.ArcTanDeg[y,x];
	sa.value ← .5*sa.value/(m0+m1);
	};

--this is angle/m0+m1/2  in radians
CCurveAngle: SampleProc = {
	v0: Vec ← Sub[sa.xy,Prev[sa].xy];
	v1: Vec ← Sub[Next[sa].xy,sa.xy];
	m0: REAL ← ABS[sa.s-Prev[sa].s];
	m1: REAL ← ABS[Next[sa].s-sa.s];
	x: REAL ← Dot[v0,v1];
	y: REAL ← Cross[v0,v1];
	sa.value ← RealFns.ArcTan[y,x];
	sa.value ← .5*sa.value/(m0+m1);
	};

signal: BOOLEAN ← FALSE;
debug: BOOLEAN ← FALSE;
Debug: PROC = {debug ← ~debug};
CCurve: SampleProc = {
	c,k: Vec;
	valid: BOOLEAN;
	r: REAL;
	[c,valid] ← FindCircleCenter[Prev[sa].xy,sa.xy,Next[sa].xy];
	IF signal AND ~valid THEN Runtime.CallDebugger["points too close"];
	k ← Sub[c,sa.xy];
	r ← Mag[k];
	CCurveAngle[sa];	--compute the sign of the curvature
	IF sa.value<0 THEN r ← -r;
	IF ~valid OR r=0 THEN sa.value ← 0 ELSE sa.value ← 1/r;
	IF debug THEN {
	  Paint: PROC[dc: Graphics.Context] = {
	     xy: Vec ← XFormSamplePt[sa.xy];
	     Graphics.MoveTo[dc,xy.x,xy.y];
	     xy ← XFormSamplePt[c];
	     Graphics.DrawTo[dc,xy.x,xy.y];
	     };
	  TJaMGraphics.Painter[Paint];
	  };
	};

FindCircleCenter: PROCEDURE[p0,p1,p2:Vector.Vec] RETURNS [o:Vector.Vec,valid:BOOLEAN]=
	{d1,d2:Vec;
	d20,d21,det:REAL;
	x20,y20,x21,y21,x22,y22:REAL;
	eps:REAL←0.001;
	d1←Mul[Sub[p0,p1],2];
	d2←Mul[Sub[p1,p2],2];
	x20←p0.x*p0.x; y20←p0.y*p0.y;
	x21←p1.x*p1.x; y21←p1.y*p1.y;
	x22←p2.x*p2.x; y22←p2.y*p2.y;
	d20←x20-x21+y20-y21; 
	d21←x21-x22+y21-y22;
	det←Det[Matrix[d1.x,d1.y,d2.x,d2.y]];
	IF ABS [det] < eps THEN
	{RETURN[Vec[0,0],FALSE]} ELSE
	{o.x←(d20*d2.y-d21*d1.y)/det; o.y←(d1.x*d21-d2.x*d20)/det;
	RETURN[o,TRUE];};
	};


ArcLength: SampleProc = {
	sa.value ← (Next[sa].s+sa.s)/2;
	};

CheckDx: SampleProc = {
	sa.value ← (sa.xy.x-Prev[sa].xy.x)/(sa.s-Prev[sa].s);
	};

CheckDy: SampleProc = {
	sa.value ← (sa.xy.y-Prev[sa].xy.y)/(sa.s-Prev[sa].s);
	};

--plotting routines for Values sequences

sScale: REAL ← 1.5;
vScale: REAL ← 60;
sOff,vOff: REAL ← 0;
saScale: REAL ← 4;
saOff: Vec ← [0,0];
nLocs: INTEGER ← 20;
drawContour: BOOLEAN ← TRUE;
absValues: BOOLEAN ← FALSE;
baseLine: BOOLEAN ← TRUE;
values,dValues,lastValues: Values ← NIL;
maxValue,minValue,maxDValue,minDValue: REAL ← 0;

GetS: PROC[val: Values,index: INTEGER] RETURNS[s: REAL] = {
	IF index >= values.l THEN RETURN[0];
	s ← sScale*val[index].s+sOff;
	IF absValues THEN s ← ABS[s];
	};

GetValue: PROC[val: Values,index: INTEGER] RETURNS[value: REAL] = {
	IF index >= values.l THEN RETURN[0];
	value ← vScale*val[index].v+vOff;
	IF absValues THEN value ← ABS[value];
	};

XFormSamplePt: PROC[xy: Vec] RETURNS[p: Vec] = {
	p.x ←saScale*xy.x+saOff.x;
	p.y ←saScale*xy.y+saOff.y;
	};


PlotValues: PROC[val: Values] = {
    Paint: PROC[dc: Graphics.Context] = {OPEN Graphics;
	length: INTEGER;
	IF val=NIL THEN RETURN;
	length ← val.l;
	MoveTo[dc,GetS[val,0],GetValue[val,0]];
	FOR i: INTEGER IN [1..length) DO DrawTo[dc,GetS[val,i],GetValue[val,i]]; ENDLOOP;
	IF baseLine THEN {MoveTo[dc,sOff, vOff]; DrawTo[dc,GetS[val,length-1],vOff]};
	};
    TJaMGraphics.Painter[Paint];
    };

DrawLocs: PROC = {
    Paint: PROC[dc: Graphics.Context] = {OPEN Graphics;
	i: INTEGER ← 0;
	v: Vec;
	drawlocs: PROC[sa: SampleHandle] = {
		IF i MOD nLocs = 0 THEN {
			v ← XFormSamplePt[sa.xy];
			MoveTo[dc,v.x,v.y];
			DrawTo[dc,GetS[lastValues,i],GetValue[lastValues,i]];
			};
		i ← i+1;
		};
	ForAllSamples[drawlocs];
	};
    TJaMGraphics.Painter[Paint];
    };

LocateVS: PROC[x,y: REAL] = {
	s: REAL ← -1;
	tol: REAL ← 2/saScale;
	findS: PROC [sa: SampleHandle] = {
		IF x IN [sa.xy.x-tol..sa.xy.x+tol] AND y IN [sa.xy.y-tol..sa.xy.y+tol] THEN s ← sa.s;
		};
        Paint: PROC[dc: Graphics.Context] = {OPEN Graphics;
	  MoveTo[dc,x,y];
--find s by searching the contour for xy, then find v by searching the values
	  x ← (x-saOff.x)/saScale;
	  y ← (y-saOff.y)/saScale;
	  ForAllSamples[findS];
	  IF s=-1 THEN RETURN;	--not found
	  FOR i: INTEGER IN [0..lastValues.l) DO
		IF lastValues[i].s=s THEN DrawTo[dc,GetS[lastValues,i],GetValue[lastValues,i]];
		ENDLOOP;
	  };
	TJaMGraphics.Painter[Paint];
	};

LocateXY: PROC[s,v: REAL] = {
	xy: Vec ← [-10000,-10000];
	pt: Vec ← [0,0];
	vs: Vec ← [s,v];
	findXY: PROC [sa: SampleHandle] = {
		IF s IN [Prev[sa].s..Next[sa].s] THEN xy ← sa.xy;
		};
        Paint: PROC[dc: Graphics.Context] = {OPEN Graphics;
	  MoveTo[dc,s,v];
	  s ← (s-sOff)/sScale;
	  v ← (v-vOff)/vScale;
	  JaMFnsDefs.PushReal[s];
	  JaMFnsDefs.PushReal[v];
--find xy by searching the contour for s
	  ForAllSamples[findXY];
	  IF xy=[-10000,-10000] THEN RETURN;
	  pt ← XFormSamplePt[xy];
	  DrawTo[dc,pt.x,pt.y];
	  PushPoint[vs];
	  PushPoint[pt];
	  };
	TJaMGraphics.Painter[Paint];
	};

DrawContour: PROC = {OPEN Graphics;
    Paint: PROC[dc: Graphics.Context] = {
	v: Vec;
	drawSample: PROC[sa: SampleHandle] = {
		v ← XFormSamplePt[sa.xy];
		DrawTo[dc,v.x,v.y];
		};
	IF sList=NIL THEN RETURN;
	v ← XFormSamplePt[sList.samples.xy];
	MoveTo[dc,v.x,v.y];
	ForAllSamples[drawSample];
	};
    TJaMGraphics.Painter[Paint];
    };

abs: BOOLEAN ← TRUE;
SetAbs: PROC = {abs ← JaMFnsDefs.PopBoolean[]};

ClearKnots: PROC = {cps.samples ← cps.lastSample ← NIL; nknots ← 0};
PickSumKnots: PROC = {
	sa: SampleHandle ← sList.samples.next;	--first point is duplicated on end
	nextknot: SampleHandle ← cps.samples;
	sum: REAL ← 0;
	addKnot: PROC[sa: SampleHandle] = {
	    NewCPoint[sa];
	    sum ← 0;
	    };
	eps: REAL ← GetReal[];	--change in angle in degrees
	cps.arcLen ← sList.arcLen;
	addKnot[sa];
	UNTIL sa=NIL DO
--may have preset knots in list
	    IF nextknot#NIL AND sa.s>nextknot.s THEN {sum ← 0; nextknot ← nextknot.next};
	    sum ← sum + (IF abs THEN ABS[sa.value] ELSE sa.value);
	    IF ABS[sum]>eps THEN addKnot[sa];
	    sa ← sa.next;
	    ENDLOOP;
	};

--pick one knot per peak.  If peak on first/last boundary, may get double knot
PickKnots: PROC = {
	sa: SampleHandle ← sList.samples.next;	--first point is duplicated on end
	sign: PROC[r: REAL] RETURNS[INTEGER] = {
	    IF r>0 THEN RETURN[1];
	    IF r<0 THEN RETURN[-1];
	    RETURN[0];
	    };
--sa is the left edge of a peak
	addKnot: PROC[sa: SampleHandle] RETURNS [SampleHandle]= {
	    rt: SampleHandle ← sa;
	    max: SampleHandle ← sa;
	    UNTIL rt=NIL OR ABS[rt.value]<eps DO
	        IF sign[rt.value]#sign[sa.value] THEN EXIT;
		IF ABS[rt.value]>ABS[max.value] THEN max ← rt;
		rt ← rt.next;
		ENDLOOP;
--rt is the right edge of a peak, max is max value
	    NewCPoint[max];
	    RETURN[rt];
	    };
	eps: REAL ← GetReal[];
	cps.arcLen ← sList.arcLen;
	UNTIL sa=NIL DO
	    IF ABS[sa.value]>eps THEN sa ← addKnot[sa] ELSE sa ← sa.next;
	    ENDLOOP;
	};
WriteKnots: PROC= {
	s: STRING ← [40];
	JaMFnsDefs.PopString[s];
	IF cps.samples = NIL THEN RETURN;
	CurveDefs.OpenFile[ConvertUnsafe.ToRope[s]];
	CurveDefs.WriteSamples[cps.samples];
	CurveDefs.CloseFile[];
	};

NewCPoint: PROC[sa: SampleHandle] = {
	s: SampleHandle ← NIL;
	new: SampleHandle ← NEW[Sample];
	new↑ ← [NIL,NIL,sa.xy,sa.s,sa.value];
	IF cps.lastSample=NIL THEN {cps.samples ← cps.lastSample ← new}
	ELSE
	  IF new.s<cps.samples.s THEN {	--head of the list
	      IF TooClose[sa, cps.samples] THEN RETURN;
	      new.next ← cps.samples;
	      cps.samples.prev ← new;
	      cps.samples ← new;
	    }
	  ELSE {
	      FOR s ← cps.samples, s.next UNTIL s=cps.lastSample OR s.next.s>new.s DO ENDLOOP;
	      IF s=cps.lastSample AND cps.lastSample.s<new.s THEN {	--tail of the list
	          IF TooClose[cps.lastSample,sa] THEN RETURN;
		  cps.lastSample.next ← new;
		  new.prev ← cps.lastSample;
		  cps.lastSample ← new;
		  }
	      ELSE {
	          IF TooClose[s,sa] THEN RETURN;
		  IF TooClose[sa,s.next] THEN RETURN;
		  new.next ← s.next;	--middle of the list
		  s.next ← new;
	          new.next.prev ← new;
	          new.prev ← s;
		  };
	      };
	nknots ← nknots+1;
	Mark[XFormSamplePt[new.xy]];
	};

TooClose: PROC[k0,k1: SampleHandle] RETURNS[BOOLEAN] = {
	sa0,sa1: SampleHandle;
	count: NAT ← 0;
	FOR sa0 ← sList.samples, sa0.next UNTIL sa0.s >= k0.s DO ENDLOOP;
	FOR sa1 ← sa0, sa1.next UNTIL sa1.s >= k1.s DO count ← count+1; ENDLOOP;
	IF count<4 THEN RETURN[TRUE] ELSE RETURN[FALSE];
	};
	
DeleteCPoint: PROC[sa: SampleHandle] = {
	s: SampleHandle ← NIL;
	IF cps.lastSample=NIL THEN RETURN;
	IF sa=cps.samples THEN {
	      cps.samples ← cps.samples.next;
	      IF cps.samples#NIL THEN cps.samples.prev ← NIL ELSE cps.lastSample←NIL;
	      }
	  ELSE {
	      FOR s ← cps.samples, s.next UNTIL s=NIL OR s=sa DO ENDLOOP;
	      IF s=NIL THEN RETURN;
	      s.prev.next ← s.next;
	      IF sa#cps.lastSample THEN sa.next.prev ← sa.prev ELSE cps.lastSample ← sa.prev;
	      };
	nknots ← nknots-1;
	EraseMark[XFormSamplePt[sa.xy]];
	};

InsertKnot: PROC = {
	s: SampleHandle ← NIL;
	tol: REAL ← 2/saScale;
	y: REAL ← GetReal[];
	x: REAL ← GetReal[];
	findS: PROC [sa: SampleHandle] = {
		IF x IN [sa.xy.x-tol..sa.xy.x+tol] AND y IN [sa.xy.y-tol..sa.xy.y+tol] THEN s ← sa;
		};
	x ← (x-saOff.x)/saScale;
	y ← (y-saOff.y)/saScale;
	ForAllSamples[findS];
	IF s#NIL THEN NewCPoint[s];
	};
	       
RemoveKnot: PROC = {
	tol: REAL ← 2/saScale;
	y: REAL ← GetReal[];
	x: REAL ← GetReal[];
	x ← (x-saOff.x)/saScale;
	y ← (y-saOff.y)/saScale;
	FOR sa: SampleHandle ← cps.samples, sa.next UNTIL sa=NIL DO
	    IF x IN [sa.xy.x-tol..sa.xy.x+tol] AND y IN [sa.xy.y-tol..sa.xy.y+tol] THEN {DeleteCPoint[sa]; EXIT};
	    ENDLOOP;
	};
	       
nknots: INTEGER ← 0;
cps: SListHandle ← NEW[SList];
lsHandle: LSFit.Handle;
InitLSFit: PROC = {ENABLE ABORTED => GOTO Quit;
	sa: SampleHandle;
	slist: SListHandle;
	scale,initlS: REAL ← 0;
	FOR sa ← sList.samples, sa.next UNTIL sa=NIL DO
	  IF cps.samples.s IN [sa.s..Next[sa].s) THEN EXIT;
	  ENDLOOP;
	IF sa=NIL THEN RETURN;
--sa now points to the break point.  Make a new slist that starts at sa
	initlS ← sa.s;	--will need this to translate the knots
	slist ← CopySamples[sList,sa];
	lsHandle ← LSFit.Create[SamplesToSequence[slist.samples]];
	lsHandle.closedCurve ← sList.closed;
	lsHandle.knots ← NEW[Seq.RealSequenceRec[nknots+1]];
	sa ← cps.samples;
	scale ← slist.arcLen/(IF cps.arcLen=0.0 THEN 1.0 ELSE cps.arcLen);
	FOR i: INTEGER IN [0..nknots) DO
	  new: REAL ← scale*sa.s-initlS;
	  IF new<0 THEN new ← slist.arcLen+new;
	  lsHandle.knots[i] ← new;
	  sa ← sa.next;
	  ENDLOOP;
	lsHandle.knots[0] ← 0;	--might want to flag this if it isn't already true
	lsHandle.knots[nknots] ← slist.arcLen;
	LSFit.InitialTValues[lsHandle];
	LSFit.FitXAndY[lsHandle];
        EXITS Quit => {}
	};

ShowKnots: PROC = {
    FOR cp: SampleHandle ← cps.samples, cp.next UNTIL cp=NIL DO
	FOR sa : SampleHandle ← sList.samples, sa.next UNTIL sa=NIL DO
	  IF cp.s IN [sa.s..Next[sa].s) THEN {Mark[XFormSamplePt[sa.xy]]; EXIT};
	  ENDLOOP;
        ENDLOOP
    };

Mark: PROC[pt: Vec] = {
    Paint: PROC[dc: Graphics.Context] = {
      Graphics.DrawBox[dc,[pt.x-2,pt.y-2,pt.x+2,pt.y+2]];
      };
    TJaMGraphics.Painter[Paint];
    };

EraseMark: PROC[pt: Vec] = {
    Paint: PROC[dc: Graphics.Context] = {
      Graphics.SetColor[dc,Graphics.white];
      Graphics.DrawBox[dc,[pt.x-2,pt.y-2,pt.x+2,pt.y+2]];
      Graphics.SetColor[dc,Graphics.black];
      };
    TJaMGraphics.Painter[Paint];
    };
    	
sList: SListHandle ← NIL;

AnalyzeContour: PUBLIC PROCEDURE[contour: ContourHandle] = {
	i: INTEGER ← 0;
--check all the null cases
	IF contour=NIL THEN RETURN ELSE sList ← CopySamples[contour.sLists];
	DoAnalysis[];
	};

CopySamples: PROC[from: SListHandle, first: SampleHandle ← NIL ] RETURNS[to: SListHandle] = {
	start,sa: SampleHandle;
	old: SListHandle ← sList;
	IF first=from.samples THEN first ← NIL;	--this is the convention
	IF first=NIL THEN start ← from.samples ELSE start ← first;
	StartSamples[start.xy.x,start.xy.y];
	FOR sa ← start.next, sa.next UNTIL sa=NIL DO
	    NewSample[sa.xy.x,sa.xy.y];
	    ENDLOOP;
	IF first#NIL THEN FOR sa ← from.samples, sa.next UNTIL sa.next=first DO
	    NewSample[sa.xy.x,sa.xy.y];
	    ENDLOOP;
	sList.closed ← from.closed;
	SetArcLen[];
	to ← sList;
	sList ← old;
	};

SetArcLen: PROC = {
	s: REAL ← 0;
	sList.nsamples ← 0;
	sList.samples.s ← 0;
	FOR sa: SampleHandle ← sList.samples.next, sa.next UNTIL sa=NIL DO
	    s ← s + Mag[Sub[sa.xy, sa.prev.xy]];
	    sa.s ← s;
	    sList.nsamples ← sList.nsamples+1;
	    ENDLOOP;
	IF sList.closed THEN s ← s + Mag[Sub[sList.samples.xy, sList.lastSample.xy]];
	sList.arcLen ← s;
	};

DoAnalysis: PROC = {
	IF sList=NIL THEN RETURN;
	IF sList.nsamples <= 3 THEN RETURN;
	ForAllSamples[sampleProc];
	SamplesToValues[];
	PushInteger[sList.nsamples];
	};
	
ForAllSamples: PROC[proc: SampleProc] = {
	FOR sa: SampleHandle ← sList.samples, sa.next UNTIL sa=NIL DO
	  proc[sa ! Real.RealException => {sa.value ← 0; CONTINUE}];
	  ENDLOOP;
	};

Dynfilter: PROCEDURE =
    BEGIN ENABLE ABORTED => GOTO Quit;
    tolerance: REAL = GetReal[];
    totalBadness: REAL ← 0;
    segments: DynFit.Segments;
    newZ: Seq.ComplexSequence;
    lsHandle: LSFit.Handle;
    len: NAT;
    IF sList=NIL THEN RETURN;
    lsHandle ← LSFit.Create[SamplesToSequence[sList.samples]];
    [segments,totalBadness] ← DynFit.FitSegments[lsHandle.z,tolerance];
    [newZ,len] ← DynFit.SampleSegments[segments,lsHandle.z];
    StartSamples[newZ[0].x,newZ[0].y];
    FOR i:NAT IN (0..len) DO NewSample[newZ[i].x,newZ[i].y] ENDLOOP;
    SetArcLen[];
    EXITS Quit => {}
    END;

StartSamples: PROC[x,y: REAL] = {
	sList ← NEW[SList];
	sList.samples ← NEW[Sample];
	sList.samples.xy ← [x,y];
	sList.lastSample ← sList.samples;
	};

sLen: REAL ← 1.0/256.0;
SetSLen: PROC = {sLen ← GetReal[]};

NewSample: PROC[x,y: REAL] = {OPEN Vector;
	new: SampleHandle;
	s: REAL ← Mag[Sub[[x,y],sList.lastSample.xy]];
	IF s>sLen THEN { 
		new ← NEW[Sample];
		sList.lastSample.next ← new;
		new.prev ← sList.lastSample;
		sList.lastSample ← sList.lastSample.next;
		new.xy ← [x,y];
		new.s ← s + new.prev.s;
		};
	};
			
FilterValues: PROC[val: Values] = {
	new,old: REAL ← 0;
	i: INTEGER ← 0;
	length: INTEGER ← val.l;
	closed: BOOLEAN ← sList.closed;
	IF val=NIL THEN RETURN;
	length ← val.l;
	IF closed THEN new ← (val[length-1].v+val[0].v+val[1].v)/3
	ELSE new ← (val[0].v+val[1].v)/2;
	FOR i IN [1..length-1) DO
		old ← new;
		new ← (val[i-1].v+val[i].v+val[i+1].v)/3;
		val[i-1].v ← old;
		ENDLOOP;
	old ← new;
	IF closed THEN val[length-1].v ← (val[length-2].v+val[length-1].v+val[0].v)/3
	ELSE val[length-1].v ← (val[length-2].v+val[length-1].v)/2;
	val[length-2].v ← old;
	};


WriteContour: PROC = {
	lastContour: ContourHandle ← NEW[Contour];
       	s: REAL;
	minX ← 10000; minY ← 10000;
	maxX ← 0; maxY ← 0;
	SetMinMax[];
	s ← maxY-minY;
	IF maxX-minX > s THEN s ← maxX-minX;
	lastContour.sLists ← sList;
	CurveDefs.WriteContour[lastContour,1/s,[-minY,-minX]]
	};
	
minX,maxX,minY,maxY: REAL ← 0;
SetMinMax: PROC [] = {
	compare: PROC [sa: SampleHandle] = {
		IF sa.xy.x>maxX THEN maxX ← sa.xy.x;
		IF sa.xy.y>maxY THEN maxY ← sa.xy.y;
		IF sa.xy.x<minX THEN minX ← sa.xy.x;
		IF sa.xy.y<minY THEN minY ← sa.xy.y;
		};
	ForAllSamples[compare];
	};

GetVS: PROC =
    BEGIN ENABLE ABORTED => GOTO Quit;
    i:INTEGER ← JaMFnsDefs.PopInteger[];
    v:Vec;
    v.x ← GetS[values,i];
    v.y ← GetValue[values,i];
    PushPoint[v];
    EXITS Quit => {}
    END;

NValues: PROC = {PushInteger[values.l]};		

PushPoint: PROC[p: Vec] = {JaMFnsDefs.PushReal[p.x]; JaMFnsDefs.PushReal[p.y]};

--the following routines are exported to AnalyzeDefs

PushKnots: PROC = {
	FOR sa: SampleHandle ← cps.samples, sa.next UNTIL sa=NIL DO
	    PushPoint[sa.xy];
		ENDLOOP;
    };

-- Initialization starts here
Register[".nvalues"L,NValues];
Register[".getvs"L,GetVS];

Register[".analyze"L,DoAnalysis];
Register[".baseline"L,BaseLine];
Register[".dcon"L,DrawContour];
Register[".drawlocs"L,RDrawLocs];
Register[".locatevs"L,RLocateVS];
Register[".locatexy"L,RLocateXY];
Register[".dfcon"L,Dynfilter];
Register[".pk"L,PickKnots];
Register[".psk"L,PickSumKnots];
Register[".insertknot"L,InsertKnot];
Register[".removeknot"L,RemoveKnot];
Register[".clearknots"L,ClearKnots];
Register[".startls"L,InitLSFit];

Register[".saoff"L,SAOff];
Register[".vsoff"L,VSOff];
Register[".sscale"L,SScale];
Register[".sascale"L,SaScale];
Register[".vscale"L,VScale];
Register[".anslen"L,SetSLen];
Register[".abs"L,SetAbs];

Register[".plotvalues"L,RPlotValues];
Register[".plotdvalues"L,RPlotDValues];
Register[".dvalues"L,DValues];
Register[".filtervalues"L,RFilterValues];
Register[".filterdvalues"L,FilterDValues];
Register[".valuestodvalues"L,ValuesToDValues];
Register[".dvaluestovalues"L,DValuesToValues];
Register[".samplestovalues"L,SamplesToValues];
Register[".valuestosamples"L,ValuesToSamples];
Register[".dx"L,SetDX];
Register[".dy"L,SetDY];
Register[".setmag"L,SetMag];
Register[".setcos"L,SetCos];
Register[".angle"L,SetAngle];
Register[".dxangle"L,SetDxAngle];
Register[".dyangle"L,SetDyAngle];
Register[".arclength"L,SetArcLength];
Register[".curvature"L,SetCurvature];
Register[".ccurve"L,SetCCurve];
Register[".ccurveangle"L,SetCCurveAngle];
Register[".andebug"L,Debug];
Register[".showknots"L,ShowKnots];
Register[".pushknots"L,PushKnots];


}.