/*Coli-Correlation.txt
version: 04.10.15 6:35
works in combination with Coli-Inspector
needs plugin  Roi_Correlator.class

For Synchronizing Overlaps.tif and Pearson Plots.tif, choose:
	Analyze>Tools>Synchronize Windows
*/

macro "Correlate Coli"{
	firstCell = 1;
	lastCell = -1;//set to -1 to evaluate all cells
	entireCell = true;//otherwise wall only
	showSumOfPlots = false;
	psfExtend = 0.1;//default: 0.1 um to fetch light outside boundary
	wallThickness = 0.2;//microns
	fraction = 0.3;//default: 0.3 = strongest 30% area of fluorescence per channel
	channelX = 2;
	channelY = 3;
	bins = 100;//for pearson plot
	ojRequires("1.03v2");
	run("Close All");
	batch = true;

	ojHideResults();
	ojInitColumn("Pearson");	
	ojInitColumn("Overlap");
	time = getTime;
	nCells = ojNObjects();
	
	setBatchMode(batch);

	newImage("Overlaps.tif", "8-bit black", 140, 140, ojNObjects());
	overlapsID = getImageID;
	run("Set...", "value=32 stack");//outside cell

	getLut(reds, green, blues);
	reds[32] = 66; green[32] = 66; blues[32] = 66;
	reds[64] = 100; green[64] = 200; blues[64] = 0;
	reds[128] = 0; green[128] = 200; blues[128] = 200;
	reds[192] = 255; green[192] = 0; blues[192] = 0;
	setLut(reds, green, blues);

	newImage("OnePlot", "32-bit black", bins, bins, 1);
	onePlotId = getImageID;
	newImage("Pearson Plots.tif", "32-bit black", bins, bins, nCells);
	pearsonPlotsID = getImageID;
	
	prevID = 0;
	if (lastCell <= 0)
		lastCell = ojNObjects();
	
	for (obj = firstCell; obj <= lastCell; obj++){	
		ojShowObject(obj);
		srcID = getImageID;
		if (srcID != prevID && isOpen(prevID)){
			selectImage(prevID);
			run("Close");
			//beep;
			//waitForUser("closing");
			selectImage(srcID);
			title = getTitle;
		}
		prevID = srcID;
		thr = ojResult("Thr", obj);
		setThreshold(0, thr);
		xstart = ojResult("startX", obj);
		ystart = ojResult("startY", obj);
		doWand(xstart, ystart);

		if (entireCell){
			run("Enlarge...", "enlarge=" + psfExtend);
		}
		else{
			run("Enlarge...", "enlarge=" + (psfExtend -wallThickness));
			run("Make Band...", "band=" + wallThickness);
		}
		doOverlap(obj);


		run("Roi Correlator", "channelx=&channelX channely=&channelY  use=&onePlotId create bins=&bins maxx=0 maxy=0");
		pearson = parseFloat(call("Roi_Correlator.getPearson")); 
		ojSetResult("Pearson", obj, pearson);	
		selectImage(onePlotId);
		run("Select All");
		run("Copy");
		changeValues(0, 1e9, 0);				
		selectImage(pearsonPlotsID);
		setSlice(obj);			
		run("Paste");
		setMetadata("label", "pearson=" + d2s(pearson, 3));			
		showProgress(obj/ojNObjects());
	
	}
	selectImage(srcID);//last linked image
	close;
	close("OnePlot");
	setBatchMode("exit and display");
	if(showSumOfPlots){
		selectImage(pearsonPlotsID);
		run("Z Project...", "projection=[Sum Slices]");
		run("Set... ", "zoom=300");
		run("In [+]");
		run("Enhance Contrast", "saturated=0.35");
	}
	
	selectImage(pearsonPlotsID);
	setSlice(1);	
	run("Enhance Contrast", "saturated=0.35");
	run("Set... ", "zoom=300");
	run("In [+]");

	selectImage(overlapsID);
	run("Set... ", "zoom=300");
	run("In [+]");

	
	ojShowResults();
	showStatus("" + ojNObjects() + " objects correlated in " + round ((getTime - time)) + " ms");
	//run("Synchronize Windows");
	run("Tile");
}


//adds an 8-bit image with 5 colors to stack of Overlaps:
//Example: fraction = 0.3
//both channels are thresholded so that 30% or the roi area is above threshold.
//BLACK:	both below thr
//GREEN:	only channel A above thr
//CYAN:		only channel B above thr
//RED:  	both channels above thr
//GRAY:		outside roi

function doOverlap(obj){
	srcID = getImageID;
	Stack.getDimensions(ww, hh, channels, slices, frames);
	Stack.getPosition(channel, slice, frame);
//	if(obj == 5) waitForUser;
	run("Duplicate...", "duplicate channels=2-3 slices=&slice frames=&frame");//exclude phC
	//setMetadata("info", "");
	rename("SingleCell");
	run("Set... ", "zoom=400");
	run("Select None");
	for (slc = 1; slc <=2; slc++){//scale down to max 255
		setSlice(slc);
		getRawStatistics(nPixels, mean, min, max);
		factor = 255/max;
		run("Multiply...", "value=&factor slice");
		setMinAndMax(0, 255);
	}
	run("8-bit");
	
	run("Restore Selection");
	run("Make Inverse");
	run("Set...", "value=16 stack");//sum outside roi will be 32
	run("Make Inverse");

	setSlice(1);
	setThresholdToFraction(fraction);
	getThreshold(lower, upper);
	changeValues(0, lower, 0);
	changeValues(lower, 1e9, 64);

	setSlice(2);

	setThresholdToFraction(fraction);
	getThreshold(lower, upper);
	changeValues(0, lower, 0);
	changeValues(lower, 1e9, 128);


	run("Z Project...", "projection=[Sum Slices]");
	run("Set... ", "zoom=400");
	setMinAndMax(0, 255);
	run("8-bit");

	run("Restore Selection");//cell outline
	getRawStatistics(nPixels, mean, min, max, std, histo);
	overlap = histo[192];// red pixels
	relOverlap = overlap/(nPixels * fraction);
	ojSetResult("Overlap", obj, relOverlap);
	run("Select All");
	run("Copy");
	close;//SUM

	selectImage("Overlaps.tif");
	setSlice(obj);
	run("Select None");
	run("Paste");
	setMetadata("label", "overlap=" + d2s(relOverlap, 3));			

	run("Select None");
	close("SingleCell");
	selectImage(srcID);
}
	
// set threshold to get fixed fraction	
function setThresholdToFraction(fraction){
	getRawStatistics(nPixels, mean, min, max, std, histo);
	sum = 0;
	thr = 0;
	for (jj= 0; jj < histo.length; jj++){
		sum += histo[jj];
		if (sum > nPixels * (1-fraction)){
			thr = jj;
			break;
		}
	}
	setThreshold(thr, 1e9);
}