-- LSSplineImpl.mesa
--last edited by Michael Plass  August 4, 1982 10:19 am

DIRECTORY Complex, Cubic, LSPiece, Seq, Vector, LSSpline;

LSSplineImpl: PROGRAM IMPORTS LSPiece EXPORTS LSSpline =
BEGIN

Handle: TYPE = LSSpline.Handle;
Rec: TYPE = LSSpline.Rec; 

Create: PUBLIC PROCEDURE [lengthHint: NAT _ 10] RETURNS [handle: Handle] = {
	handle _ NEW[Rec];
	handle.z _ NEW[Seq.ComplexSequenceRec[lengthHint]];
	handle.tan _ NEW[Seq.ComplexSequenceRec[lengthHint]];
	handle.t _ NEW[Seq.RealSequenceRec[lengthHint]];
	handle.length _ 0;
	};

Reset: PUBLIC PROCEDURE [handle: Handle] = {handle.length _ 0};

Sample: PUBLIC PROCEDURE [handle: Handle, samplePoint, tangent: Complex.Vec _ [0,0]] = {
	WHILE handle.length >= handle.z.length DO
		temp: Handle _ Create[handle.length + handle.length/2 + 1];
		FOR i:NAT IN [0..handle.length) DO
			Sample[temp, handle.z[i], handle.tan[i]];
			ENDLOOP;
		handle^ _ temp^;
		ENDLOOP;
	handle.z[handle.length] _ samplePoint;
	handle.tan[handle.length] _ tangent;
	handle.length _ handle.length + 1;
	}; 

CubicProc: TYPE = PROCEDURE [Cubic.Bezier];

Subdivide: PUBLIC PROCEDURE [handle: Handle, cubicProc: CubicProc, tolerance: REAL, from: NAT _ 0, thru: NAT _ LAST[NAT]] = {
	b: Cubic.Bezier;
	err: REAL;
	iterations: NAT;
	maxDev: REAL;
	IF handle.length = 0 THEN RETURN;
	thru _ MIN[thru, handle.length - 1];
	IF from >= thru THEN RETURN;
	IF thru = from+1 THEN {
		b0: Vector.Vec _ handle.z[from];
		b3: Vector.Vec _ handle.z[thru];
		cubicProc[[b0, b0, b3, b3]];
		RETURN;
		};
	[b, err, iterations, maxDev] _ LSPiece.FitPiece[
		z: handle.z,
		t: handle.t,
		from: from,
		thru: thru,
		maxd: tolerance,
		maxit: firstTryIterationLimit,
		initTangent: handle.tan[from],
		finalTangent: handle.tan[thru]
		];
	stats.firstTryCount _ stats.firstTryCount + 1;
	stats.firstTryIterations _ stats.firstTryIterations + iterations;
	IF maxDev <= tolerance*firstTryToleranceFactor THEN {
		[b, err, iterations, maxDev] _ LSPiece.FitPiece[
			z: handle.z,
			t: handle.t,
			from: from,
			thru: thru,
			maxd: tolerance/secondTryOptimistFactor,
			maxit: secondTryIterationLimit,
			initTangent: handle.tan[from],
			finalTangent: handle.tan[thru],
			useOldTValues: TRUE
			];
		stats.secondTryCount _ stats.secondTryCount + 1;
		stats.secondTryIterations _ stats.secondTryIterations + iterations;
		};
	IF maxDev <= tolerance THEN cubicProc[b]
	ELSE {
		mid: NAT _ (from+thru)/2;
		Subdivide[handle, cubicProc, tolerance, from, mid];
		Subdivide[handle, cubicProc, tolerance, mid, thru];
		}; 
	};

firstTryIterationLimit: NAT _ 3;
firstTryToleranceFactor: REAL _ 3.0;
secondTryIterationLimit: NAT _ 20;
secondTryOptimistFactor: REAL _ 2.0;

stats: RECORD [
	firstTryCount: INT _ 0,
	firstTryIterations: INT _ 0,
	secondTryCount: INT _ 0,
	secondTryIterations: INT _ 0
	];

END.