/*
Align_Fluor_Channels 15
Norbert Vischer
18-jul-2014  19:18

Macros to correct misaligned fluorescent channels in a hyperstack with respect to
first channel that contains phase contrast bacteria.
You can choose manual or semi-automatic methods. 
Make sure to have a back-up of your images in case alignment has erratic behaviour.
These macros work without ObjectJ.


a) Manual Alignment: 
  - open hyperstack to be aligned
  - select a fluorescence channel
  - choose "Show PhC Outlines [1]" to show outlines derived from channel 1,
  - preferably zoom in
  - click with a selection tool inside one of the contours and drag outlines to correct position 
  - choose "Apply Manual Alignment [2]" to make alignment permanent
  - repeat this for all fluorescence channels, then save the aligned stack

b) Automatic Alignment:

- doubling fftSize will quadruple processing time and double precision
- choose "Create Offsets Table... [5]"
	this asks to select a folder in order to align all tiff files found there.
	A text file "Offsets.txt" will be saved  in the same folder
	
- choose "Apply Offsets Table [6]"
	asks to select a folder containing file "Offsets.txt".
	This file will be opened and each listed frame is adjusted after user confirmation.
	For safety, "Offsets.txt" will be renamed to "AppliedOffset.txt" so 
	it will not be used twice.

c) Fixed Alignment:
- This is not implemented yet. If (after a test phase) it appears 
that same channels have always the same misalignment, it may 
be more safe and efficient to find a mean offset per channel and 
apply this constant value to that channel in all images
	
*/

var 
selX, selY, selW, selH,//remember manual outlines

fftSize = 4096, //choose from 1024, 2048, 4096, ...
row = 0,
currentDir,//current directory
k40 = 40,

;



macro "Show PhC Outlines [1]"{
	Stack.getPosition(channel, slice, frame);
	Stack.setPosition(1, slice, frame);
	run("Select None");
	setAutoThreshold("Default");

	run("Create Selection");
	roiManager("reset");
	roiManager("add");
	getSelectionBounds(selX, selY, selW, selH);
	resetThreshold();
	Stack.setPosition(channel, slice, frame);
}

function showOutlinesOverlay(){
	Stack.getPosition(channel, slice, frame);
	Stack.setPosition(1, slice, frame);
	run("Select None");
	setAutoThreshold("Default");


	run("Create Selection");
	getSelectionBounds(selX, selY, selW, selH);
	resetThreshold();
	Stack.setPosition(channel, slice, frame);
	run("Remove Overlay");
	run("Add Selection...");
	run("Select None");
}

macro "Apply Manual Alignment [2]"{
	getSelectionBounds(left, top, ww, hh);
	if (ww != selW || hh != selH || roiManager("count") != 1)
		exit("You need to refresh PhC outlines");
	Stack.getPosition(channel, slice, frame);
	if (channel == 1)
		exit("channel 1 cannot be aligned");
	dx = selX - left;
	dy = selY - top;
	run("Select All");
	run("Translate...", "x=" + dx +" y=" + dy +" interpolation=Bilinear slice");
	roiManager("select", 0);
	Stack.setPosition(channel, slice, frame);
}

macro "Remove Outlines [3]"{
	run("Select None");
	run("Remove Overlay");
}


//checks all tiff image in folder and records offsets of mis-aligned channels
//Images are not changed yet, and
//results are stored in  Offsets.txt in same folder 

macro "Create Offsets Table... [5]"{
	requires("1.48m");
	setBatchMode(true);
	print("\\Clear");
	run("Clear Results");
	while (isOpen("Offsets.txt")){
		selectWindow("Offsets.txt");
		run("Close");
	}
	currentDir = getDirectory("");
	allFiles = getFileList(currentDir);
	files = newArray(allFiles.length);
	nTiffs = 0;
	print("directory:  ", currentDir);
	print("---");
	
	for (jj = 0; jj < allFiles.length; jj++){
		name = allFiles[jj];
		if (endsWith(name, ".tif") || endsWith(name, ".tiff")){
			files[nTiffs++] = name;
			print(name);
		}
	}
	files = Array.trim(files, nTiffs);
	selectWindow("Log");
	waitForUser("Files to be processed are shown in Log Window \nClick OK to continue or Esc to stop");
 	found = 0;
	run("Close All");
	newImage("Observer", "8-bit white", k40, k40, 1);
	makeLine(0, k40/2, k40, k40/2);
	run("Add Selection...");
	makeLine(k40/2, 0, k40/2, k40);
	run("Add Selection...");
	setLocation(6, 24);
	setBatchMode("show");
	run("Set... ", "zoom=400 x=20 y=20");


	img = 0;
	imgCount = files.length;
	for (jj = 0; jj < files.length; jj++){
		name = files[jj];
		fullPath = currentDir+name;
		open(fullPath);
		Stack.getDimensions(ww, hh, channels, slices, frames);
		ok = (channels > 1 && !(slices > 1 && frames > 1));
		if (!ok) close;
		else {
			img++;
			calcDisplacement(img, imgCount);
			
			selectImage("Observer");
			close("\\Others");
			found++;
		}
	}		
	showMessage("Number of files processed: " + found + " of " + files.length);	
	selectWindow("Results");
	run("Input/Output...", "jpeg=75 gif=-1 file=.txt save_column");

	saveAs("Results", currentDir + "Offsets.txt");
}


function calcDisplacement(thisImg, imgCount){
	run("Line Width...", "line=1");
	enough = false;
	checkedFrames = 0;
	hStackID = getImageID;//hyperstackID
	title = getTitle;

	selectImage(hStackID);
	run("Select None");
	xFactor = getWidth/fftSize;
	yFactor = getHeight/fftSize;
	
	Stack.getDimensions(width, height, nChanns, slices1, frames1);
	if (nChanns == 1) exit("Hyperstack with at least 2 channels expected");
	if (slices1 !=1 && frames1 !=1) exit("Either nSlices or nFrames must be =1");
	if (slices1 > 1)	
		run("Stack to Hyperstack...", "order=xyczt(default) channels=&nChanns slices=1 frames=&slices1 display=Color");
	Stack.getDimensions(width, height, nChanns, slices1, frames1);
	Stack.setDisplayMode("color");
	selectImage(hStackID);

	chNames = newArray(nChanns+1);
	for (ch = 1; ch <= nChanns; ch++)//one-based
		chNames[ch] = "Ch_" + ch;		 
	for (thisFrame = 1; thisFrame <= frames1 && !enough; thisFrame++){
		showProgress(thisFrame, frames1);
		//allOffsets += "" + thisFrame;
		for (chn = 1; chn <= nChanns; chn++){
			selectImage(hStackID);
			Stack.setPosition(chn, 1, thisFrame);
			run("Duplicate...", "title=Tmp slice");	
			run("8-bit");
			run("Scale...", "x=- y=- width=" + fftSize + " height=" + fftSize + " interpolation=Bilinear average create title=" + chNames[chn]);
		}
		close("Tmp");
		

		// now we have e.g. 3 channels called Ch_1, Ch_2, Ch_3, all are squares of fftSize
		name = chNames[1];
		selectImage(name);//convert phase contrast channel into mask
		setAutoThreshold("Default");
		run("Convert to Mask");
		run("Grays");
		for (channel = 2; channel <= nChanns; channel++){

			leftName = chNames[1];
			rightName = chNames[channel];
			run("FD Math...", "image1=" + leftName +" operation=Correlate image2="+ rightName + " result=Result do");	
			run("8-bit");
			halfFFT = fftSize/2;
			showProgress(img + thisFrame/frames1, imgCount);
			makeLine(halfFFT - 10, halfFFT, halfFFT + 10, halfFFT);
			run("Add Selection...");
			makeLine(halfFFT, halfFFT - 10, halfFFT, halfFFT + 10);
			run("Add Selection...");
			setLocation(10, 80, 200, 200);
			run("Set... ", "zoom=400 x=&halfFFT y=&halfFFT");

			//setBatchMode("show");
			close("tmp");
			rename("tmp");
			makeRectangle(fftSize * 0.4, fftSize * 0.4, fftSize * 0.2, fftSize * 0.2);
			run("Find Maxima...", "noise=100 output=[Point Selection]");
			getSelectionCoordinates(xx, yy);

			makeRectangle(halfFFT - k40/2, halfFFT - k40/2, k40, k40);
			run("Copy");
			selectImage("Observer");
			run("Paste");
			run("Enhance Contrast", "saturated=0.35");
			makePoint(xx[0] - halfFFT + k40/2, yy[0] - halfFFT + k40/2);
			

			dx = round(10 * xFactor * (xx[0] -fftSize/2))/10 ;//precision 0.1px is enough
			dy = round(10 * yFactor * (yy[0] -fftSize/2))/10;
	
			setResult("file", row, title);
			setResult("frm", row, thisFrame);
			setResult("dx_chn"+ channel, row, dx);
			setResult("dy_chn"+ channel, row, dy);
			updateResults;
	
			close(rightName);
			close("Result");
		}
		row++;
		
		checkedFrames++;
		//allOffsets += "\n";
		close(leftName);
		//close;
	}
}



macro "Apply Offsets Table [6]"{
//debug;
	run("Close All");
		currentDir = getDirectory("Choose a Directory");
	offsetPath =  currentDir + "Offsets.txt";
	if (!File.exists(offsetPath)) 
		exit("Offsets.txt not found");
	run("Results... ", "open=[" + offsetPath +"]");
	selectWindow("Results");
	prevName = "";
	for (row = 0; row < nResults; row++){
		fname = getResultLabel(row);
		if (fname != prevName && prevName!= ""){
			if (getBoolean("Save and Close?")){
				run("Save");
				close;
			}
			
		}
		if (!isOpen(fname)){
			open(currentDir+fname);
			getDimensions(ww, hh, channels, slices, frames);
			run("Set... ", "zoom=200");
			prevName = fname;
		}
		frame = getResult("frm", row);
		for (chn = 2; chn <= channels; chn++){
			Stack.setPosition(chn, 1, frame);
			columnX = "dx_chn" + chn;
			columnY = "dy_chn" + chn;
			dx = getResult(columnX, row);
			dy = getResult(columnY, row);
			showOutlinesOverlay();
			waitForUser("Translate dx="+ dx + " dy=" + dy + " ?");
			run("Select None");
			run("Translate...", "x=" + dx +" y=" + dy +" interpolation=Bilinear slice");
			waitForUser;
		}
	}
	if (getBoolean("Save and Close?")){
			run("Save");
			close;
		}
	File.rename(offsetPath, currentDir + "AppliedOffset.txt");
}