{Global declarations}
var TimeUnit,
				ScaleUnit,
				xLabel,yLabel,
				Equation     : string;
				TimeInt,
				i,
				x,y,
				Time0,
				xOrigin,
				yOrigin,
				wdt,hgt
				marginL,marginR,			{left and right margins for graphs}
				marginT,marginB,			{top and bottom margins for graphs}
				frameL,frameR,					{left and right coordinates for frame edge}
				frameT,frameB,					{top and bottom coordinates for frame edge}
				frameW,frameH,					{width and height of graph frame}
				nGridLines,								{number of grid lines on graph - approximate}
				txOffset,tyOffset,	{text coordinate offsets}
				cRad,
				FitType,
				Width,Height,
				theStack   : integer;
				Scale,
				xMin,yMin,
				xMax,yMax,
				tMin,tMax,
				a,b,c       : real;
    EverySlice:boolean;

procedure Error(s:string);
begin
	PutMessage(s);
	exit;
end; {Error}

{=================================================================}
procedure CurveFit(FitType:integer, ncPoints:integer);

var a11,a12,a13,a21,a22,a23,a31,a32,a33,
b1,b2,b3,c1,c2,c3,sumX,sumXX,sumXXX,sumXXXX,
sumY,sumYY,sumXY,sumXXY,x,y,DetOfSys,rfit : real;

begin
{}ShowMessage('Fitting data to curve...');
{initialize all of the variables}
	sumX:=0;
	sumXX:=0;
	sumXXX:=0;
	sumXXXX:=0;
	sumY:=0;
	sumYY:=0;
	sumXY:=0;
	sumXXY:=0;
	for i:=1 to ncPoints do begin
		x:=rUser1[i]; y:=rUser2[i];
		sumX:=sumX+x;
		sumXX:=sumXX+sqr(x);
		sumXXX:=sumXXX+x*sqr(x);
		sumXXXX:=sumXXXX+sqr(x)*sqr(x);
		sumY:=sumY+y;
		sumYY:=sumYY+sqr(y);
		sumXY:=sumXY+x*y;
		sumXXY:=sumXXY+sqr(x)*y;
	end; {for}

{build coefficient matrices}
	a11:=ncPoints;
	a12:=sumX;
	a21:=a12;
	a22:=sumXX;
	if (FitType=0) then begin
		a13:=0;
		a23:=0;
		a33:=1;
	end; {if}
	if (FitType=1) then begin
		a13:=a22;
		a23:=sumXXX;
		a33:=sumXXXX;
	end; {if}
	a31:=a13;
	a32:=a23;
	c1:=sumY;
	c2:=sumXY;
	c3:=sumXXY;

{Calculate the determinant of the system}
	DetOfSys:=
a11*(a22*a33-a23*a32)-a12*(a21*a33-a23*a31)+a13*(a21*a32-a22*a31);
	if DetOfSys=0 then begin
		PutMessage('Whoa, something is wrong here!');
		exit;
	end; {if}

{calculate coefficients}
	b1:=
(c1*(a22*a33-a23*a32)-a12*(c2*a33-a23*c3)+a13*(c2*a32-a22*c3))/DetOfSys;
	b2:=
(a11*(c2*a33-a23*c3)-c1*(a21*a33-a23*a31)+a13*(a21*c3-c2*a31))/DetOfSys;	
	if FitType='Line' then b3:=0 else 
		b3:=
(a11*(a22*c3-c2*a32)-a12*(a21*c3-c2*a31)+c1*(a21*a32-a22*a31))/DetOfSys;

{calculate correlation factor}
	rfit:=
(ncPoints*sumXY-sumX*sumY)/sqrt((ncPoints*sumXX-sqr(sumX))*(ncPoints*sumYY-sqr(sumY)));
	a:=b1; b:=b2; c:=b3;
	{}ShowMessage('');
end; {procedure CurveFit}

{=================================================================}
procedure MakeGraph(gWidth,gHeight,nPoints:integer;
																				gxMin,gxMax,gyMin,gyMax:real;
																			 gTitle:string);
{
Build a graph. Expects to find data in user arrays.
globals: nGridLines
									frame*, margin*
									x,y
									txOffset,tyOffset
									xLabel,yLabel
procedures: NewSize
												DrawPoint
}
var gxOrigin,gyOrigin,			{coordinates of graph origin}
				i,j,ndx,xp,yp
						:integer;
				xSpacing,ySpacing,			{spacing between grid lines}
				n, dx	:real;
				DoLine, done
						:boolean;
begin
 {¥begin open plot window}
	SaveState; SetBackgroundColor(0); SetForegroundColor(255);
	SetNewSize(gHeight,gWidth);
	MakeNewWindow(gTitle);
	SetText('Center Justified, No Background');
	SetFontSize(10);
	MoveTo(gHeight/2,tyOffset); write(yLabel);
	RotateLeft(true);
	MoveTo(gWidth/2,gHeight-4*tyOffset); write(xLabel);
	MoveTo(gWidth/2,gHeight-2*tyOffset); write(Equation);
	RestoreState;

	{¥Calculate frame size and location}
	frameL:=marginL; frameR:=gWidth-marginR;
	frameT:=marginT; frameB:=gHeight-marginB;
	frameW:=frameR-frameL;
	frameH:=frameB-frameT;

	{¥calculate coordinates of origin}
	gxOrigin:=round((0-gxMin)/(gxMax-gxMin)*frameW+frameL);
	gyOrigin:=round(frameB-(0-gyMin)/(gyMax-gyMin)*frameH);

	{¥calculate spacing between grid lines}
	xSpacing:={round}(frameW/nGridLines);
	ySpacing:={round}(frameH/nGridLines);

	{¥draw the grid using the background color}
	SaveState;
	SetForegroundColor(255);
	SetBackgroundColor(250);
	SetText('Center Justified, No Background');
	SetFontSize(9);
	y:=frameB+tyOffset;
	x:=gxOrigin;
	while (x<=frameR) do begin
		if (x>=frameL) then begin
			MakeLineRoi(x,frameT,x,frameB); Clear; KillRoi;
			n:=(x-frameL)*(gxMax-gxMin)/frameW+gxMin;
			MoveTo(x,y); write(n:1:2);
		end; {if}
		x:=x+xSpacing;
	end; {while}
	x:=gxOrigin;
	while (x>=frameL) do begin
		if (x<=frameR) then begin
			MakeLineRoi(x,frameT,x,frameB); Clear; KillRoi;
			n:=(x-frameL)*(gxMax-gxMin)/frameW+gxMin;
			MoveTo(x,y); write(n:1:2);
		end; {if}
		x:=x-xSpacing;
	end; {while}

	SetText('Right Justified, No Background');
	x:=frameL-txOffset;
	y:=gyOrigin;
	while (y<=frameB) do begin
		if (y>=frameT) then begin
			MakeLineRoi(frameL,y,frameR,y); Clear; KillRoi;
			n:=(frameB-y)*(gyMax-gyMin)/frameH+gyMin;
			MoveTo(x,y); write(n:1:2);
		end; {if}
		y:=y+ySpacing;
	end; {while}
	y:=gyOrigin;
	while (y>=frameT) do begin
		if (y<=frameB) then begin
			MakeLineRoi(frameL,y,frameR,y); Clear; KillRoi;
			n:=(frameB-y)*(gyMax-gyMin)/frameH+gyMin;
			MoveTo(x,y); write(n:1:2);
		end; {if}
		y:=y-ySpacing;
	end; {while}

	{¥draw frame}
	MoveTo(frameL,frameT);
	LineTo(frameR,frameT);
	LineTo(frameR,frameB);
	LineTo(frameL,frameB);
	LineTo(frameL,frameT);

	{¥draw axis}
	if (gxMin*gxMax)<0 then begin
		MoveTo(gxOrigin,frameB); LineTo(gxOrigin,frameT);
	end; {if};
	if (gyMin*gyMax)<0 then begin
		MoveTo(frameL,gyOrigin); LineTo(frameR,gyOrigin);
	end; {if}

	{¥plot the points}
	{}ShowMessage('Plotting points...');
	SetForegroundColor(215);
	for i:= 1 to rCount do begin
		x:=rUser1[i]; y:=rUser2[i];
		x:=round((rUser1[i]-gxMin)/(gxMax-gxMin)*frameW+frameL);
		y:=round(frameB-(rUser2[i]-gyMin)/(gyMax-gyMin)*frameH);
		if i=1 then MoveTo(x,y) else LineTo(x,y);
		MakeOvalRoi(x-cRad,y-cRad,2*cRad+1,2*cRad+1);fill;	
		DrawBoundary; KillRoi;
	end;{for i}

	{*Graph the curve}
	SetForegroundColor(235);
	ndx:=100;						{increase this to change the precision of the graph}
	DoLine:=false;
	done:=false;
	i:=0;
	dx:=(gxMax-gxMin)/ndx;
	repeat
		x:=gxMin+i*dx;
		if x>gxMax then begin
			x:=gxMax;
			done:=true;
		end; {if}
		y:=a + b*x + c*x*x;
		if not((y<=gyMax) and (y>=gyMin)) then DoLine:=false 
			else begin
				xp:=round((x-gxMin)/(gxMax-gxMin)*frameW+frameL);
				yp:=round(frameB-(y-gyMin)/(gyMax-gyMin)*frameH);
				if DoLine then LineTo(xp,yp) else MoveTo(xp,yp);
				DoLine:=true;
		end;
		i:=i+1;
	until done;
	RestoreState;

	{}ShowMessage('');

end; {MakeGraph}
{======================================================================}

macro '[Z] Analyze Motion...';

begin

	RequiresVersion(1.55);
{check for stack}
	if nPics=0 then Error('This operation requires a stack.');
	if nSlices=0 then Error('This window is not a stack.');
	theStack:=pidNumber;
	GetPicSize(Width,Height);

{Initialize some stuff}
	marginL:=50;	marginR:=10;					{default margin sizes for graphs}
	marginB:=50;	marginT:=10;
	nGridLines:=10;															{five horizontal/vertical grid lines on graphs}
	txOffset:=5;
	tyOffset:=7;
	wdt:= {434}380; hgt:={340}300;          {width and height of graph windows}
	cRad:=4;
	InvertY(true);

{Set the time info first}
	TimeUnit:=GetString('Units of time?','Seconds');
	repeat
	 TimeInt:=GetNumber('Time interval between slices?',0.10);
	 if not(TimeInt>0) then PutMessage('Enter a positive number!');
 until (TimeInt>0);

{Set the origin}
	xOrigin:=0; yOrigin:=Height;

	MoveWindow(84,40);
	SetOptions('X-Y Center');
	ResetCounter;
	ShowResults; MoveWindow(Width+88,40);
	SelectPic(theStack);

{Gather position data}
	GetScale(Scale,ScaleUnit);
 PutMessage('Click "OK" then click on ball in each frame.');
	SelectAll; KillRoi;
	for i:= 1 to nSlices do begin
		SelectSlice(i);
	 repeat
			ShowMessage('Click to select point ',i:1,'\X: ',x:1,'  Y: ',y:1);
  	SetCursor('Cross');
   GetMouse(x,y);
  until button;
	 repeat until not(button);
		x:=(x-xOrigin)/Scale;
		y:=(yOrigin-y)/Scale;	
		SetCounter(rCount+1);
  rX[rCount]:=x; rY[rCount]:=y;
		UpdateResults;
	end; {for}

{Find min/max}
	xMin:= rX[1]; yMin:=rY[1];
	xMax:= rX[1]; yMax:=rY[1];
	for i:= 2 to rCount do begin
		if rX[i]<xMin then xMin:=rX[i];
		if rX[i]>xMax then xMax:=rX[i];
		if rY[i]<yMin then yMin:=rY[i];
		if rY[i]>yMax then yMax:=rY[i];	
	end; {for}		

 tMin:=0; tMax:=(rCount-1)*TimeInt;

	ShowResults;

{graph X vs T}
	xLabel:=concat('Time (',TimeUnit,')');
	yLabel:=concat('X (',ScaleUnit,')');
	for i:= 1 to rCount do begin
		rUser2[i]:=rX[i];
		rUser1[i]:=(i-1)*TimeInt;
	end; {for i}
 CurveFit(0,rCount); c:=0;
	Equation :=concat('x = ',b:1:2,' t + ',a:1:2);
	MakeGraph(wdt,hgt,rCount,tMin,tMax,xMin,xMax,'X vs. T');
	MoveWindow(99,60);

{graph Y vs T}
	xLabel:=concat('Time (',TimeUnit,')');
	yLabel:=concat('Y (',ScaleUnit,')');
	for i:= 1 to rCount do begin
		rUser2[i]:=rY[i];
		rUser1[i]:=(i-1)*TimeInt;
	end; {for i}
	CurveFit(1,rCount);
	Equation := concat('y = ',c:1:2,' t^2 + ',b:1:2,' t + ',a:1:2);
	MakeGraph(wdt,hgt,rCount,tMin,tMax,yMin,yMax,'Y vs. T');
	MoveWindow(114,80);

{graph X vs Y}
	xLabel:=concat('X (',ScaleUnit,')');
	yLabel:=concat('Y (',ScaleUnit,')');
	for i:= 1 to rCount do begin
		rUser2[i]:=rY[i];
		rUser1[i]:=rX[i];
	end; {for i}
	CurveFit(1,rCount);
	Equation := concat('y = ',c:1:2,' x^2 + ',b:1:2,' x + ',a:1:2);
	MakeGraph(wdt,hgt,rCount,xMin,xMax,yMin,yMax,'Y vs. X');
	MoveWindow(129,100);

end;

{---------- +++++++++++++ ----------}


procedure CheckForPic;
{This procedure checks that there is at least one open image to draw on.}
begin
	if nPics=0 then Error('There are no open images.');
end;
{}
procedure sDrawBoundary(Filled:boolean);
{This procedure outlines the current selection with the current foreground (paintbrush) color. If Filled=true, it will also fill the inside of the boundary with the current background (eraser) color. If the current image is a stack and the global variable EverySlice=true drawing will be duplicated on every slice, otherwise it will only occur on the current slice.}

var i: integer;
begin
	if (nSlices>1) and EverySlice then begin
		for i:= 1 to nSlices do begin
			SelectSlice(i); 
			if Filled then Clear;
			DrawBoundary;
		end;
	end else begin
		if Filled then Clear;
		DrawBoundary;
	end;
	KillRoi;
end;

macro '(-'; begin end;
macro 'Set to draw on:'; begin end;
macro '		[S]	All Slices, or';
begin
	EverySlice:=true;
end;

macro '		[T]	This Slice Only';
begin
	EverySlice:=false;
end;
{}
macro 'Select Object to Draw:';
begin
	Error('Select an object from the list.'); 
end;
{}
macro '		[A]	Arrow';
{This macro will draw an arrow using the current line selection.}
var 	sScale,sLen,sAngle,Slope,Angle,sign:real;
							i,x1,y1,x2,y2,xHead,yHead,xTail,yTail,lw:integer;
begin
 RequiresVersion(1.55);
	CheckForPic;
	GetLine(xTail,yTail,xHead,yHead,lw);
	if (xTail=-1) then begin
		Error('Make a line selection first.');
	end; {if}
 sScale:=.2; sAngle:=.75;
	sLen:=10;			{sScale*sqrt(sqr(xTail-xHead)+sqr(yTail-yHead));}
	if xHead=xTail then begin
		Angle:=1.571;
		sign:=(yTail-yHead)/abs(yTail-yHead);
	end else begin
		sign:=(xTail-xHead)/abs(xHead-xTail);
		Slope:=(yHead-yTail)/(xHead-xTail);
		Angle:=arctan(Slope);
	end;
	x1:=xHead+sign*sLen*cos(Angle+sAngle);
	x2:=xHead+sign*sLen*cos(Angle-sAngle);
	y1:=yHead+sign*sLen*sin(Angle+sAngle);
	y2:=yHead+sign*sLen*sin(Angle-sAngle);
	if EverySlice then begin
		for i:= 1 to nSlices do begin
			SelectSlice(i);
			MoveTo(xTail,yTail); LineTo(xHead,yHead);
			LineTo(x1,y1); LineTo(x2,y2); LineTo(xHead,yHead);
		end; 
	end else begin
			MoveTo(xTail,yTail); LineTo(xHead,yHead);
			LineTo(x1,y1); LineTo(x2,y2); LineTo(xHead,yHead);
	end;
end;
{}
macro '		[C]	Circle From Center'
{Draws a circle from center out}
var X,Y,CX,CY,Dia,Rad :integer;
			scale:real;
			units:string;
begin
	RequiresVersion(1.55);
	GetScale(scale,units);
	SetCursor('cross');
	{Get Center}
	Repeat
		GetMouse(CX,CY);
		ShowMessage('Click at center  \(',CX:1,',',CY:1,')');
		{Wait(.1);}
	Until Button;
	Repeat until not(Button);
	Repeat
		GetMouse(X,Y);
		Rad:=sqrt(sqr(CX-X)+sqr(CY-Y));
		MakeOvalRoi(CX-Rad,CY-Rad,2*Rad,2*Rad);
		DrawBoundary;
		ShowMessage('Click when done... \Radius: ',Rad/scale:1:2,' ',units);
		{Wait(.08)};
		Undo;
	Until Button;
	repeat until not(Button);
	sDrawBoundary(false);
	KillRoi;
end;

macro '		[O]	Outline';
{This macro draws the current region of interest (marching ants selection)
in the current foreground (paintbrush) color using the current line width.}
begin
	RequiresVersion(1.55);
	CheckForPic;
	sDrawBoundary(false);
end;
{}
macro '		[F]	Outline and Fill';
{
This macro draws the current region of interest (marching ants selection)
in the current foreground (paintbrush) color using the current line width.
It also fills in the region with the current background (eraser) color.}
begin
 RequiresVersion(1.55);
	CheckForPic;
	sDrawBoundary(true);
end;
{}
macro '		[H]	Horizontal Line';
{This macro draws a horizontal line the point selected with the mouse.
It uses the current foreground (paintbrush) color and line width.}
var x,y,h,w:integer;
begin
 RequiresVersion(1.55);
	CheckForPic;
	SetCursor('Cross');
	Repeat
		GetMouse(x,y);
		ShowMessage('Click on point on line...','\\X: ',x:1,'\Y: ',y:1);
	until Button;
	Repeat until not(Button);
	GetPicSize(w,h);
	MakeLineRoi(0,y,w,y);
	sDrawBoundary(false);
end;
{}
macro '		[V]	Vertical Line';
{This macro draws a vertical line the point selected with the mouse.
It uses the current foreground (paintbrush) color and line width.}
var x,y,h,w:integer;
begin
 RequiresVersion(1.55);
	CheckForPic;
	SetCursor('Cross');
	Repeat
		GetMouse(x,y);
		ShowMessage('Click on point on line...','\\X: ',x:1,'\Y: ',y:1);
	until Button;
	Repeat until not(Button);
	GetPicSize(w,h);
	MakeLineRoi(x,0,x,h);
	sDrawBoundary(false);
end;
{}
macro '		[X]	Axes';
{This macro draws both a horizontal and a vertical line the point selected
with the mouse. It uses the current foreground (paintbrush) color and line width.}
var x,y,w,h:integer;
begin
 RequiresVersion(1.55);
	CheckForPic;
	SetCursor('Cross');
	ShowMessage('Click on origin...');
	Repeat
		GetMouse(x,y);
	until Button;
	Repeat until not(Button);
	GetPicSize(w,h);
	MakeLineRoi(x,0,x,h); sDrawBoundary(false);
	MakeLineRoi(0,y,w,y);	sDrawBoundary(false);
end;
{}
procedure ClearOutside;
begin
  RequiresVersion(1.55);
		Copy; SelectAll; Clear;
  RestoreRoi; Paste; 
end;
{}
procedure DrawGrid
begin
  RequiresVersion(1.55);
 	x:=0; y:=0;
  repeat
    x:=x+xinc;  y:=y+yinc;
    moveto(0,y); lineto(width,y);
    moveto(x,0); lineto(x,height);
  until (x>width) and (y>height);
end;
{}
macro '		[G]	GridÉ';
var
  x,y,xinc,yinc,width,height,i:integer;
begin
  RequiresVersion(1.55);
  GetPicSize(width,height);
  xinc:=GetNumber('Horizontal Spacing:',16);
  yinc:=GetNumber('Vertical Spacing:',xinc);
	if (nSlices>1) and EverySlice then begin
		for i:=1 to nSlices do begin
			SelectSlice(i);	DrawGrid;
		end;
	end else DrawGrid;
end;
{}
macro '		[L]	Clear Outside';
var i:integer;
 {Erase region outside current selection to background color.}
begin
	RequiresVersion(1.55);
	if (nSlices>1) and EverySlice then begin
		for i:=1 to nSlices do begin
			SelectSlice(i); RestoreRoi;
			ClearOutside; 
		end;
	end else ClearOutside;
	KillRoi; 
end;