/*
Align_Fluor_Channels 17
Norbert Vischer
25.07.17 14:53

Macros to correct misalignment of fluorescence channels 
with respect to phase contrast in channel #1. We use this for bacteria such as e. Coli.
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.
	Keep Caps lock key down to observe or abort the alignment process.
	For safety, "Offsets.txt" will be renamed to "AppliedOffset.txt" so 
	it will not be used twice.
	
*/

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

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



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);
	setTool(0);
}

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.51p");
	setBatchMode(true);
	print("\\Clear");
	run("Clear Results");
	while (isOpen("Offsets.txt")){
		selectWindow("Offsets.txt");
		run("Close");
	}
	currentDir = getDirectory("Choose directory containing misaligned stacks");
	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;
	row = 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;
	firstTime = true;
	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 = getResultString("file", row);//25.10.2015
		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();
			caps = is("Caps Lock Set");
			resetMinAndMax();//25.10.2015
			msg ="Translate dx="+ dx + " dy=" + dy + " ?\n ";
			inspect(msg, firstTime);
			firstTime = false;
			run("Select None");
			run("Translate...", "x=" + dx +" y=" + dy +" interpolation=Bilinear slice");
			inspect(msg, firstTime);
		}
	}
	if (getBoolean("Save and Close?")){
			run("Save");
			close;
		}
	File.rename(offsetPath, currentDir + "AppliedOffset.txt");
}

//queries user to continue depending on caps lock and first time
function inspect(msg, firstTime){
	caps = is("Caps Lock Set");
	if (firstTime || (caps)){
		msg = msg + "\nKeep Caps Lock key down to remain in Inspection Mode";
		msg = msg + "\nUse Esc to abort";
		msg = msg + "\n---\nClick OK to Continue";
		waitForUser("Inspection Mode", msg);
	}
}