-- CurveTracer.mesa
-- Michael Plass  August 3, 1982 2:41 pm

-- This provides an alternative method for getting a sampled curve from an AIS image, appropriate if the image consists of a line drawing.  The basic idea is to first blur the image so that the backround just reaches the center of lines of the desired width, and then to trace out the 'valley' thus defined.
-- CurveTracer is driven by JaM.  Here are the calls:
-- .ctraceopen	<string> => .	Reads the AIS file.
-- .ctraceclose => .		Closes the AIS file.
-- .csamp x y w => intensity.  Returns the intensity of the blurred image at the given point.
-- .ctrace	x1 y1 x2 y2 w t n => .  Traces a curve.  The expected line width is w.  The trace will start at (x1,y1), and will end somewhere.  One possible end condition is that the trace comes close to (x2,y2) or (x1,y1), another is that n points have been sampled. A third possibility is that the sampled value exceeds the threshold t. The trace starts out in the approximate direction of (x2,y2).

DIRECTORY
  AIS USING [Error, FRef, WRef, OpenFile, CloseFile, OpenWindow, ReadSample, CloseWindow],
  Complex,
  Curve,
  Inline USING [LowHalf],
  Real USING [RoundI],
  JaMFnsDefs,
  IO
  ;

CurveTracer: PROGRAM IMPORTS AIS, Complex, Curve, Inline, Real, JaMFnsDefs, IO =
BEGIN

-- This stuff is borrowed from ContourMain
ais: AIS.FRef _ NIL;
window: AIS.WRef _ NIL;

AISToScreen: PROC [x,y: REAL] RETURNS[nx,ny: REAL] ={
	nx _ x; ny _ window.lastScan-y;
	};

ScreenToAIS: PROC [x,y: REAL] RETURNS[nx,ny: REAL] ={
	nx _ x; ny _ window.lastScan-y;
	};

OpenAIS: PROCEDURE = {
	ENABLE AIS.Error =>{JaMFnsDefs.JaMExec["(AIS file error) .print"]; CONTINUE};
	s:STRING_[64];
 	JaMFnsDefs.PopString[s];
 	ais_AIS.OpenFile[s,FALSE];
	window _ AIS.OpenWindow[ais];
	pieceLine _ piecePixel _ -LAST[NAT]; -- invalidate cache
	};

Free: PROCEDURE = {
	ENABLE AIS.Error =>{JaMFnsDefs.JaMExec["(AIS file error) .print"]; CONTINUE};
	AIS.CloseWindow[window]; window _ NIL;
	AIS.CloseFile[ais];  ais _ NIL;
	};
-- End of borrowed code

NarrowToNat: PROC [n: INT] RETURNS [NAT] =
    BEGIN
    IF n<0 THEN n_0;
    IF n>LAST[NAT] THEN n_LAST[NAT];
    RETURN[Inline.LowHalf[n]];
    END;

-- Here is a cache of a part of the blurred image.
pieceSize: NAT = 100;
BlurredLine: TYPE = ARRAY [0..pieceSize) OF NAT;
blurredPiece:  ARRAY [0..pieceSize) OF REF BlurredLine;
pieceLine, piecePixel: INTEGER;
blurriness: NAT _ 0;
GetPiece: PROC [line,pixel: INT] =
	BEGIN
	FOR i:NAT IN [0..pieceSize) DO
		FOR j:NAT IN [0..pieceSize) DO
			blurredPiece[i][j] _ 8*
				(IF i+line IN [window.firstScan..window.lastScan]
				AND j+pixel IN [window.firstPixel..window.lastPixel]
				THEN AIS.ReadSample[window, Inline.LowHalf[line+i], Inline.LowHalf[pixel+j]]
				ELSE 255);
			ENDLOOP;
		ENDLOOP;
	pieceLine _ line;
	piecePixel _ pixel;
	END;

BlurPiece: PROC = 
	BEGIN
	ENABLE AIS.Error =>{JaMFnsDefs.JaMExec["(AIS error) .print"]; CONTINUE};
	prevline,line,nextline: BlurredLine;
	FOR j:NAT IN [0..pieceSize) DO prevline[j] _ blurredPiece[0][j] ENDLOOP;
	FOR i:NAT IN [1..pieceSize-1) DO
		FOR j:NAT IN [0..pieceSize) DO line[j] _ blurredPiece[i][j] ENDLOOP;
		FOR j:NAT IN [0..pieceSize) DO nextline[j] _ blurredPiece[i+1][j] ENDLOOP;
		FOR j:NAT IN [1..pieceSize-1) DO
			blurredPiece[i][j] _ (
				prevline[j-1] + prevline[j] + prevline[j+1] +
				    line[j-1] +     line[j] +     line[j+1] +
				nextline[j-1] + nextline[j] + nextline[j+1])/9;
			ENDLOOP;
		prevline _ line;
		ENDLOOP;
	END;

GetBlurredPiece: PROC[x,y: REAL, b: NAT] = 
	BEGIN -- x and y are in screen coordinates, b is the blurriness
	-- The piece is centered around (x,y)
	line,pixel: REAL;
	[pixel,line] _ ScreenToAIS[x,y];
	line _ line - pieceSize/2;
	pixel _ pixel - pieceSize/2;
	GetPiece [Real.RoundI[line], Real.RoundI[pixel]];
	THROUGH [1..b] DO
		BlurPiece;
		ENDLOOP; 
	blurriness _ b;
	END;

BlurredSample: PROC [a: Complex.Vec, b: NAT] RETURNS [intensity: REAL] = 
	BEGIN OPEN a; -- a is in screen coordinates, b is the blurriness
	line, pixel: REAL;
	[pixel,line] _ ScreenToAIS[x,y];
	IF b # blurriness
	OR NOT (     line IN [pieceLine+b..pieceSize+pieceLine-b)
	        AND pixel IN [piecePixel+b..pieceSize+piecePixel-b))
		THEN GetBlurredPiece[x,y,b];
	intensity _ blurredPiece[Real.RoundI[line]-pieceLine][Real.RoundI[pixel]-piecePixel]/8.0;
	IF FALSE THEN {IO.Put[Curve.TTY[],IO.real[a.x],IO.char[' ]];
		IO.Put[Curve.TTY[],IO.real[a.y],IO.char[' ],IO.real[intensity]];
		IO.Put[Curve.TTY[],IO.string["\n"]]};
	END;

Sample: PROC = {
	x,y,b: REAL;
	b _ MIN[pieceSize/2 - 1, ABS[JaMFnsDefs.GetReal[]]];
	y _ JaMFnsDefs.GetReal[];
	x _ JaMFnsDefs.GetReal[];
	JaMFnsDefs.PushReal[BlurredSample[[x,y],Real.RoundI[b]]];
	};

DistSqr: PROC[a,b: Complex.Vec] RETURNS [REAL] =
	{RETURN[Complex.SqrAbs[Complex.Sub[a,b]]]};

Direction: TYPE =  [0..8);
delta: ARRAY  Direction OF Complex.Vec = [[1,0],[1,1],[0,1],[-1,1],[-1,0],[-1,-1],[0,-1],[1,-1]];

AwayFrom: PROC[i,j: INTEGER] RETURNS [BOOLEAN] =
	BEGIN
	i _ i-j;
	WHILE i<-2 DO i_i+8 ENDLOOP;
	WHILE i>5 DO i_i-8 ENDLOOP;
	RETURN[i>2];
	END;

CurveTrace: PROC = {
	x1,y1,x2,y2,w,t: REAL;
	b, n: NAT;
	closestApproach: REAL _ 10.0E+20;
	stopper: Complex.Vec;
	CloseEnough: PROC [a: Complex.Vec] RETURNS [BOOLEAN] = 
		BEGIN OPEN a;
		IF MAX[ABS[a.x-stopper.x],ABS[a.y-stopper.y]] <= .001 THEN RETURN[TRUE];
		IF x IN [x2-w/2..x2+w/2] AND y IN [y2-w/2..y2+w/2] THEN
			BEGIN
			dSqr: REAL _ DistSqr[a,[x2,y2]];
			IF dSqr>closestApproach THEN RETURN[TRUE];
			closestApproach _ dSqr;
			END;
		RETURN[FALSE];
		END;

	current: Complex.Vec;
	cameFrom: Direction _ FIRST[Direction];

	n _ JaMFnsDefs.PopInteger[];
	t _ JaMFnsDefs.GetReal[];
	w _ ABS[JaMFnsDefs.GetReal[]];
	b _ Real.RoundI[MIN[pieceSize/2.0 - 1, w]];
	y2 _ JaMFnsDefs.GetReal[];
	x2 _ JaMFnsDefs.GetReal[];
	y1 _ JaMFnsDefs.GetReal[];
	x1 _ JaMFnsDefs.GetReal[];
	current _ [x1,y1];
	stopper _ [x2,y2];
	FOR d: Direction IN Direction DO 
		IF DistSqr[[x2,y2],Complex.Add[current,delta[d]]] >
			DistSqr[[x2,y2],Complex.Add[current,delta[cameFrom]]] THEN
			cameFrom _ d
		ENDLOOP; 
	FOR i:NAT IN [1..n] UNTIL CloseEnough[current] DO
		high: REAL _ 0;
		low: REAL _ 1000000;
		new: Complex.Vec;
		newd: Direction;
		IF i>=b THEN Curve.AddSample[Curve.defaultHandle, current.x, current.y];
		IF i=b THEN stopper _ current;
		FOR d: Direction IN Direction DO 
			IF AwayFrom[cameFrom,d] THEN
				BEGIN
				p: Complex.Vec _ Complex.Add[current,delta[d]];
				value: REAL _ BlurredSample[p,b];
				IF value < low THEN {new _ p; low _ value; newd _ d};
				IF value > high THEN high _ value;
				END;
			ENDLOOP;
		cameFrom _ IF newd<4 THEN newd+4 ELSE newd-4;
		IF high-low < .1 THEN EXIT;
		IF low > t THEN EXIT;
		IF JaMFnsDefs.GetJaMBreak[] THEN EXIT;
		current _ new;
		ENDLOOP; 
	};

FOR i:NAT IN [0..pieceSize) DO
	blurredPiece[i] _ NEW[BlurredLine];
	ENDLOOP;

JaMFnsDefs.Register[".ctraceopen"L,OpenAIS];
JaMFnsDefs.Register[".ctraceclose"L,Free];
JaMFnsDefs.Register[".csamp"L,Sample];
JaMFnsDefs.Register[".ctrace"L,CurveTrace];



END.