////////////////////////////////////////////////////
//  ImageJ macros for flagellar beat analysis
//
//  Scripted by Bill Wickstead
//  v3.0 March 2007


////////////////////////////////////////////////////
//   GLOBAL VARIABLES   //
//  constants
var pi=3.1415926536;

//  button-click variables
var leftButton=16;
var rightButton=4;
var shift=1, ctrl=2, alt=8;
var altClick=24, shiftClick=17;
var pencilWidth=1, brushWidth=10;

////////////////////////////////////////////////////
//   TIME STAMPER VARIABLES   //
var TSpositionX="right", TSpositionY="bottom";
var TSmarginX=10, TSmarginY=10;
var TSfontSize=12, TSstartTime=0, TSstep=5.208333;
var TSsuffix=' ms', TSdecimals=1;
var TSfont="SansSerif";

////////////////////////////////////////////////////
//   ARROW THRU STACK VARIABLES   //
var arrowSize=12,  arrowAngle=0;
var arrowSelectSize=false, arrowSelectAngle=true;
var arrowAllSlice=true, arrowDupStack=false;

////////////////////////////////////////////////////
//   SWIMMING ON THE SPOT VARIABLES   //
var sliceStep=3;
var trackingPoints=2;
var centreAlign=3;		// centreAlign options [3=do both]
var doTicks=false;		// make cell displacement ruler
var tickY=100;			// position of ruler (below centre)
var tickHeight=10;
var tickSpace=266.7;		// pixels
var keepGraphs=false		// keep fitting graphs for further use

////////////////////////////////////////////////////
//   CURVATURE VARIABLES   //
var curvCal=0;		// um/pix value; 0=none
var curvSmooth=2;	// smoothing; 0=none
var curvPlot=true;	// create plots
var curvPlotFix;		// fixed y axis
var curvYmin, curvYmax;
var curvClearResults=true;

////////////////////////////////////////////////////
//   BEAT ANALYSIS VARIABLES   //
var beatStep=20;
var beatSpeed=20;
var markerSize=5;
var sampleSpline=4;


////////////////////////////////////////////////////
//   MATHS FUNCTIONS   //
function int(x) {  y=round( x - (x % 1) );  return y;  }
function roundUp(x) { 
	y=int(x); 
	if (x % y >0) { y++; } else if (x % y <0) { y--; }
	return y; }
function sign(x) { if (x<0) { return -1; } else if (x>0) { return 1; } else { return 0; } }
function magnitude(x) {  // returns magnitude of x
	x=abs(x);
	if (x==0) { return 0; }
	else if (x<1) { y=int(-1+log(x)/log(10));  return y; }
	else { y=int(log(x)/log(10));  return y; }  }
function magScale(x) {  // returns rounded up 1st sig.fig. of x
	y=roundUp( x/ pow(10, magnitude(x)) );  return y; }  }

macro '--- to use shortcuts hit alt+[key] ---' {}

////////////////////////////////////////////////////
//  DrawingTools
//
//  This is a set of drawing tools similar to the pencil, paintbrush, eraser and paint bucket tools in NIH Image. 
//  Pencil and paintbrush draw in the current foreground color, eraser draws in the current background color.
//  Paintbucket tool fills the selected area using the foreground color.
//  Alt key + pencil, paintbrush or paintbucket tool uses background color.
//  Double-click on the paintbucket or eye-dropper tool to set colors.
//  Double-click on the pencil, paintbrush or eraser tool to set the drawing width for that tool.
//
//  Icons contributed by Tony Collins.

  macro "Unused Tool-1 - " {}  // leave slot between text tool and magnifying glass unused
  macro "Unused Tool-2 - " {}  // leave slot between dropper and pencil unused

  macro "Pencil Tool - C037L494fL4990L90b0Lc1c3L82a4Lb58bL7c4fDb4L5a5dL6b6cD7b" {
        getCursorLoc(x, y, z, flags);
        if (flags&alt!=0) setColorToBackgound();
        draw(pencilWidth); }

   macro "Paintbrush Tool - C037F6036F3699CfffD71F4771D5eD7eD9e" {
        getCursorLoc(x, y, z, flags);
        if (flags&alt!=0) setColorToBackgound();
        draw(brushWidth); }

   macro "Eraser Tool - C037R0aa4 P0a61f1aa0Pbef5f1" { 
        setColorToBackgound();
        draw(eraserWidth); }

    macro "Flood Fill Tool -C037B21P085373b75d0L4d1aL3135L4050L6166D57D77D68La5adLb6bcD09D94" {
        requires("1.34j");
        setupUndo();
        getCursorLoc(x, y, z, flags);
        if (flags&alt!=0) setColorToBackgound();
        floodFill(x, y); }

   
function draw(width) {
        requires("1.32g");
        setupUndo();
        getCursorLoc(x, y, z, flags);
        setLineWidth(width);
        moveTo(x,y);
        x2=-1; y2=-1;
        while (true) {
            getCursorLoc(x, y, z, flags);
            if (flags&leftButton==0) exit();
            if (x!=x2 || y!=y2) lineTo(x,y);
            x2=x; y2 =y;
            wait(10); } }

   function setColorToBackgound() {
       savep = getPixel(0, 0);
       makeRectangle(0, 0, 1, 1);
       run("Clear");
       background = getPixel(0, 0);
       run("Select None");
       setPixel(0, 0, savep);
       setColor(background); }


macro 'Extra tool options...' {
	choice=newArray("pencil", "paintbrush", "eraser", "floodfill");
	Dialog.create("tool options");
	Dialog.addChoice("Tool: ", choice);
	Dialog.show();
	choice=Dialog.getChoice();
	if (choice=="pencil") { pencilWidth = getNumber("Pencil Width (pixels):", pencilWidth); }
	else if (choice=="paintbrush") { brushWidth = getNumber("Brush Width (pixels):", brushWidth); }
	else if (choice=="eraser") { eraserWidth = getNumber("Eraser Width (pixels):", eraserWidth); }
	else if (choice=="floodfill") { requires("1.34j"); restorePreviousTool(); run("Color Picker..."); }
}	
  // Runs when the user double-clicks on the pencil tool icon
 // macro 'Pencil Tool Options...' { pencilWidth = getNumber("Pencil Width (pixels):", pencilWidth); }
  // Runs when the user double-clicks on the paint brush tool icon
 // macro 'Paintbrush Tool Options...' { brushWidth = getNumber("Brush Width (pixels):", brushWidth); }
  // Runs when the user double-clicks on the eraser tool icon
 // macro 'Eraser Tool Options...' { eraserWidth = getNumber("Eraser Width (pixels):", eraserWidth); }
  // Runs when the user double-clicks on the flood fill tool icon
 // macro 'Flood Fill Tool Options...' { requires("1.34j"); restorePreviousTool(); run("Color Picker..."); }


////////////////////////////////////////////////////
//  General macros tools
//
//  Written by Bill Wickstead Jan 2006

macro '         --general--' {}

macro 'Dispose All [d]' {
	if (nImages<1) { exit("no images open"); }
	if (nImages>1) ok= getBoolean("OK to dispose of " + nImages + ' images without saving?');
	else ok= getBoolean('OK to dispose of image without saving?');
	if (!ok) { exit (); }
	i=nImages;
	while (i>0) { selectImage(i); close; i--; }
}

macro 'Enhance constrast [h]' { run ("Enhance Contrast", "Saturated Pixels=0 Normalize"); }
macro 'Apply LUT [a]' { run("Apply LUT"); }

macro 'Capture points from spline [S]' {
	s=sampleSpline;
	Dialog.create("Spline options");
	Dialog.addNumber("Sample every n-th points:", s);
	Dialog.addCheckbox("Clear results:", true);
	Dialog.show;
	s=Dialog.getNumber();
	clearRes=Dialog.getCheckbox();
	if (s<1 || s>50) { showMessage("Out of range"); exit; }
	if (clearRes) { run("Clear Results"); }  NR=nResults;
	run("Fit Spline");
	getSelectionCoordinates(x,y);	
	for (n=0+NR; n<x.length+NR; n+=s) {
		setResult("x", NR+(n-NR)/s, x[n-NR]);  setResult("y", NR+(n-NR)/s, getHeight-y[n-NR]);
	}
	updateResults();
}

macro 'Crop to... [c]' {
  width=round(getWidth()/2);
  height=round(getHeight()/2);
  continue=1;
  while (continue==1) {
    Dialog.create("Crop to");
    Dialog.addNumber("Width:", width);
    Dialog.addNumber("Height:", height);
    Dialog.show();
    width=Dialog.getNumber();
    height=Dialog.getNumber();
    if (width>getWidth() || height>getHeight()) {
      message="Out of bounds\nPicture size: "+getWidth()+"x"+getHeight();
      showMessage(message);
    } else { continue=0; }
  }

  left=(getWidth() - width)/2;
  top=(getHeight() - height)/2;
  makeRectangle(left, top, width, height);
  run("Crop");
}
macro 'Draw centre line' {
	choice=newArray("black", "white");
	Dialog.create("Centre line");
	Dialog.addNumber("Line width:", 1);
	Dialog.addChoice ("line colour: ", choice);
	Dialog.show();
	lw=Dialog.getNumber();
	col=Dialog.getChoice();
	if (col=="black") { setColor(0); }
	else if (col=="white") { setColor(255); }
	else { showMessage("out of bounds"); exit; }
	w=getWidth();  h=getHeight();  NS=nSlices;
	setLineWidth(lw);
	for (n=1; n<=NS; n++) {
		setSlice(n);
		drawLine(0,h/2,w-1,h/2);
		drawLine(w/2,0,w/2,h-1);		
	}
	setSlice(1);
}

////////////////////////////////////////////////////
//  Grayscale manipulation macros
//
//  Written by Bill Wickstead April 2006

macro '         --levels--' {}
macro 'Lighten / Darken...' {
	choice=newArray("lighten", "darken");
	Dialog.create("Lighten / Darken");
	Dialog.addChoice("Lighten / Darken", choice);
	Dialog.addNumber("All pixels by:", 0);
	Dialog.show();
	choice=Dialog.getChoice();
	AddVal=Dialog.getNumber();
	if (choice=="lighten") {
		if (nSlices>1) { run("Add...", "stack value="+AddVal); }
		else { run("Add...", "value="+AddVal); }
	} else {
		if (nSlices>1) { run("Subtract...", "stack value="+AddVal); }
		else { run("Subtract...", "value="+AddVal); }
	}
}
macro 'Remove non-black' {
	dupl=true;
	tempArray=split(getTitle(), '.');
	name=tempArray[0]+'_nonblack';
	if (dupl) { run("Duplicate...", "title="+name+" duplicate"); } else { rename(name); }
	if (nSlices>1) { run("Multiply...", "stack value=255"); }
	else { run("Subtract...", "value=255"); }
}


////////////////////////////////////////////////////
//  Stack manipulation macros
//
//  Written by Bill Wickstead Jan 2006

function NoStackError() {
  if (nImages<1) { exit("image required"); }
  if (nSlices<1) { exit("stack required"); }
}

function NoSelectionError() {
  if (selectionType()==-1) { exit("selection required"); }
}

function PictureType() {
  info=getImageInfo();
  index1=indexOf(info, 'Bits per pixel: ');
  index1 += 16;
  index2=indexOf(info, '(', index1);
  type=substring(info, index1, index2-1)+'-bit';
  print(substring(info, index2+1, index2+4));
  if (substring(info, index2+1, index2+4)=='RGB') { type='RGB'; }
  return type;
}

function sureSelect(ID, name) {
// Makes sure selectImage() command goes to the correct image
// selectImage() error appears to be a bug in versions 1.3* of ImageJ
// function should become redundant
	selectImage(ID);
	checkID=getImageID();
	checkname=getTitle();
	if (ID != checkID || checkname != name) {
		exit ("Image select failure\nRequested: ID: "+ID+"   named: "+name+"\nGot:     ID: "+checkID+"   named: "+checkname);
	}
}

macro '         --stacks--' {}
macro 'Goto first slice [f]' { setSlice(1); }
macro 'Goto last slice [l]' { setSlice(nSlices); }

macro 'Goto slice... [g]' {
	NoStackError();
	n=round(nSlices/2);
	continue=1;
	while (continue==1) {
		Dialog.create("Select slice");
		Dialog.addMessage("Number of slices: "+nSlices);
		Dialog.addNumber("Goto slice:", n);
		Dialog.show();
		n=Dialog.getNumber();
		if (n<1 || n>nSlices()) {
			message="Out of bounds\nNumber of Slices: "+nSlices;
			showMessage(message);
		} else { continue=0; }
	}
	setSlice(n);
}

macro 'Slice to new window [n]' {
  run ("Select All");
  run ("Copy");
  name="slice"+getSliceNumber();
  type=PictureType();
  newImage(name, type, getWidth(), getHeight(), 1);
  run ("Paste");
}

macro 'Sample slices...' {
	NoStackError();
	NS=nSlices;  continue=false;
	while (!continue) {
		Dialog.create("Sample slices");
		Dialog.addNumber("Sample every n-th slice:", 5);
		Dialog.addNumber("Start slice", 1);
		Dialog.show();
		sampleFreq=Dialog.getNumber();
		startSlice=Dialog.getNumber();
		if (sampleFreq<1 || sampleFreq>NS) {
			message="Out of bounds\nNumber of Slices: "+NS;
			showMessage(message);
		} else if (startSlice<1 || startSlice>NS) {
			message="Out of bounds\nNumber of Slices: "+NS;
			showMessage(message);
		} else { continue=true; }
	}
	
	tempArray=split(getTitle(), '.');
	name=tempArray[0]+'_sampled';
	run("Duplicate...", "title="+name+" duplicate");
	setSlice(1);
	for (i=1; i<startSlice; i++) { run("Delete Slice"); }
	NS=nSlices;
	i=0;  sliceNum=1;
	while (i<NS) {
		setSlice(sliceNum);
		if ( i % sampleFreq ==0) { sliceNum++;  i++; }
		else { run("Delete Slice");  i++; }
	}
	setSlice(1);
}

macro 'Stack trimmer...' {
	dupl=true;  firstSlice=1;  lastSlice=nSlices;
	Dialog.create("Trim stack");
	Dialog.addNumber("First slice: ", firstSlice);
	Dialog.addNumber("Last slice: ", lastSlice);
	Dialog.addCheckbox("create duplicate stack ", dupl);
	Dialog.show();
	firstSlice=Dialog.getNumber();
	lastSlice=Dialog.getNumber();
	dupl=Dialog.getCheckbox();
	if (firstSlice<1 || lastSlice>nSlices) { exit("Error: selected slices out of bounds"); }
	tempArray=split(getTitle(), '.');
	name=tempArray[0]+'_trimmed';
	if (dupl) { run("Duplicate...", "title="+name+" duplicate"); } else { rename(name); }
	i=nSlices;  setSlice(i);
	while (i>lastSlice) { run("Delete Slice");  i--; }
	i=1;  setSlice(i);
	while (i<firstSlice) { run("Delete Slice");  i++; }
}

macro 'Split stack...' {
	splitAt=round(nSlices/2);
	splitAt=getNumber("Split after slice: ", splitAt);
	if (splitAt >nSlices) { exit("Error: split > nSlices"); }
	run("Select All");
	tempArray=split(getTitle(), '.');
	name1=tempArray[0]+'_part1';
	name2=tempArray[0]+'_part2';
	rename(name1);
	part1=getImageID();
	run("Duplicate...", "title="+name2+" duplicate");
	part2=getImageID();
	for (i=1; i<=splitAt; i++) { setSlice(1);  run("Delete Slice"); }
	selectImage(part1);
	for (i=nSlices; i>splitAt; i--) { setSlice(i);  run("Delete Slice"); }
	selectImage(part1);  setSlice(1);
	selectImage(part2);  setSlice(1);
}		// End of macro

macro 'Time stamper... [t]' {
	NoStackError();
	Dialog.create("Time stamper");
	Dialog.addNumber("Start time:",TSstartTime);
	Dialog.addNumber("Time between frames:",TSstep);
	tempArray=newArray("SansSerif", "Monospaced", "Serif");
	choice=newArray(tempArray.length);  choice[0]=TSfont;
	n=0;  for (i=1; i<tempArray.length; i++) { 
		if (tempArray[n]==TSfont) { n++; }
		choice[i]=tempArray[n]; n++; }
	Dialog.addChoice("Font:",choice);
	Dialog.addNumber("Font size:", TSfontSize);
	Dialog.addNumber("Decimal places:", TSdecimals);
	Dialog.addString("Suffix:", TSsuffix);
	tempArray=newArray("left", "center", "right");
	choice=newArray(tempArray.length);  choice[0]=TSpositionX;
	n=0;  for (i=1; i<tempArray.length; i++) { 
		if (tempArray[n]==TSpositionX) { n++; }
		choice[i]=tempArray[n]; n++; }
	Dialog.addChoice("Position X:", choice);
	tempArray=newArray("top", "center", "bottom");
	choice=newArray(tempArray.length);  choice[0]=TSpositionY;
	n=0;  for (i=1; i<tempArray.length; i++) { 
		if (tempArray[n]==TSpositionY) { n++; }
		choice[i]=tempArray[n]; n++; }
	Dialog.addChoice("Position Y:", choice);
	Dialog.addNumber("X margin:", TSmarginX);
	Dialog.addNumber("Y margin:", TSmarginY);
	Dialog.show();

	TSstartTime=Dialog.getNumber();
	TSstep=Dialog.getNumber();
	TSfont=Dialog.getChoice();
	TSfontSize=Dialog.getNumber();
	TSdecimals=Dialog.getNumber();
	TSsuffix=Dialog.getString();
	TSpositionX=Dialog.getChoice();
	TSpositionY=Dialog.getChoice();
	TSmarginX=Dialog.getNumber();
	TSmarginY=Dialog.getNumber();

	setFont(TSfont, TSfontSize, 'Bold');
	setJustification(TSpositionX);
	if (TSpositionX=="left") { x=TSmarginX; }
		else if (TSpositionX=="center") { x=getWidth()/2; }
		else if (TSpositionX=="right") { x=getWidth()-TSmarginX; }
		else { showMessage("Error #201");  exit; }
	if (TSpositionY=="top") { y=TSmarginY+TSfontSize; }
		else if (TSpositionY=="center") { y=getHeight()/2+TSfontSize/2; }
		else if (TSpositionY=="bottom") { y=getHeight()-TSmarginY; }
		else { showMessage("Error #202");  exit; }

	time=TSstartTime;
	for (i=1; i<=nSlices; i++) {
		setSlice(i); 
		if (TSdecimals>0) {
			str=toString(int(time))+'.';  rem=time % int(time);
			for (n=1; n<TSdecimals; n++) { rem=rem*10;  str=str+toString( int(rem) );  rem=rem % int(rem); }
			rem=round(rem*10);  if (rem==10) { str=str+toString(1); } else { str=str+toString(rem); }
		} else { str=toString( round(time) ); }
		str=str+TSsuffix;
		drawString(str, x, y);
		time+=TSstep;
	}
	setSlice(1);
}		// End of macro


////////////////////////////////////////////////////
//  Stack transformation macros
//
//  Written by Bill Wickstead March 2006

macro '         --stack transformations--' {}

macro 'Fix point(s) thru stack... [F]' {
	NoStackError();
	choice=newArray("1", "2");
	Dialog.create("Fix point(s) thru stack");
	Dialog.addChoice("Number of fixed points per slice", choice);
	Dialog.addCheckbox("Centre:", true);
	Dialog.addCheckbox("Align to x-axis:", false);
	Dialog.addCheckbox("Interpolate:", true);
	Dialog.addCheckbox("Create new window:", true);
	Dialog.show();
	choice=Dialog.getChoice();
	centreNew=Dialog.getCheckbox();
	alignNew=Dialog.getCheckbox();
	interpol=Dialog.getCheckbox();
	dupl =Dialog.getCheckbox();

	NS=nSlices;  	tempArray=split(getTitle(), '.');  name=tempArray[0]+'_fixed';
	if (dupl) { run("Duplicate...", "title="+name+" duplicate"); } else { rename(name); }
	w=getWidth();  h=getHeight();

	x1=newArray(NS);  y1=newArray(NS);  x2=newArray(NS);  y2=newArray(NS);
	theta=newArray(NS);  setTool(7);
	if (choice=="1") { showMessageWithCancel("Please track point to fix thru stack"); }
		else { showMessageWithCancel("Please track 1st point to fix thru stack"); }
// Capture points
	setTool(7);
// proximal
	flagsOld=-1;
	i=0; while (i<NS) {
		setSlice(i+1);
		getCursorLoc(x, y, z, flags);
		if (flags != flagsOld) {
			if (flags == altClick) { 
				if (i>0) { 
					i--;  pointX=newArray(i);  pointY=newArray(i);
					for (n=0; n<i; n++) { pointX[n]=x1[n];  pointY[n]=y1[n]; }
				}
			} else if (flags == leftButton) {
				x1[i]=x;  y1[i]=y;  i++;
				pointX=newArray(i);  pointY=newArray(i);
				for (n=0; n<i; n++) { pointX[n]=x1[n];  pointY[n]=y1[n]; }
			}
		}  flagsOld=flags;
	}
// distal
	if (choice=="2") {
		showMessageWithCancel("Please track 2nd point to fix thru stack");  
		flagsOld=-1;
		i=0; while (i<NS) {
			setSlice(i+1);
			getCursorLoc(x, y, z, flags);
			if (flags != flagsOld) {
			if (flags == altClick) { 
				if (i>0) { 
					i--;  pointX=newArray(i);  pointY=newArray(i);
					for (n=0; n<i; n++) { pointX[n]=x2[n];  pointY[n]=y2[n]; }
				}
			} else if (flags == leftButton) {
				x2[i]=x;  y2[i]=y;  i++;
				pointX=newArray(i);  pointY=newArray(i);
				for (n=0; n<i; n++) { pointX[n]=x2[n];  pointY[n]=y2[n]; }
			}
			}  flagsOld=flags;
		}
	}

	if (choice=="2") {
		for (n=0; n<NS; n++) { theta[n]=atan2(y1[n]-y2[n], x2[n]-x1[n]); }
// swell image such that no piece leaves picture on rotate
		bigW=0;  bigH=0;
		for (n=0; n<NS; n++) {
			if (alignNew) {  testW=abs(w*cos(theta[n]))+abs(h*sin(theta[n])); 
				testH=abs(h*cos(theta[n]))+abs(w*sin(theta[n])); }
			else { testW=abs(w*cos(theta[n]-theta[0]))+abs(h*sin(theta[n]-theta[0]));
				testH=abs(h*cos(theta[n]-theta[0]))+abs(w*sin(theta[n]-theta[0])); }
			if (testW>bigW) { bigW=testW; }
			if (testH>bigH) { bigH=testH; }
		}
		for (n=0; n<NS; n++) {
			x1[n]+=(bigW-w)/2;  y1[n]+=(bigH-h)/2;
			x2[n]+=(bigW-w)/2;  y2[n]+=(bigH-h)/2; }
		run("Canvas Size...", "width="+bigW+" height="+bigH+" position=Center");
// Rotate
		for (n=0; n<NS; n++) {
			if (alignNew) { angle=180*theta[n]/pi; } else { angle=180*(theta[n]-theta[0])/pi; }
			setSlice(n+1);  run("Select All");
			if (interpol) { run("Arbitrarily...", "slice angle="+angle+" interpolate"); } else { run("Arbitrarily...", "slice angle="+angle); }
		}
// swell again such that no piece leaves on move
		minX=0;  maxX=0;  minY=0;  maxY=0;
		newX=newArray(NS);  newY=newArray(NS);
		for (n=0; n<NS; n++) {
			if (alignNew) {
				newX[n]=(x1[n]-bigW/2)*cos(theta[n]) + (bigH/2-y1[n])*sin(theta[n]) + bigW/2;
				newY[n]=bigH/2-(bigH/2-y1[n])*cos(theta[n]) + (x1[n]-bigW/2)*sin(theta[n]);
			} else {
				newX[n]=(x1[n]-bigW/2)*cos(theta[n]-theta[0]) + (bigH/2-y1[n])*sin(theta[n]-theta[0]) +bigW/2;
				newY[n]=bigH/2-(bigH/2-y1[n])*cos(theta[n]-theta[0]) + (x1[n]-bigW/2)*sin(theta[n]-theta[0]);
			}
			if (centreNew) { 
				delx=abs(-newX[n]+bigW/2);  dely=abs(-newY[n]+bigH/2);
				if (delx>maxX) { maxX=delx; }
				if (dely>maxY) { maxY=dely; }
			} else { 
				delx=x1[0]-newX[n];  dely=y1[0]-newY[n];
				if (delx<minX) { minX=delx; } else if (delx>maxX) { maxX=delx; }
				if (dely<minY) { minY=dely; } else if (dely>maxY) { maxY=dely; }
			}
		}
		if (centreNew) { minX=-maxX;  minY=-maxY; }
		run("Canvas Size...", "width="+bigW-minX+" height="+bigH-minY+" position=Bottom-Right");
		run("Canvas Size...", "width="+bigW-minX+maxX+" height="+bigH-minY+maxY+" position=Top-Left");
// Move
		newW=getWidth();  newH=getHeight();
		for (n=0; n<NS; n++) {
			newX[n]-=minX;  newY[n]-=minY;
			setSlice(n+1);
			if (centreNew) { delx=-newX[n]+newW/2;  dely=-newY[n]+newH/2; }
				else { delx=x1[0]-minX-newX[n];  dely=y1[0]-minX-newY[n]; }
			run("Select All");  run("Cut");
			makeRectangle(delx, dely,newW, newH);  run("Paste");
		}
	} else {
// swell such that no piece leaves on move
		minX=0;  maxX=0;  minY=0;  maxY=0;
		if (centreNew) {
			for (n=0; n<NS; n++) {
				delx=abs(-x1[n]+w/2);  dely=abs(-y1[n]+h/2); 
				if (delx>maxX) { maxX=delx; }
				if (dely>maxY) { maxY=dely; }
			}
			minX=-maxX;  minY=-maxY;
		} else {
			for (n=0; n<NS; n++) {
				delx=x1[0]-x1[n];  dely=y1[0]-y1[n];
				if (delx<minX) { minX=delx; } else if (delx>maxX) { maxX=delx; }
				if (dely<minY) { minY=dely; } else if (dely>maxY) { maxY=dely; }
			}
		}
		run("Canvas Size...", "width="+w-minX+" height="+h-minY+" position=Bottom-Right");
		run("Canvas Size...", "width="+w-minX+maxX+" height="+h-minY+maxY+" position=Top-Left");
// Move
		newW=getWidth();  newH=getHeight();
		for (n=0; n<NS; n++) {
			x1[n]-=minX;  y1[n]-=minY;
			setSlice(n+1);
			if (centreNew) { delx=-x1[n]+newW/2;  dely=-y1[n]+newH/2; }
				else { delx=x1[0]-x1[n];  dely=y1[0]-y1[n]; }
			run("Select All");  run("Cut");
			makeRectangle(delx, dely,newW,newH);  run("Paste");
		}  
	}
	setSlice(1);  run("Select None");
}  			// end of macro

macro 'Align 2 stacks...' {
	names=newArray(nImages);
	for (i=0; i<nImages; i++) { selectImage(i+1);  names[i]=getTitle(); }
	Dialog.create("Choose stacks to align");
	Dialog.addChoice("First stack: ", names);
	Dialog.addChoice("Second stack: ", names);
	Dialog.show();
	name1=Dialog.getChoice();
	name2=Dialog.getChoice();
	showMessage("Select 2 points on each image for alignment");
	x1=newArray(2);  y1=newArray(2);  x2=newArray(2);  y2=newArray(2);  setTool(7);

	selectImage(name1); 
	NS1=nSlices;   w1=getWidth();   h1=getHeight();  type1=PictureType();  ID1=getImageID();
	setSlice(NS1);  flagsOld=-1;
	i=0; while (i<2) {
		getCursorLoc(x, y, z, flags);
		if (flags != flagsOld) {
			if (flags == leftButton) { x1[i]=x;  y1[i]=y;  i++; }
		}  flagsOld=flags;
	}
	selectImage(name2);
	NS2=nSlices;   w2=getWidth();   h2=getHeight();  type2=PictureType();  ID2=getImageID();
	setSlice(1);  flagsOld=-1;
	i=0; while (i<2) {
		getCursorLoc(x, y, z, flags);
		if (flags != flagsOld) {
			if (flags == leftButton) { x2[i]=x;  y2[i]=y;  i++; }
		}  flagsOld=flags;
	}

// swell smaller image
	if (w1>w2) { w3=w1; } else { w3=w2; }
	if (h1>h2) { h3=h1; } else { h3=h2; }
	selectImage(ID1);  run("Canvas Size...", "width="+w3+" height="+h3+" position=Center");
	selectImage(ID2);  run("Canvas Size...", "width="+w3+" height="+h3+" position=Center");
	for (n=0; n<2; n++) { 
		x1[n]+=(w3-w1)/2;  y1[n]+=(h3-h1)/2;
		x2[n]+=(w3-w2)/2;  y2[n]+=(h3-h2)/2;
	}

	theta1=atan2(y1[0]-y1[1], x1[1]-x1[0]);
	theta2=atan2(y2[0]-y2[1], x2[1]-x2[0]);
	angle=180*(theta2-theta1)/pi;

	bigW=abs(w3*cos(theta2-theta1))+abs(h3*sin(theta2-theta1));
	bigH=abs(h3*cos(theta2-theta1))+abs(w3*sin(theta2-theta1));
	for (n=0; n<2; n++) { x2[n]+=(bigW-w3)/2;  y2[n]+=(bigH-h3)/2; }
	run("Canvas Size...", "width="+bigW+" height="+bigH+" position=Center");

// Rotate
	for (i=1; i<=NS2; i++) {
		setSlice(i);  run("Select All");
		run("Arbitrarily...", "slice angle="+angle+" interpolate");
	}
// Move
	newX= (x2[0]-bigW/2)*cos(theta2-theta1) + (bigH/2-y2[0])*sin(theta2-theta1) +bigW/2;
	newY= bigH/2 - (bigH/2-y2[0])*cos(theta2-theta1) + (x2[0]-bigW/2)*sin(theta2-theta1);
	delx=x1[0]-newX;  dely=y1[0]-newY;
	for (i=1; i<=NS2; i++) {
		setSlice(i);  run("Select All");  run("Cut");
		makeRectangle(delx, dely, bigW, bigH);  run("Paste");
	}

	makeRectangle(0,0,w3,h3);  run("Crop");
	setSlice(1);	
}


////////////////////////////////////////////////////
//  Swimming on the spot macros
//
//  Written by Bill Wickstead March 2006

macro '         --swimming on the spot--' {}

macro 'Transform to stationary [T]' {
// Gathers info on movement of 1 or 2 points thru a stack, fits to polynomial
// and transforms to stationary swimming
	NoStackError();
	NS=nSlices;   width=getWidth();   height=getHeight();
	type1=PictureType();   name1=getTitle();  ID1=getImageID();

// Capture points
	setTool(7);
	len=round(0.5+NS/sliceStep);
	if (NS % sliceStep == 0) { len--; }
// proximal
	if (centreAlign==0) { showMessageWithCancel("Cell displacement calc:\nTrack single point thru stack"); }
	else {
		if (trackingPoints==1) { showMessageWithCancel("Track single point thru stack"); }
		else if (trackingPoints==2) { showMessageWithCancel("Track PROXIMAL point thru stack"); }
		else { exit("Error in trackingPoints value"); }
	}
	proxX=newArray(len);  proxY=newArray(len);  flagsOld=-1;
	i=0; while (i<len) {
		setSlice(i*sliceStep+1);
		getCursorLoc(x1, y1, z1, flags);
		if (flags != flagsOld) {
			if (flags == altClick) { 
				if (i>0) { 
					i--;  pointX=newArray(i);  pointY=newArray(i);
					for (n=0; n<i; n++) { pointX[n]=proxX[n];  pointY[n]=proxY[n]; }
					makeSelection("point", pointX, pointY);
				}
			} else if (flags == leftButton) {
				proxX[i]=x1;  proxY[i]=y1;  i++;
				pointX=newArray(i);  pointY=newArray(i);
				for (n=0; n<i; n++) { pointX[n]=proxX[n];  pointY[n]=proxY[n]; }
				makeSelection("point", pointX, pointY);
			}
		}  flagsOld=flags;
	}
// distal
	if (trackingPoints==2) {
		showMessageWithCancel("Track DISTAL point thru stack"); 
		distX=newArray(len);  distY=newArray(len);  flagsOld=-1;
		i=0; while (i<len) {
			setSlice(i*sliceStep+1);
			getCursorLoc(x1, y1, z1, flags);
			if (flags != flagsOld) {
			if (flags == altClick) { 
				if (i>0) { 
					i--;  pointX=newArray(i);  pointY=newArray(i);
					for (n=0; n<i; n++) { pointX[n]= distX[n];  pointY[n]=distY[n]; }
					makeSelection("point", pointX, pointY);
				}
			} else if (flags == leftButton) {
				distX[i]=x1;  distY[i]=y1;  i++;
				pointX=newArray(i);  pointY=newArray(i);
				for (n=0; n<i; n++) { pointX[n]=distX[n];  pointY[n]=distY[n]; }
				makeSelection("point", pointX, pointY);
			}
			}  flagsOld=flags;
		}
	}
//Output to log1
print("");  
if (centreAlign==0) { print("Cell displacement calc.: "+name1); }
else {
	print("Transform to stationary: "+name1);
	if (trackingPoints==1) { print("Tracking single point"); }
	else { print("Tracking proximal and distal points"); }
}

// slice array - ARRAY CALLED 'x', also spare arrays 'exy' and 'exx' for plots //
	x=newArray(len);  exy=newArray(5*len);  exx=newArray(5*len);
	for (i=0; i<len; i++) { x[i]=i*sliceStep+1; }
	for (i=0; i<5*len; i++) { exx[i]=(i*sliceStep+0.5)/5; }

// Best fit lines a=parameters for x, b=parameters for y
	a1=newArray(2);  a2=newArray(3);  a3=newArray(4);  a4=newArray(5);
	b1=newArray(2);  b2=newArray(3);  b3=newArray(4);  b4=newArray(5);
	resX=newArray(5);  resY=newArray(5);

// proxX parameters
	coeffs_proxX=newArray(5);  coeffs_proxY=newArray(5);
// setting up matrices
	c0=newArray(5);  c1=newArray(5);  c2=newArray(5);  c3=newArray(5);  c4=newArray(5);  d=newArray(5);
	for (i=0; i<5; i++) {
		c0[i]=0;  c1[i]=0;  c2[i]=0;  c3[i]=0;  c4[i]=0;  d[i]=0;
		for (n=0; n<len; n++) {
			c0[i] += pow(x[n], i);
			c1[i] += pow(x[n], i+1);
			c2[i] += pow(x[n], i+2);
			c3[i] += pow(x[n], i+3);
			c4[i] += pow(x[n], i+4);
			d[i] += proxX[n]*pow(x[n], i);
		}
	}
// 1st order
	Denom=Det2(c0, c1); 
	a1[0]=Det2(d, c1) / Denom;
	a1[1]=Det2(c0, d) / Denom;
// 2nd order
	Denom=Det3(c0, c1, c2); 
	a2[0]=Det3(d, c1, c2) / Denom;
	a2[1]=Det3(c0, d, c2) / Denom;
	a2[2]=Det3(c0, c1, d) / Denom;
// 3rd order
	Denom=Det4(c0, c1, c2, c3); 
	a3[0]=Det4(d, c1, c2, c3) / Denom;
	a3[1]=Det4(c0, d, c2, c3) / Denom;
	a3[2]=Det4(c0, c1, d, c3) / Denom;
	a3[3]=Det4(c0, c1, c2, d) / Denom;
// 4th order
	Denom=Det5(c0, c1, c2, c3, c4); 
	a4[0]=Det5(d, c1, c2, c3, c4) / Denom;
	a4[1]=Det5(c0, d, c2, c3, c4) / Denom;
	a4[2]=Det5(c0, c1, d, c3, c4) / Denom;
	a4[3]=Det5(c0, c1, c2, d, c4) / Denom;
	a4[4]=Det5(c0, c1, c2, c3, d) / Denom;
// residuals
	resX[1]=Residual(proxX, len, a1, 1);
	resX[2]=Residual(proxX, len, a2, 2);
	resX[3]=Residual(proxX, len, a3, 3);
	resX[4]=Residual(proxX, len, a4, 4);
// chose model
	Plot.create("Proximal X", "Slice #", "X");
	Plot.setLimits(0, 2*sliceStep*len, 0, width);
	Plot.setColor("red");  Plot.add("circles", x, proxX);
	for (i=0; i<5*len; i++) { exy[i]=a4[0]+a4[1]*exx[i]+a4[2]*pow(exx[i], 2)+a4[3]*pow(exx[i], 3)+a4[4]*pow(exx[i], 4); }
	Plot.setColor("magenta");  Plot.add("line", exx, exy);  Plot.addText("4th order: residual, "+resX[4], 0.5, 0.4);
	for (i=0; i<5*len; i++) { exy[i]=a3[0]+a3[1]*exx[i]+a3[2]*pow(exx[i], 2)+a3[3]*pow(exx[i], 3); }
	Plot.setColor("cyan");  Plot.add("line", exx, exy);  Plot.addText("3rd order: residual, "+resX[3], 0.5, 0.3);
	for (i=0; i<5*len; i++) { exy[i]=a2[0]+a2[1]*exx[i]+a2[2]*pow(exx[i], 2); }
	Plot.setColor("blue");  Plot.add("line", exx, exy);  Plot.addText("2nd order: residual, "+resX[2], 0.5, 0.2);
	for (i=0; i<5*len; i++) { exy[i]=a1[0]+a1[1]*exx[i]; }
	Plot.setColor("black");  Plot.add("line", exx, exy);  Plot.addText("1st order: residual, "+resX[1], 0.5, 0.1);
	Plot.setColor("red");  Plot.add("circles", x, proxX);
	Plot.show();

	Dialog.create("Choose model");
	Dialog.addNumber("Order:", 1);
	Dialog.show();
	order_proxX=Dialog.getNumber();

	for (i=0; i<=order_proxX; i++) { 
		if (order_proxX==1) { coeffs_proxX[i]=a1[i]; }
		else if (order_proxX==2) { coeffs_proxX[i]=a2[i]; }
		else if (order_proxX==3) { coeffs_proxX[i]=a3[i]; }
		else if (order_proxX==4) { coeffs_proxX[i]=a4[i]; }
	}
	if ( !keepGraphs ) { selectWindow("Proximal X");  close(); }

// proxY parameters
// setting up matrices
	c0=newArray(5);  c1=newArray(5);  c2=newArray(5);  c3=newArray(5);  c4=newArray(5);  d=newArray(5);
	for (i=0; i<5; i++) {
		c0[i]=0;  c1[i]=0;  c2[i]=0;  c3[i]=0;  c4[i]=0;  d[i]=0;
		for (n=0; n<len; n++) {
			c0[i] += pow(x[n], i);
			c1[i] += pow(x[n], i+1);
			c2[i] += pow(x[n], i+2);
			c3[i] += pow(x[n], i+3);
			c4[i] += pow(x[n], i+4);
			d[i] += 	proxY[n]*pow(x[n], i);
		}
	}
// 1st order
	Denom=Det2(c0, c1); 
	b1[0]=Det2(d, c1) / Denom;
	b1[1]=Det2(c0, d) / Denom;
// 2nd order
	Denom=Det3(c0, c1, c2); 
	b2[0]=Det3(d, c1, c2) / Denom;
	b2[1]=Det3(c0, d, c2) / Denom;
	b2[2]=Det3(c0, c1, d) / Denom;
// 3rd order
	Denom=Det4(c0, c1, c2, c3); 
	b3[0]=Det4(d, c1, c2, c3) / Denom;
	b3[1]=Det4(c0, d, c2, c3) / Denom;
	b3[2]=Det4(c0, c1, d, c3) / Denom;
	b3[3]=Det4(c0, c1, c2, d) / Denom;
// 4th order
	Denom=Det5(c0, c1, c2, c3, c4); 
	b4[0]=Det5(d, c1, c2, c3, c4) / Denom;
	b4[1]=Det5(c0, d, c2, c3, c4) / Denom;
	b4[2]=Det5(c0, c1, d, c3, c4) / Denom;
	b4[3]=Det5(c0, c1, c2, d, c4) / Denom;
	b4[4]=Det5(c0, c1, c2, c3, d) / Denom;
// residuals
	resY[1]=Residual(proxY, len, b1, 1);
	resY[2]=Residual(proxY, len, b2, 2);
	resY[3]=Residual(proxY, len, b3, 3);
	resY[4]=Residual(proxY, len, b4, 4);
// chose model
	Plot.create("Proximal Y", "Slice #", "Y");
	Plot.setLimits(0, 2*sliceStep*len, 0, height);
	Plot.setColor("red");  Plot.add("circles", x, proxY);
	for (i=0; i<5*len; i++) { exy[i]=b4[0]+b4[1]*exx[i]+b4[2]*pow(exx[i], 2)+b4[3]*pow(exx[i], 3)+b4[4]*pow(exx[i], 4); }
	Plot.setColor("magenta");  Plot.add("line", exx, exy);  Plot.addText("4th order: residual, "+resY[4], 0.5, 0.4);
	for (i=0; i<5*len; i++) { exy[i]=b3[0]+b3[1]*exx[i]+b3[2]*pow(exx[i], 2)+b3[3]*pow(exx[i], 3); }
	Plot.setColor("cyan");  Plot.add("line", exx, exy);  Plot.addText("3rd order: residual, "+resY[3], 0.5, 0.3);
	for (i=0; i<5*len; i++) { exy[i]=b2[0]+b2[1]*exx[i]+b2[2]*pow(exx[i], 2); }
	Plot.setColor("blue");  Plot.add("line", exx, exy);  Plot.addText("2nd order: residual, "+resY[2], 0.5, 0.2);
	for (i=0; i<5*len; i++) { exy[i]=b1[0]+b1[1]*exx[i]; }
	Plot.setColor("black");  Plot.add("line", exx, exy);  Plot.addText("1st order: residual, "+resY[1], 0.5, 0.1);
	Plot.setColor("red");  Plot.add("circles", x, proxY);
	Plot.show();

	Dialog.create("Choose model");
	Dialog.addNumber("Order:", 1);
	Dialog.show();
	order_proxY=Dialog.getNumber();

	for (i=0; i<=order_proxY; i++) { 
		if (order_proxY==1) { coeffs_proxY[i]=b1[i]; }
		else if (order_proxY==2) { coeffs_proxY[i]=b2[i]; }
		else if (order_proxY==3) { coeffs_proxY[i]=b3[i]; }
		else if (order_proxY==4) { coeffs_proxY[i]=b4[i]; }
	}
	if ( !keepGraphs ) { selectWindow("Proximal Y");  close(); }

if (trackingPoints==2) {
// distX parameters
	coeffs_distX=newArray(5);  coeffs_distY=newArray(5);
// setting up matrices
	c0=newArray(5);  c1=newArray(5);  c2=newArray(5);  c3=newArray(5);  c4=newArray(5);  d=newArray(5);
	for (i=0; i<5; i++) {
		c0[i]=0;  c1[i]=0;  c2[i]=0;  c3[i]=0;  c4[i]=0;  d[i]=0;
		for (n=0; n<len; n++) {
			c0[i] += pow(x[n], i);
			c1[i] += pow(x[n], i+1);
			c2[i] += pow(x[n], i+2);
			c3[i] += pow(x[n], i+3);
			c4[i] += pow(x[n], i+4);
			d[i] += distX[n]*pow(x[n], i);
		}
	}
// 1st order
	Denom=Det2(c0, c1); 
	a1[0]=Det2(d, c1) / Denom;
	a1[1]=Det2(c0, d) / Denom;
// 2nd order
	Denom=Det3(c0, c1, c2); 
	a2[0]=Det3(d, c1, c2) / Denom;
	a2[1]=Det3(c0, d, c2) / Denom;
	a2[2]=Det3(c0, c1, d) / Denom;
// 3rd order
	Denom=Det4(c0, c1, c2, c3); 
	a3[0]=Det4(d, c1, c2, c3) / Denom;
	a3[1]=Det4(c0, d, c2, c3) / Denom;
	a3[2]=Det4(c0, c1, d, c3) / Denom;
	a3[3]=Det4(c0, c1, c2, d) / Denom;
// 4th order
	Denom=Det5(c0, c1, c2, c3, c4); 
	a4[0]=Det5(d, c1, c2, c3, c4) / Denom;
	a4[1]=Det5(c0, d, c2, c3, c4) / Denom;
	a4[2]=Det5(c0, c1, d, c3, c4) / Denom;
	a4[3]=Det5(c0, c1, c2, d, c4) / Denom;
	a4[4]=Det5(c0, c1, c2, c3, d) / Denom;
// residuals
	resX[1]=Residual(distX, len, a1, 1);
	resX[2]=Residual(distX, len, a2, 2);
	resX[3]=Residual(distX, len, a3, 3);
	resX[4]=Residual(distX, len, a4, 4);
// chose model
	Plot.create("Distal X", "Slice #", "X");
	Plot.setLimits(0, 2*sliceStep*len, 0, width);
	Plot.setColor("red");  Plot.add("circles", x, distX);
	for (i=0; i<5*len; i++) { exy[i]=a4[0]+a4[1]*exx[i]+a4[2]*pow(exx[i], 2)+a4[3]*pow(exx[i], 3)+a4[4]*pow(exx[i], 4); }
	Plot.setColor("magenta");  Plot.add("line", exx, exy);  Plot.addText("4th order: residual, "+resX[4], 0.5, 0.4);
	for (i=0; i<5*len; i++) { exy[i]=a3[0]+a3[1]*exx[i]+a3[2]*pow(exx[i], 2)+a3[3]*pow(exx[i], 3); }
	Plot.setColor("cyan");  Plot.add("line", exx, exy);  Plot.addText("3rd order: residual, "+resX[3], 0.5, 0.3);
	for (i=0; i<5*len; i++) { exy[i]=a2[0]+a2[1]*exx[i]+a2[2]*pow(exx[i], 2); }
	Plot.setColor("blue");  Plot.add("line", exx, exy);  Plot.addText("2nd order: residual, "+resX[2], 0.5, 0.2);
	for (i=0; i<5*len; i++) { exy[i]=a1[0]+a1[1]*exx[i]; }
	Plot.setColor("black");  Plot.add("line", exx, exy);  Plot.addText("1st order: residual, "+resX[1], 0.5, 0.1);
	Plot.setColor("red");  Plot.add("circles", x, distX);
	Plot.show();

	Dialog.create("Choose model");
	Dialog.addNumber("Order:", 1);
	Dialog.show();
	order_distX=Dialog.getNumber();

	for (i=0; i<=order_distX; i++) { 
		if (order_distX==1) { coeffs_distX[i]=a1[i]; }
		else if (order_distX==2) { coeffs_distX[i]=a2[i]; }
		else if (order_distX==3) { coeffs_distX[i]=a3[i]; }
		else if (order_distX==4) { coeffs_distX[i]=a4[i]; }
	}
	if ( !keepGraphs ) { selectWindow("Distal X");   close(); }

// distY parameters
// setting up matrices
	c0=newArray(5);  c1=newArray(5);  c2=newArray(5);  c3=newArray(5);  c4=newArray(5);  d=newArray(5);
	for (i=0; i<5; i++) {
		c0[i]=0;  c1[i]=0;  c2[i]=0;  c3[i]=0;  c4[i]=0;  d[i]=0;
		for (n=0; n<len; n++) {
			c0[i] += pow(x[n], i);
			c1[i] += pow(x[n], i+1);
			c2[i] += pow(x[n], i+2);
			c3[i] += pow(x[n], i+3);
			c4[i] += pow(x[n], i+4);
			d[i] += 	distY[n]*pow(x[n], i);
		}
	}
// 1st order
	Denom=Det2(c0, c1); 
	b1[0]=Det2(d, c1) / Denom;
	b1[1]=Det2(c0, d) / Denom;
// 2nd order
	Denom=Det3(c0, c1, c2); 
	b2[0]=Det3(d, c1, c2) / Denom;
	b2[1]=Det3(c0, d, c2) / Denom;
	b2[2]=Det3(c0, c1, d) / Denom;
// 3rd order
	Denom=Det4(c0, c1, c2, c3); 
	b3[0]=Det4(d, c1, c2, c3) / Denom;
	b3[1]=Det4(c0, d, c2, c3) / Denom;
	b3[2]=Det4(c0, c1, d, c3) / Denom;
	b3[3]=Det4(c0, c1, c2, d) / Denom;
// 4th order
	Denom=Det5(c0, c1, c2, c3, c4); 
	b4[0]=Det5(d, c1, c2, c3, c4) / Denom;
	b4[1]=Det5(c0, d, c2, c3, c4) / Denom;
	b4[2]=Det5(c0, c1, d, c3, c4) / Denom;
	b4[3]=Det5(c0, c1, c2, d, c4) / Denom;
	b4[4]=Det5(c0, c1, c2, c3, d) / Denom;
// residuals
	resY[1]=Residual(distY, len, b1, 1);
	resY[2]=Residual(distY, len, b2, 2);
	resY[3]=Residual(distY, len, b3, 3);
	resY[4]=Residual(distY, len, b4, 4);
// chose model
	Plot.create("Distal Y", "Slice #", "Y");
	Plot.setLimits(0, 2*sliceStep*len, 0, height);
	Plot.setColor("red");  Plot.add("circles", x, distY);
	for (i=0; i<5*len; i++) { exy[i]=b4[0]+b4[1]*exx[i]+b4[2]*pow(exx[i], 2)+b4[3]*pow(exx[i], 3)+b4[4]*pow(exx[i], 4); }
	Plot.setColor("magenta");  Plot.add("line", exx, exy);  Plot.addText("4th order: residual, "+resY[4], 0.5, 0.4);
	for (i=0; i<5*len; i++) { exy[i]=b3[0]+b3[1]*exx[i]+b3[2]*pow(exx[i], 2)+b3[3]*pow(exx[i], 3); }
	Plot.setColor("cyan");  Plot.add("line", exx, exy);  Plot.addText("3rd order: residual, "+resY[3], 0.5, 0.3);
	for (i=0; i<5*len; i++) { exy[i]=b2[0]+b2[1]*exx[i]+b2[2]*pow(exx[i], 2); }
	Plot.setColor("blue");  Plot.add("line", exx, exy);  Plot.addText("2nd order: residual, "+resY[2], 0.5, 0.2);
	for (i=0; i<5*len; i++) { exy[i]=b1[0]+b1[1]*exx[i]; }
	Plot.setColor("black");  Plot.add("line", exx, exy);  Plot.addText("1st order: residual, "+resY[1], 0.5, 0.1);
	Plot.setColor("red");  Plot.add("circles", x, distY);
	Plot.show();

	Dialog.create("Choose model");
	Dialog.addNumber("Order:", 1);
	Dialog.show();
	order_distY=Dialog.getNumber();

	for (i=0; i<=order_distY; i++) { 
		if (order_distY==1) { coeffs_distY[i]=b1[i]; }
		else if (order_distY==2) { coeffs_distY[i]=b2[i]; }
		else if (order_distY==3) { coeffs_distY[i]=b3[i]; }
		else if (order_distY==4) { coeffs_distY[i]=b4[i]; }
	}
	if ( !keepGraphs ) { selectWindow("Distal Y");   close(); }
}	// end of distal track if statement!!

//Output to log2
if (trackingPoints==1) { print("Single point:"); }
else { print("Proximal point:"); }
print("  X order="+order_proxX);
for (i=0; i<=order_proxX; i++) { print("    x"+i+"  "+coeffs_proxX[i]); }
print("  Y order="+order_proxY);
for (i=0; i<=order_proxY; i++) { print("    y"+i+"  "+coeffs_proxY[i]); }
if (trackingPoints==2) {
	print("Distal point"); 
	print("  X order="+order_distX);
	for (i=0; i<=order_distX; i++) { print("    x"+i+"  "+coeffs_distX[i]); }
	print("  Y order="+order_distY);
	for (i=0; i<=order_distY; i++) { print("    y"+i+"  "+coeffs_distY[i]); }
}

// TRANSFORMATIONS //
	if (centreAlign!=0) {	// skip if cell displacement calc.
	sureSelect(ID1, name1);	// make sure I'm still in the right image!
	newSize=sqrt(pow(width, 2) + pow(height, 2));
	run("Canvas Size...", "width="+newSize+" height="+newSize+" position=Center");
	for (i=1; i<=NS; i++) {
		x1=0;  for (n=0; n<=order_proxX; n++) { x1 += coeffs_proxX[n]*pow(i, n); }
		y1=0;  for (n=0; n<=order_proxY; n++) { y1 += coeffs_proxY[n]*pow(i, n); }
// Move
		if (centreAlign>1) {
			delx=round(width/2 - x1);  dely=round(height/2 - y1);
			setSlice(i);  run("Select All");  run("Cut");
			makeRectangle(delx, dely, newSize, newSize);  run("Paste");
		}
// Rotate
		if (centreAlign!=2) {
			if (trackingPoints==1) {
				x2=0;    for (n=0; n<=order_proxX; n++) { x2 += coeffs_proxX[n]*pow(i+1, n); }
				y2=0;    for (n=0; n<=order_proxY; n++) { y2 += coeffs_proxY[n]*pow(i+1, n); }
			} else {
				x2=0;    for (n=0; n<=order_distX; n++) { x2 += coeffs_distX[n]*pow(i, n); }
				y2=0;    for (n=0; n<=order_distY; n++) { y2 += coeffs_distY[n]*pow(i, n); }			}
			if (x2-x1==0) { 
				if (y1-y2<0) { angle=90; } else { angle=270; }
			} else { 
				angle=180*atan((y1-y2)/(x2-x1))/pi;
				if (x2-x1<0) { angle += 180; }
			}
			setSlice(i);  run("Select All");
			run("Arbitrarily...", "slice angle="+angle+" interpolate");
		}
	}
	setSlice(1);
	} // end of cell displacement skip

// CALCULATE CELL DISPLACEMENT //
	dist=0;
	for (i=1; i<NS; i++) {
		x1=0;   for (n=0; n<=order_proxX; n++) { x1 += coeffs_proxX[n]*pow(i, n); }
		y1=0;   for (n=0; n<=order_proxY; n++) { y1 += coeffs_proxY[n]*pow(i, n); }		x2=0;   for (n=0; n<=order_proxX; n++) { x2 += coeffs_proxX[n]*pow(i+1, n); }
		y2=0;   for (n=0; n<=order_proxY; n++) { y2 += coeffs_proxY[n]*pow(i+1, n); }		dist += sqrt(pow(x2-x1,2)+pow(y2-y1,2));
	}
	displacement=dist/(NS-1);
// ADD CELL DISPLACEMENT RULER //
	if (doTicks==true && centreAlign!=0) {
		tickN=roundUp((dist+newSize)/tickSpace)+1;  tickX=newArray(tickN);
		for (n=0; n<tickN; n++) {  tickX[n]=n*tickSpace; }
		name=getTitle()+'_ruler';
		run("Duplicate...", "title="+name+" duplicate"); }
	//slice1
		setSlice(1);  Y=tickY+newSize/2;  setLineWidth(pencilWidth);
		for (n=0; n<tickN; n++) {
			if (tickX[n]>=0 && tickX[n]<newSize) { drawLine(tickX[n],Y,tickX[n],Y+tickHeight); } }
	//other slices
		for (i=1; i<NS; i++) {
			setSlice(i+1);  setLineWidth(pencilWidth);
			x1=0;   for (n=0; n<=order_proxX; n++) { x1 += coeffs_proxX[n]*pow(i, n); }
			y1=0;   for (n=0; n<=order_proxY; n++) { y1 += coeffs_proxY[n]*pow(i, n); }			x2=0;   for (n=0; n<=order_proxX; n++) { x2 += coeffs_proxX[n]*pow(i+1, n); }
			y2=0;   for (n=0; n<=order_proxY; n++) { y2 += coeffs_proxY[n]*pow(i+1, n); }
			for (n=0; n<tickN; n++) {  tickX[n]-=sqrt(pow(x2-x1,2)+pow(y2-y1,2)); }
			for (n=0; n<tickN; n++) {
				if (tickX[n]>=0 && tickX[n]<newSize) { drawLine(tickX[n],Y,tickX[n],Y+tickHeight); } }
		}
		setSlice(1);
	} // end of doTicks
//Output to log3
	print("Cell displacement= "+displacement+" pix/frame");
	showMessage("Av. cell displacement= "+displacement+" pix/frame");
}		//End of macro

function Residual(test, len, a, order) {
	Sr=0;
	for (i=0; i<len; i++) {
		x=0;
		for (j=0; j<=order; j++) {
			x += a[j]*pow(i*sliceStep+1, j);
		}	
		Sr += pow( (x - test[i]), 2);
	}
	return Sr;
}		//End of function Residual

macro 'Transform to stationary options...' {
	Dialog.create("Transform to stationary options");
	Dialog.addNumber("Sample every n slices: n:", sliceStep);
	choice=newArray("2", "1");
	Dialog.addChoice("Number of points to track:", choice);
	choice=newArray("Centre+Align", "Centre only", "Align only", "Cell displacement calc.");
	Dialog.addChoice("Transformation:", choice);
	Dialog.addCheckbox("Keep fit graphs:", keepGraphs);
	Dialog.addCheckbox("Add cell displacement ruler:",doTicks); 
	Dialog.addNumber("Position of ruler below centre [pix]:", tickY);
	Dialog.addNumber("Height of ticks [pix]:",tickHeight);
	Dialog.addNumber("Tick spacing [pix]:", tickSpace);
	Dialog.show();
	sliceStep=Dialog.getNumber();
	trackingPoints =Dialog.getChoice();
	choice=Dialog.getChoice();
	if (choice=="Centre+Align") { centreAlign =3; }
	else if (choice=="Centre only") { centreAlign =2; }
	else if (choice=="Align only") { centreAlign =1; }
	else if (choice=="Cell displacement calc.") { centreAlign=0;  trackingPoints=1; }
	keepGraphs=Dialog.getCheckbox();
	doTicks=Dialog.getCheckbox();
	tickY=Dialog.getNumber();
	tickHeight=Dialog.getNumber();
	tickSpace=Dialog.getNumber();
}

// Det method functions //
function Det5(a0, a1, a2, a3, a4) {
	b1=newArray(4);  b2=newArray(4);  b3=newArray(4);  b4=newArray(4);  val=0;
	for (i=0; i<5; i++) {
		j=0;
		for (n=0; n<4; n++) {
			if (n==i) { j++; }
			b1[n]=a1[j];  b2[n]=a2[j];  b3[n]=a3[j];  b4[n]=a4[j];
			j++;
		}
		val += a0[i]*pow(-1, i)*Det4(b1, b2, b3, b4);
	}
	return val; }

function Det4(a0, a1, a2, a3) {
	b1=newArray(3);  b2=newArray(3);  b3=newArray(3);  val=0;
	for (i=0; i<4; i++) {
		j=0;
		for (n=0; n<3; n++) {
			if (n==i) { j++; }
			b1[n]=a1[j];  b2[n]=a2[j];  b3[n]=a3[j];
			j++;
		}
		val += a0[i]*pow(-1, i)*Det3(b1, b2, b3);
	}
	return val; }

function Det3(a0, a1, a2) {
	b1=newArray(2);  b2=newArray(2);  val=0;
	for (i=0; i<3; i++) {
		j=0;
		for (n=0; n<2; n++) {
			if (n==i) { j++; }
			b1[n]=a1[j];  b2[n]=a2[j];
			j++;
		}
		val += a0[i]*pow(-1, i)*Det2(b1, b2);
	}
	return val; }

function Det2(a0, a1) {
	val = a0[0]*a1[1] - a0[1]*a1[0];
	return val; }
// End Det method functions //

////////////////////////////////////////////////////
//  Beat analysis macros
//
//  Written by Bill Wickstead June 2006

macro '         --beat analysis--' {}

function playbackSub() {
	t=1000/beatSpeed;
	for (i=1; i<beatStep; i++) {
		wait(t);
		setSlice(getSliceNumber()+1);
	}
}

macro 'Wavelength and amplitude [w]' {
	NoStackError();
	NS=nSlices;
	showMessageWithCancel("Wavelength and amplitude estimated from wave min/max positions thru stack.\n(assumes wave is aligned along x axis)\nPlease click on min/max for each slice\n \n<shift><click> ends macro and outputs results");

	setTool(7);  setSlice(1);
	continue=1;  count=0;  flagsOld=-1;
	x1=newArray(NS);  y1=newArray(NS);  x2=newArray(NS);  y2=newArray(NS);
	frameN=newArray(NS);
while (continue==1 && count<NS) {		// continue loop
	i=0;
	while (i<2) {
		getCursorLoc(x, y, z, flags);
		if (flags != flagsOld) {
			if (flags == leftButton) {
				if (i<1) {
					x1[count]=x;  y1[count]=y;  i++;
					pointX=newArray(i);  pointY=newArray(i);
					pointX[0]=x1[count];  pointY[0]=y1[count];
					makeSelection("point", pointX, pointY);
				} else {
					x2[count]=x;  y2[count]=y;  i++;
					pointX=newArray(i);  pointY=newArray(i);
					pointX[0]=x1[count];  pointY[0]=y1[count];  pointX[1]=x2[count];  pointY[1]=y2[count];
					makeSelection("point", pointX, pointY);
				}
			} else if (flags == shiftClick) { continue=0;  i=2; 
			} else if (flags == rightButton) { print("exit signal");   exit; }
		}  flagsOld=flags;
	}
	frameN[count]=getSliceNumber();
	if (getSliceNumber()<NS) { setSlice(getSliceNumber()+1); }
		else { continue=0; }
	if (continue==1) { count++; }
}	// end continue loop
	
// OUTPUT
	run("Clear Results");
	allAmp=0;  allLambda=0;
	for (n=0; n<count; n++) {
		if (curvCal==0) { amp=abs(y1[n]-y2[n])/2;  lambda=2*abs(x1[n]-x2[n]); } 
			else { amp=curvCal*abs(y1[n]-y2[n])/2;  lambda=curvCal*2*abs(x1[n]-x2[n]);  }
		setResult("frame #", n, frameN[n]);
		if (curvCal==0) { setResult("amplitude [pix]", n, amp);  setResult("wavelength [pix]", n, lambda);  }
			else { setResult("amplitude [um]", n, amp);  setResult("wavelength [um]", n, lambda);  }
		allAmp+=amp;  allLambda+=lambda; 
	}
	allAmp=allAmp/count;  allLambda=allLambda/count;
	temp='total';
	setResult("frame #", count, temp);
	if (curvCal==0) { setResult("amplitude [pix]", count, allAmp);  setResult("wavelength [pix]", count, allLambda);  }
		else { setResult("amplitude [um]", count, allAmp);  setResult("wavelength [um]", count, allLambda);  }
	updateResults();

	print(' ');  name=getTitle();  print('Wavelength and amplitude measurements for: '+name);
	if (curvCal==0) { print("amplitude: "+allAmp+" [pix]    wavelength: "+allLambda+" [pix]");  }
		else { print("(calibration factor="+curvCal+" um/pix)");  print("amplitude: "+allAmp+" [um]    wavelength: "+allLambda+" [um]");  }
} 	// end wavelength and amplitude macro

macro 'Wavelength and amplitude options...' {
	Dialog.create("Wavelength and amplitude  options");
	Dialog.addNumber("Calibration (um/pix) [0=none]:", curvCal);
	Dialog.show();
	curvCal =Dialog.getNumber();	
}

macro 'Curvature [C]' {
	if (curvClearResults) { run("Clear Results"); }  NR=nResults;
	run("Fit Spline");
	getSelectionCoordinates(x,y);

	len=(x.length)/sampleSpline;
	xx=newArray(len+1);  yy=newArray(len+1);
	for (n=0; n<len; n++) {
		xx[n]=x[sampleSpline*n];  yy[n]= getHeight-y[sampleSpline*n];
	}

	theta=newArray(len);  dist=newArray(len+1);  curv=newArray(len-1);  d=newArray(len-1);
	cumulative=0;
	quad=3*pi/4;  piCount=0;
	for (n=0; n<len-1; n++) {
		theta[n]=2*piCount*pi+atan2( yy[n+1]-yy[n], xx[n+1]-xx[n] );
		if (n>0) {
			if (theta[n]-2*piCount*pi>quad && theta[n-1]-2*piCount*pi<-quad) {
				piCount--;  theta[n]+=2*piCount*pi;  
			} else if (theta[n]-2*piCount*pi<-quad && theta[n-1]-2*piCount*pi>quad) {
				piCount++;  theta[n]+=2*piCount*pi;
			}
		}

		if (curvCal==0) { dist[n]=cumulative; } else { dist[n]=cumulative*curvCal; } 
		cumulative+=sqrt( (yy[n+1]-yy[n])*(yy[n+1]-yy[n]) + (xx[n+1]-xx[n])*(xx[n+1]-xx[n]) );
	}
	if (curvCal==0) { dist[n]=cumulative; } else { dist[n]=cumulative*curvCal; }

	for (n=0; n<len-2; n++) { 
		curv[n]=2*(theta[n]-theta[n+1])/(dist[n+2]-dist[n]);
		if (curvCal==0) { setResult("distance [pix]", n+NR, dist[n+1]); }
			else { setResult("distance [um]", n+NR, dist[n+1]); }
	}

	if (curvSmooth==0) {
		for (n=0; n<len-2; n++) { 
			if (curvCal==0) { setResult("curvature [rad/pix]", n+NR, curv[n]); }
				else { setResult("curvature [rad/um]", n+NR, curv[n]); }
		}
	} else {
		thetaSmooth=newArray(len);
		for (n=0; n<len-2; n++) { 
			if (n>=curvSmooth && len-2-n>=curvSmooth) {
				thetaAve=0;
				for (i=-curvSmooth; i<=curvSmooth; i++) { thetaAve+=curv[n+i]; }
				thetaSmooth[n]=thetaAve/(2*curvSmooth+1);
			} else { 
				thetaSmooth[n]=curv[n];
			}
		}
		for (n=0; n<len-2; n++) {
			curv[n]=thetaSmooth[n];
			if (curvCal==0) { setResult("curvature [rad/pix]", n+NR, curv[n]); }
				else { setResult("curvature [rad/um]", n+NR, curv[n]); }
		}
	}
	updateResults();

	curvmax=0;  curvmin=0;
	for (n=0; n<len-2; n++) { 
		if (curv[n]>curvmax) { curvmax=curv[n]; }
			else if (curv[n]<curvmin) { curvmin=curv[n]; }
		d[n]=dist[n+1];
	}
	if (curvPlot==true) {
		if (curvCal==0) { Plot.create("Curvature", "distance [pix]", "curvature [rad/pix]"); }
			else { Plot.create("Curvature", "distance [um]", "curvature [rad/um]"); }
		if (curvCal==0) { Plot.setLimits(0, cumulative, curvmin, curvmax); }
			else { Plot.setLimits(0, cumulative*curvCal, curvmin, curvmax); }
		Plot.setLineWidth(1);
		Plot.setColor("red");  Plot.add("crosses", d, curv );
		centrelineX=newArray(0, cumulative); centrelineY=newArray(0,0);
		Plot.setColor("black");  Plot.add("line", centrelineX, centrelineY );
		Plot.show;
	}
}

macro 'Shear angle [H]' {
	if (curvClearResults) { run("Clear Results"); }  NR=nResults;
	run("Fit Spline");
	getSelectionCoordinates(x,y);

	len=(x.length)/sampleSpline;
	xx=newArray(len+1);  yy=newArray(len+1);
	for (n=0; n<len; n++) {
		xx[n]=x[sampleSpline*n];  yy[n]= getHeight-y[sampleSpline*n];
	}

	theta=newArray(len);  dist=newArray(len);
	cumulative=0;  thetamax=0;  thetamin=0;
	quad=3*pi/4;  piCount=0;
	for (n=0; n<len-1; n++) {
		theta[n]=2*piCount*pi+atan2( yy[n+1]-yy[n], xx[n+1]-xx[n] );
		if (n>0) {
			if (theta[n]-2*piCount*pi>quad && theta[n-1]-2*piCount*pi<-quad) {
				piCount--;  theta[n]+=2*piCount*pi;  
			} else if (theta[n]-2*piCount*pi<-quad && theta[n-1]-2*piCount*pi>quad) {
				piCount++;  theta[n]+=2*piCount*pi;
			}
		}

		if (theta[n]>thetamax) { thetamax=theta[n]; }
			else if (theta[n]<thetamin) { thetamin=theta[n]; }
		cumulative+=0.5*sqrt( (yy[n+1]-yy[n])*(yy[n+1]-yy[n]) + (xx[n+1]-xx[n])*(xx[n+1]-xx[n]) );
		if (curvCal==0) { dist[n]=cumulative; } else { dist[n]=cumulative*curvCal; } 
		if (curvCal==0) { setResult("distance [pix]", n+NR, dist[n]); }
			else { setResult("distance [um]", n+NR, dist[n]); }
	}
	if (curvSmooth==0) { 
		for (n=0; n<len-1; n++) { setResult("shear [rad]", n+NR, theta[n]); }
	} else {
		thetaSmooth=newArray(len);
		for (n=0; n<len-1; n++) { 
			if (n>=curvSmooth && len-1-n>=curvSmooth) {
				thetaAve=0;
				for (i=-curvSmooth; i<=curvSmooth; i++) { thetaAve+=theta[n+i]; }
				thetaSmooth[n]=thetaAve/(2*curvSmooth+1);
			} else { 
				thetaSmooth[n]=theta[n];
			}
		}
		for (n=0; n<len-1; n++) {
			theta[n]=thetaSmooth[n];
			setResult("shear [rad]", n+NR, theta[n]); 
		}
	}
	updateResults();

	if (curvPlot==true) {
		if (curvCal==0) { Plot.create("Shear angle", "distance [pix]", "angle [rad]"); }
			else { Plot.create("Shear angle", "distance [um]", "angle [rad]"); }
		if (curvCal==0) { Plot.setLimits(0, cumulative, thetamin, thetamax); }
			else { Plot.setLimits(0, cumulative*curvCal, thetamin, thetamax); }
		Plot.setLineWidth(1);
		Plot.setColor("red");  Plot.add("crosses", dist, theta );
		centrelineX=newArray(0, cumulative); centrelineY=newArray(0,0);
		Plot.setColor("black");  Plot.add("line", centrelineX, centrelineY );
		Plot.show;
	}
}

macro 'Curvature/Shear angle options...' {
	Dialog.create("Curvature/Shear angle options");
	Dialog.addNumber("Sample every n-th point:", sampleSpline);
	Dialog.addNumber("Calibration (um/pix) [0=none]:", curvCal);
	Dialog.addNumber("Smoothing (points) [0=none]:", curvSmooth);
	Dialog.addCheckbox("Clear results:", curvClearResults);
	Dialog.addCheckbox("Plot results:", curvPlot);
	Dialog.show;
	
	sampleSpline=Dialog.getNumber();
	curvCal=Dialog.getNumber();
	curvSmooth= Dialog.getNumber();
	curvClearResults=Dialog.getCheckbox();
	curvPlot=Dialog.getCheckbox();
	if (sampleSpline <1 || sampleSpline >50) { showMessage("Out of range"); exit; }
}
//  THE END //