////////////////////////////////////////////////////
//  ImageJ macros for non-linear LSM style correction of ellipical deformation.
//  Note: function assumes regular nonagon deformed by 1 (only) symmetrical ellipse.
//  Any other deformation: forget it!
//
//  Largely replaces functions scripted in NIH.
//
//  Scripted by Bill Wickstead

////////////////////////////////////////////////////
//   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;

// variables associated with ellipse correction / polar transformation
var T=0.000001;
var axoScale=500;
var thetaStep=pi/30;
var alphaStep=pi/30;
var minHit=4;
var itMax=50;
var firstUse=1;
var cropImage=true;
var cropSize=800;

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;
}

////////////////////////////////////////////////////
//  Spectral addition tools
//
//  Written by Bill Wickstead Jan 2006

macro '         --spectral addition--' {}

macro 'Nine-fold rotation [n]' {
	run("Select All");
	run("Copy");
	if (getWidth()>getHeight()) { size=getWidth(); } else { size=getHeight(); }
	type=PictureType();
 	newImage('temp', type, size, size, 9);
	for (i=1; i<=9; i++) { setSlice(i); run("Paste"); }
	for (i=2; i<=9; i++) {
		setSlice(i);
		angle=(i-1)*40;
		run("Arbitrarily...", "slice angle="+angle+" interpolate");
	}
	oldPic=getImageID();
	run("Z Project...", "start=1 stop=9 projection=[Average Intensity]");
	selectImage(oldPic);
	close;
}		// End of 'Nine-fold rotation' macro


macro 'Ellipse correction [e]' {
// Capture 9-fold points
//	if (firstUse==1) { 
//		showMessage('Select 9 points for ellipse fitting\n(alt+click to Undo a point)');  firstUse=0; }
	setTool(7);	
	x=newArray(9); y=newArray(9);	// set arrays
	flagsOld=-1;
	i=0; while (i<9) {
		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]=x[n];  pointY[n]=getHeight-1-y[n]; }
					makeSelection("point", pointX, pointY);
				}
			} else if (flags == leftButton) {
				x[i]=x1;  y[i]=getHeight-1-y1;  i++;
				pointX=newArray(i);  pointY=newArray(i);
				for (n=0; n<i; n++) { pointX[n]=x[n];  pointY[n]=getHeight-1-y[n]; }
				makeSelection("point", pointX, pointY);
			}
		} flagsOld=flags;
	}

// initial values
	tempArray=split(getTitle(), '.');
	name=tempArray[0];
	x0=0;  y0=0;
	for (i=0; i<9; i++) { x0=x0+x[i];  y0=y0+y[i]; }
	x0=x0/9;  y0=y0/9;

	if (y[0]-y0>0) { theta=atan((x[0]-x0)/(y[0]-y0)); }
	else if (y[0]-y0<0) { theta=pi+atan((x[0]-x0)/(y[0]-y0)); }
	else { if (x[0]-x0)>0) { theta=pi/2; } else { theta=3*pi/2; } }

	maxN=0;  minN=0;
	for (i=0; i<9; i++) {
		d=sqrt(pow(x[i]-x0, 2) + pow(y[i]-y0, 2));
		if (i==0) { a=d;  b=d; }
		else { 
			if (d>b) { b=d;  maxN=i; }
			if (d<a) { a=d;  minN=i; }
		}
	}
	if (y[maxN]-y0>0) { alpha=atan((x[maxN]-x0)/(y[maxN]-y0)); }
	else if (y[maxN]-y0<0) { alpha=pi+atan((x[maxN]-x0)/(y[maxN]-y0)); }
	else { if (x[maxN]-x0)>0) { alpha=pi/2; } else { theta=3*pi/2; } }

	aStep=a/20;  bStep=b/20;  xStep=a/500;  yStep=xStep; 
	alStep=alphaStep;  thStep=thetaStep;			// place permenant options into temp holders

// output to Log window
	print("Ellipse Correction: ", name);
	print("Datapoints:");
	for (i=0; i<9; i++) { print('   ',x[i], '   ', y[i]); }
	print('\nInitial values:');
	print('x0: ', x0, '   y0: ', y0, '\ntheta: ', theta, '(', theta*180/pi, ')   alpha: ', alpha, '(', alpha*180/pi, ')');
	print('a: ', a, ' (dblt: ', minN+1, ')', '   b: ', b, ' (dblt: ', maxN+1, ')');
	print('\nFitting...');

// fitting procedure
	SrOld=Nono(x0, y0, theta, alpha, a, b, x, y);
	hitCount=0;
	iterationCount=1;
	while (hitCount<minHit) {
		SrInt=SrOld;
	//x0
		better=0;
		for (n=-2; n<=2; n++) {
			test=x0+n*xStep;
			Sr=Nono(test, y0, theta, alpha, a, b, x, y);
			if (Sr<SrInt) {  better=1;  SrInt=Sr;  betterVal=test; }
		}
		if (better==1) { x0=betterVal; } else { xStep=xStep/2; }
	//y0
		better=0;
		for (n=-2; n<=2; n++) {
			test=y0+n*yStep;
			Sr=Nono(x0, test, theta, alpha, a, b, x, y);
			if (Sr<SrInt) {  better=1;  SrInt=Sr;  betterVal=test; }
		}
		if (better==1) { y0=betterVal; } else { yStep=yStep/2; }
	//a
		better=0;
		for (n=-3; n<=3; n++) {
			test=a+n*aStep;
			Sr=Nono(x0, y0, theta, alpha, test, b, x, y);
			if (Sr<SrInt) {  better=1;  SrInt=Sr;  betterVal=test; } 
		}
		if (better==1) { a=betterVal; } else { aStep=aStep/2; }
	//b
		better=0;
		for (n=-3; n<=3; n++) {
			test=b+n*bStep;
			Sr=Nono(x0, y0, theta, alpha, a, test, x, y);
			if (Sr<SrInt) {  better=1;  SrInt=Sr;  betterVal=test; }
		}
		if (better==1) { b=betterVal; } else { bStep=bStep/2; }
	//theta
		better=0;
		for (n=-3; n<=3; n++) {
			test=theta+n*thStep;
			Sr=Nono(x0, y0, test, alpha, a, b, x, y);
			if (Sr<SrInt) {  better=1;  SrInt=Sr;  betterVal=test; } else { thStep=thStep/2; }
		}
		if (better==1) { theta=betterVal; }
	//alpha
		better=0;
		for (n=-15; n<=15; n++) {
			test=alpha+n*alStep;
			Sr=Nono(x0, y0, theta, test, a, b, x, y);
			if (Sr<SrInt) {  better=1;  SrInt=Sr;  betterVal=test; }
		}
		if (better==1) { alpha=betterVal; } else { alStep=alStep/2; }
	//prepare for next iteration
		print(iterationCount, ':   start: ', SrOld, '    stop: ', SrInt);
		if ( (SrOld-SrInt)/SrOld < T ) { hitCount++; } else { hitCount=0; }
		SrOld=SrInt;
		iterationCount++;
		if (iterationCount>itMax) { exit("No convergence within itMax"); }
	}
	//end parameters output
	if (b<a) { test=b;  b=a;  a=test; alpha=alpha+pi/2; }
	print('\nFinal values:');
	print('x0: ', x0, '   y0: ', y0, '\ntheta: ', theta, '(', theta*180/pi, ')   alpha: ', alpha, '(', alpha*180/pi, ')');
	print('a: ', a, '   b: ', b);

// do tranformation
	//centre
	print('\nTransforms:\n    centring...');
	run("Select All");
	run("Copy");
	setForegroundColor(255,255,255);
	run("Fill");
	makeRectangle(round(-x0+(getWidth()-1)/2), round( y0-(getHeight()-1)/2), getWidth(), getHeight());
	run("Paste");

	//rotate to put ellipse long axis along x
	print('    rotate 1...');
	angle=90-180*alpha/pi;
	run("Select All");
	run("Arbitrarily...", "angle="+angle+" interpolate");

	//correct for skew
	print('    skew correction...');
	ratio=a/b;
	run("Scale...", "x="+ratio+" y=1 interpolate");

	//rotate such that dblt1 is at top
	print('    rotate 2...');
	angle=-90+180*(alpha-theta)/pi;
	run("Select All");
	run("Arbitrarily...", "angle="+angle+" interpolate");

	//normalise
	print('    normalise...');
	ratio=axoScale*0.5/a;
	oldPic=getImageID();
	name=name+'_ellipse_corr';
	run("Scale...", "x="+ratio+" y="+ratio+" interpolate create title=["+name+"]");
	selectImage(oldPic);
	close;

	//crop
	if (cropImage) {
		print('    cropping...');
		left=(getWidth() - cropSize)/2;
		top=(getHeight() - cropSize)/2;
		makeRectangle(left, top, cropSize, cropSize);
		run("Crop");
	}
	print('\nFinished!');
}		//End of 'Ellipse correction' macro

function Nono(x0, y0, theta, alpha, a, b, x, y) {
	Sr=0;
	nonoX=newArray(9);  nonoY=newArray(9);
	for (i=0; i<9; i++) {
		f=2*pi*i/9 + theta + pi/2 - alpha;
		nonoX[i]= b*sin(f)*cos(alpha - pi/2) + a*cos(f)*sin(alpha - pi/2) + x0;
		nonoY[i]= a*cos(f)*cos(alpha - pi/2) - b*sin(f)*sin(alpha - pi/2) + y0;
		Sr=Sr + pow(x[i]-nonoX[i], 2) + pow(y[i]-nonoY[i], 2);
	}
	return Sr;
}		//End of function

macro 'Ellipse correction options...' {
	Dialog.create("Ellipse correction options");
	Dialog.addNumber("Diameter of normalised axonemes [pix]: ", axoScale);
	Dialog.addNumber("Threshold: ", T);
	Dialog.addNumber("Initial stepsize for theta: ", thetaStep);
	Dialog.addNumber("Initial stepsize for alpha: ", alphaStep); 
	Dialog.addNumber("Iterations at threshold for convergence: ", minHit);
	Dialog.addNumber("Quit if no convergence within [iterations]: ", itMax);
	Dialog.addCheckbox("Crop after transform ", cropImage);
	Dialog.addNumber("Crop to size [pix]: ", cropSize);
	Dialog.show();
	axoScale=Dialog.getNumber();
	T=Dialog.getNumber();
	thetaStep=Dialog.getNumber();
	alphaStep=Dialog.getNumber();
	minHit=Dialog.getNumber();
	itMax=Dialog.getNumber();
	cropImage=Dialog.getCheckbox();
	cropSize=Dialog.getNumber();
}		//End of Macro