-- 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];


}.