-- GrowPatches.mesa
--last edited by Maureen Stone November 7, 1983 3:02 pm
--last edited by Michael Plass  August 30, 1982 1:42 pm

DIRECTORY Complex, Cubic, Curve, JaMFnsDefs, LSPiece, Seq, Vector, Real, RealFns,  FitDisplayUtilities;

GrowPatches: PROGRAM IMPORTS Cubic, Curve, JaMFnsDefs, LSPiece, Vector, Real, RealFns, FitDisplayUtilities =
BEGIN OPEN JaMFnsDefs;

myEps: REAL _ .0001;	--something very small.  We're using maximum deviation as test
MyEps: PROC = {myEps _ GetReal[]};

FitPatches: PROC = {
	z: Seq.ComplexSequence _ Curve.CurrentSamples[Curve.defaultHandle];
	n: NAT _ z.length;
	t: Seq.RealSequence _ NEW[Seq.RealSequenceRec[n]];
	node,cusp: Seq.NatSequence;
	tangent,ctangents: Seq.ComplexSequence;
	maxit: NAT;
	b: Cubic.Bezier;
	err,maxd: REAL;
	iterations: NAT;
	nextCusp,cIndex: INT _ 0;
	[node,tangent] _ Curve.CurrentNodes[Curve.defaultHandle];
	[cusp,ctangents] _ Curve.CurrentCusps[Curve.defaultHandle];
	maxit _ PopInteger[];
	maxd _ GetReal[];
	Curve.ResetLinks[Curve.defaultHandle];
	nextCusp _ IF cusp.length=0 THEN -1 ELSE cusp[0];
	FOR i:NAT IN [1..node.length) DO
		t1,t2: Complex.Vec;
		IF node[i-1]=nextCusp THEN {
			t1 _ ctangents[2*cIndex+1];
			cIndex _ cIndex+1;
			IF cIndex<cusp.length THEN nextCusp _ cusp[cIndex];  --stop on end cusp
			}
		ELSE t1 _ tangent[i-1];
		IF node[i]=nextCusp THEN t2 _ ctangents[2*cIndex] 
		ELSE t2 _ tangent[i];
		IF GetJaMBreak[] THEN EXIT;
		[b,err,iterations] _
			FitPiece[z: z, t: t,
			from: node[i-1], thru: node[i],
			eps: myEps, maxd: maxd, maxit: maxit,
			initFree: FALSE, finalFree: FALSE,
			initTangent: t1, finalTangent: t2];
		Curve.AddLink[Curve.defaultHandle, b];
		ENDLOOP;
	};
	
GrowPatches: PROC = {
	pstart,pend: NAT;
	z: Seq.ComplexSequence _ Curve.CurrentSamples[Curve.defaultHandle];
	n: NAT _ z.length;
	t: Seq.RealSequence _ NEW[Seq.RealSequenceRec[n]];
	node,cusp: Seq.NatSequence;
	tangent,ctangents: Seq.ComplexSequence;
	maxit: NAT;
	b,prevB: Cubic.Bezier;
	err,maxd: REAL;
	iterations,npieces: NAT;
	maxdev: REAL _ 0;
	nextCusp,cIndex: INT _ 0;
	t1,t2: Complex.Vec;
	[node,tangent] _ Curve.CurrentNodes[Curve.defaultHandle];
	[cusp,ctangents] _ Curve.CurrentCusps[Curve.defaultHandle];
	maxit _ PopInteger[];
	maxd _ GetReal[];
	pstart _ npieces _ 0;
	Curve.ResetLinks[Curve.defaultHandle];
	nextCusp _ IF cusp.length=0 THEN -1 ELSE cusp[0];
	WHILE pstart < node.length-1 DO
		IF GetJaMBreak[] THEN EXIT;
		IF pstart=nextCusp THEN {
			t1 _ ctangents[2*cIndex+1];
			cIndex _ cIndex+1;
			IF cIndex<cusp.length THEN nextCusp _ cusp[cIndex];  --stop on end cusp
			}
		ELSE t1 _ tangent[pstart];
		pend _ pstart+1;
		WHILE pend < node.length DO
			IF pend=nextCusp THEN t2 _ ctangents[2*cIndex]
			ELSE t2 _ tangent[pend];
			[b,err,iterations,maxdev] _ FitPiece[z: z, t: t,
					from: node[pstart], thru: node[pend],
					eps: myEps, maxd: maxd, maxit: maxit,
					initFree: FALSE, finalFree: FALSE,
					initTangent: t1, finalTangent: t2];
			IF pend=nextCusp THEN EXIT;
			IF maxdev>maxd THEN {
				IF pend-pstart>1 THEN {pend _ pend - 1; b _ prevB};
				EXIT};
			prevB _ b;
			pend _ pend+1;
			ENDLOOP;
		npieces _ npieces+1;
		Curve.AddLink[Curve.defaultHandle, b];
		pstart _ pend;
		ENDLOOP;
	PushInteger[npieces];
	};

FitPiece: PROCEDURE
	[z: Seq.ComplexSequence, t: Seq.RealSequence _ NIL,
	from, thru: NAT, -- fit points in range [from..thru] modulo z.length
	eps: REAL _ .00001, -- stop iterating if the sum of squares gets below this amount
    maxd: REAL _ .00001,  --stop iterating if the maximum deviation gets below this amount
	maxit: NAT _ 500, -- limit on number of iterations
	deltat: REAL _ 0.000005, -- stop when the t values each change by less than this
	initFree, finalFree: BOOLEAN _ FALSE,
	initTangent, finalTangent: Complex.Vec _ [0,0],
	useOldTValues: BOOLEAN _ FALSE,
	unitDeviation: REAL _ 1,
	exponent: NAT _ 2
	]
	RETURNS [b: Cubic.Bezier, err: REAL, iterations: NAT, maxDev: REAL, powerError: REAL] = {

--this makes sure there are enough points to give a good fit, then calls LSPiece.FitPiece
	npoints: NAT _ IF thru>from THEN thru-from+1 ELSE thru-from+z.length+1;
	minPoints: NAT _ 6;
	IF npoints>=minPoints THEN {
		[b,err,iterations, maxDev] _ LSPiece.FitPiece[z: z, t: t, from: from, thru: thru,
			eps: eps, maxd: maxd, maxit: maxit, deltat: deltat,
			initFree: initFree, finalFree: finalFree,
			initTangent: initTangent,
			finalTangent: finalTangent, useOldTValues: useOldTValues];
		powerError _ Error[b, z, t, from, thru, unitDeviation, exponent];
		}
	ELSE {	--else interpolate and then call
		newZ: Seq.ComplexSequence _ NEW[Seq.ComplexSequenceRec[minPoints]];
		newT: Seq.RealSequence  _ NEW[Seq.RealSequenceRec[minPoints]];
		ninterp: INT _ minPoints-npoints;
		zi: NAT _ 0;
		modZ: PROC[index: NAT] RETURNS [NAT] = {RETURN[index MOD z.length]};
		interpolate: PROC[p0,p1: Vector.Vec, npts: NAT] = {
			dv: Vector.Vec _ Vector.Mul[Vector.Sub[p1,p0],1.0/(npts+1)]; 
			FOR j: NAT IN [0..npts) DO
				newZ[zi] _ Vector.Add[newZ[zi-1],dv];
				zi _ zi+1;
				ENDLOOP;
			};
		newZ[0] _ z[from];
		IF npoints=2 THEN {zi _ 1; interpolate[p0: z[from], p1: z[thru], npts: ninterp]}
		ELSE {
			d: Seq.RealSequence _ NEW[Seq.RealSequenceRec[npoints-1]];
			minD: REAL _ 100000000;
			maxD: REAL _ 0;
			totalD: REAL _ 0;
			np: NAT _ 0;
			i,n: NAT _ 0;
			adjust: BOOLEAN _ FALSE;
			FOR i _ from, modZ[i+1] UNTIL i=thru DO
				d[n] _ Vector.Mag[Vector.Sub[z[i],z[modZ[i+1]]]];
				totalD _ totalD+d[n];
				IF d[n]<minD THEN minD _ d[n];
				IF d[n]>maxD THEN maxD _ d[n];
				n _ n+1;
				ENDLOOP;
			FOR i IN [0..d.length) DO np _ np+Real.RoundI[(d[i]/totalD)*ninterp]; ENDLOOP;
			IF np#ninterp THEN adjust _ TRUE;
			FOR i IN [0..d.length) DO
				p: INTEGER _ Real.RoundI[(d[i]/totalD)*ninterp];
				newZ[zi] _ z[modZ[from+i]]; zi _ zi+1;
				IF adjust THEN{
					IF np>ninterp AND d[i]=minD THEN {p_p-1; np _ np-1};
					IF np<ninterp AND d[i]=maxD THEN {p_p+1; np _ np+1};
					IF np=ninterp THEN adjust _ FALSE};	--should only be off by one 
				IF p>0 THEN interpolate[p0: z[modZ[from+i]], p1: z[modZ[from+i+1]], npts: p];
				ENDLOOP;
			};
		newZ[newZ.length-1] _ z[thru];
		[b,err,iterations, maxDev] _ LSPiece.FitPiece[z: newZ, t: newT, from: 0, thru: newZ.length-1,
			eps: eps, maxd: maxd, maxit: maxit, deltat: deltat,
			initFree: initFree, finalFree: finalFree,
			initTangent: initTangent,
			finalTangent: finalTangent, useOldTValues: useOldTValues];
		powerError _ Error[b, newZ, newT, 0, newZ.length-1, unitDeviation, exponent];
		};
	RETURN[b,err,iterations, maxDev, powerError];
	};

ChopPatches: PROC = {
	pstart,pend,lastbad,next: NAT;
	z: Seq.ComplexSequence _ Curve.CurrentSamples[Curve.defaultHandle];
	n: NAT _ z.length;
	t: Seq.RealSequence _ NEW[Seq.RealSequenceRec[n]];
	node: Seq.NatSequence;
	tangent: Seq.ComplexSequence;
	maxit: NAT;
	b: Cubic.Bezier;
	err,maxd,maxdev: REAL;
	iterations,npieces: NAT;
	Curve.DeleteNode[Curve.defaultHandle, n-2];
	Curve.DeleteNode[Curve.defaultHandle, n-3];
	[node,tangent] _ Curve.CurrentNodes[Curve.defaultHandle];
	maxit _ PopInteger[];
	maxd _ GetReal[];
	pstart _ npieces _ 0;
	Curve.ResetLinks[Curve.defaultHandle];
	WHILE pstart < node.length-1 DO
	  pend _ lastbad _ node.length-1;
	  IF GetJaMBreak[] THEN EXIT;
	  DO
		IF node[pend]-node[pstart]<3 THEN EXIT;
		[b,err,iterations,maxdev] _ FitPiece[z: z, t: t, from: node[pstart], thru: node[pend], eps: myEps, maxd: maxd,maxit: maxit, initFree: FALSE, finalFree: FALSE, initTangent: tangent[pstart], finalTangent: tangent[pend]];
		IF maxdev<=maxd
		  THEN {
		   next _ (pend+lastbad)/2;
			IF next=lastbad OR next=pend
			 THEN EXIT  --this is best solution as the next bigger one has already been marked bad
			 ELSE pend _ next; --try a bigger one
		   }
		   ELSE {lastbad _ pend; pend _ (pend+pstart)/2};  --subdivide and try a smaller one.
		ENDLOOP;
	  WHILE node[pend]-node[pstart] < 3 DO pend _ pend + 1 ENDLOOP;
	  [b,err,iterations,maxdev] _ FitPiece[z: z, t: t, from: node[pstart], thru: node[pend], eps: myEps, maxd: maxd, maxit: maxit, initFree: FALSE, finalFree: FALSE, initTangent: tangent[pstart], finalTangent: tangent[pend]];
	  npieces _ npieces+1;
	  Curve.AddLink[Curve.defaultHandle, b];
	  pstart _ pend;
	  ENDLOOP;
	PushInteger[npieces];
	};

useMaxD: BOOLEAN _ TRUE;
DynSpline: PROC =
	BEGIN
	z: Seq.ComplexSequence _ Curve.CurrentSamples[Curve.defaultHandle];
	t: Seq.RealSequence _ NEW[Seq.RealSequenceRec[z.length]];
	maxit, exponent: NAT;
	bezier: Cubic.Bezier;
	unitDeviation: REAL;
	penalty: REAL;
	bestPrev: Seq.NatSequence;
	segCount: Seq.NatSequence;
	subtotalCost: Seq.RealSequence;
	OutputPatchesUpto: PROC [m: NAT] =
		BEGIN
		IF m>0 THEN
			BEGIN
			OutputPatchesUpto[bestPrev[m]];
			[bezier] _ FitPiece[
				z: z, t: t,
				from: node[bestPrev[m]], thru: node[m],
				eps: myEps,
				maxd: unitDeviation/4,
				maxit: maxit*2,
				initFree: FALSE, finalFree: FALSE,
				initTangent: tangent[bestPrev[m]],
				finalTangent: tangent[m]
				];
			Curve.AddLink[Curve.defaultHandle, bezier];
			END;
		END;
	n: NAT;
	node: Seq.NatSequence;
	tangent: Seq.ComplexSequence;
	oldbezier: Cubic.Bezier _ [[-1,-1],[-1,-1],[-1,-1],[-1,-1]];
	[node,tangent] _ Curve.CurrentNodes[Curve.defaultHandle];
	n _ node.length;
	exponent _ PopInteger[]; 
	maxit _ PopInteger[];
	unitDeviation _ GetReal[];
	penalty _ GetReal[];
	bestPrev _ NEW[Seq.NatSequenceRec[n]];
	segCount _ NEW[Seq.NatSequenceRec[n]];
	subtotalCost _ NEW[Seq.RealSequenceRec[n]];
	bestPrev[0] _ 0;
	segCount[0] _ 0;
	subtotalCost[0] _ 0.0;
	FOR m:NAT IN [1..n) DO
		badness: REAL;
		lowbad:REAL _ LAST[INT]/2;
		bestl:NAT _ 0;
		err, err2, thisErr:REAL _ 0;
		FitDisplayUtilities.HiLight[z[node[m]].x, z[node[m]].y];
		FOR l:NAT DECREASING IN [0..m-1] WHILE thisErr+penalty<lowbad DO
			IF GetJaMBreak[] THEN RETURN;
			IF node[m]-node[l] >= 1 THEN
				BEGIN
				FitDisplayUtilities.HiLight[z[node[l]].x, z[node[l]].y];
				[b: bezier, err: err2, powerError: err] _
				FitPiece[z: z, t: t,
						from: node[l], thru: node[m],
						eps: myEps,
						maxd: unitDeviation/2,
						maxit: maxit,
						initFree: FALSE, finalFree: FALSE,
						initTangent: tangent[l],
						finalTangent: tangent[m],
						unitDeviation: unitDeviation,
						exponent: exponent
						];
				thisErr _ IF useMaxD THEN err2 ELSE err;
				badness _ thisErr+subtotalCost[l] + penalty;
					IF badness < lowbad THEN
					BEGIN
					lowbad _ badness;
					bestl _ l;
					END;
				FitDisplayUtilities.HiLight[z[node[l]].x, z[node[l]].y];
				FitDisplayUtilities.ShowBezierInverted[oldbezier];
				FitDisplayUtilities.ShowBezierInverted[bezier];
				oldbezier _ bezier;
				END;
			ENDLOOP;
		subtotalCost[m] _ lowbad;
		segCount[m] _ segCount[bestl] + 1;
		bestPrev[m] _ bestl;
		FitDisplayUtilities.HiLight[z[node[m]].x, z[node[m]].y];
		ENDLOOP;
	FitDisplayUtilities.ShowBezierInverted[oldbezier];
	Curve.ResetLinks[Curve.defaultHandle];
	OutputPatchesUpto[n-1];
	END;		

Error: PROC [
	bezier: Cubic.Bezier,
	z: Seq.ComplexSequence,
	t: Seq.RealSequence,
	from, thru: NAT,
	unitDeviation: REAL,
	exponent: NAT
	] RETURNS [error: REAL] = {-- finds sum of powers of distances
	Wrap: PROC[i: NAT] RETURNS [imod: NAT] = INLINE
		{imod _ IF i>=z.length THEN i-z.length ELSE i};
	c: Cubic.Coeffs _ Cubic.BezierToCoeffs[bezier];
	unitDevSqrRecip: REAL _ 1/(unitDeviation*unitDeviation);
	error _ 0;
	FOR i: NAT _ from, Wrap[i+1] DO
		ti: REAL _ t[i];
		x: REAL _ ((c.c3.x*ti+c.c2.x)*ti+c.c1.x)*ti+c.c0.x;
		y: REAL _ ((c.c3.y*ti+c.c2.y)*ti+c.c1.y)*ti+c.c0.y;
		zi: Complex.Vec _ z[i];
		dSqr: REAL _ ((zi.x-x)*(zi.x-x) + (zi.y-y)*(zi.y-y))*unitDevSqrRecip;
		dPower: REAL _ 1;
		e: NAT _ exponent;
		WHILE e>1 DO dPower _ dPower*dSqr; e _ e-2 ENDLOOP;
		IF e=1 THEN dPower _ dPower*RealFns.SqRt[ABS[dSqr]];
		error _ error + dPower;
		IF i=thru THEN EXIT
		ENDLOOP;
	};

JaMFnsDefs.Register[".fit"L,FitPatches];		-- maxdev maxit => . Makes links between every pair of nodes
JaMFnsDefs.Register[".grow"L,GrowPatches];		-- maxdev maxit => . Makes links
JaMFnsDefs.Register[".chop"L,ChopPatches];		-- maxdev maxit => . Makes links
JaMFnsDefs.Register[".dynspline"L,DynSpline];		-- penalty unitdev maxit exponent => . Makes links using dynamic programming
JaMFnsDefs.Register[".myeps"L,MyEps];		--  myEps is sum of the distance squared

END.

Michael Plass  August 19, 1982 9:04 am. Added QuickTangents.
Michael Plass  August 19, 1982 2:15 pm. Enhanced dynspline error calculations.
Michael Plass  August 26, 1982 12:56 pm. Moved OpenCurve, DynNodes, and QuickTangents into NodeImpl.
Michael Plass  August 30, 1982 1:41 pm. Moved SetTangents to NodeImpl.
�Ê;–
"Mesa" style˜IprocšùÏc„
œ
Ïk	œmžœ
žœQžœ
žœžœ
�?ÏnœžœŸ
œžœNžœ$žœhžœžœžœžœÜ
žœžœžœžœžœ
žœžœžœžœ;žœžœ�œžœžœžœžœžœžœ
žœŠ
žœ
žœTžœŸœžœžœKžœ$žœhžœ%žœžœžœžœˆžœžœžœžœžœžœžœ
žœžœžœ;žœžœ�œžœ,žœžœžœžœžœ®
žœ
žœ.žœžœ
žœžœ
žœžœžœ#žœ&žœUžœŸœž	œ1žœžœ�4œžœ�?œ
žœ�Aœžœ�!œ	žœ
�8œžœžœCžœžœžœžœ	žœžœžœ
žœžœ�Zœ
žœžœžœ
žœ#žœžœžœßžœ�!œžœAžœ-žœžœžœžœžœžœžœžœžœžœJžœžœ
žœžœ>žœžœžœBžœžœ+žœžœžœ
žœžœžœžœžœžœžœUžœžœžœžœžœžœžœžœ-žœžœžœ
žœžœžœžœžœVžœžœžœžœžœžœžœžœžœžœ
žœ�œžœžœGžœžžœ-Ÿœžœ žœKžœ$žœYžœ&žœžœ�žœžœ'žœžœ
žœžœžœžœ
žœŠ
žœ
žœ@žœžœ$žœžœžœ
žœ�Kœžœ�œžœ,�#œžœžœžœžœŒ
žœ
žœ•
žœ'žœžœŸ	œžœžœ]žœ3žœ)žœžœ\ŸœžœžœžœžœžœžœÅ
žœ
žœžœžœžœ´žœ%žœ)žœWžœžœ
žœžœžœžœžœ
žœ
žœžœCžœžœ
ž
œ
žœ
žœžœžœžœ
žœžœžœžœü
žœ
žœ•
žœ	žœžœ:žœžœžœ-žœÆ
žœžœ—
žœxžœŸœžœVžœžœ
žœžœ	žœ�$œ
Ÿœžœžœžœžœžœžœ
žœžœHžœ1žœžœžœžœžœ2žœOžœFžœžœžœžœžœžœžœžœ%žœžœ
žœžœ1�=œ,�!œ,�!œ/�Mœ'�)œ
žœº˜¾b—…—1@6�