//Coli-Correlation.txt
//version: 30.09.15 1:24  vischer@uva.nl
//works in combination with Coli-Inspector
//Measures Pearson's coeff. per cell and overlap per cell

macro "Correlate Coli"{
	
	doCells = 0;//stop measuring after 'doCells' (use 0 to measure all cells)
	psfExtend = 0.1;//default: 0.1 um to fech light outside boundary
	entireCell = true;//otherwise wall only
	wallThickness = 0.2;//microns
	fraction = 0.3;//default: 0.3 = strongest 30% area of fluorescence per channel
	channelX = 2;//fluorescence channel
	channelY = 3;//fluorescence channel
	bins = 100;//for pearson plot
	showSumPlot = false; //SUM projection of pearson plot
	
	
	
	
	run("Close All");
	batch = true;
	withPlot = true;
	withPearson = true;
	withOverlap= true;
	
	ojHideResults();
	ojDeleteColumn("Overlap");
	ojDeleteColumn("Pearson");
	if (withPearson){
		ojInitColumn("Pearson");
		
	}		
	close("*Correlation Plots*");
	close("OnePlot");
	ojShowImage(1);

	time = getTime;
	
	if (doCells <= 0 || doCells > ojNObjects())
		doCells = ojNObjects();
	if (withOverlap){//prepare
		ojInitColumn("Overlap");
		close("Overlaps.tif");
		newImage("Overlaps.tif", "8-bit black", 140, 140, doCells);
		overlapsID = getImageID;
		run("Set...", "value=32 stack");//outside cell

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

		run("Set... ", "zoom=400");
	}

	setBatchMode(batch);
	onePlotId = 0;
	
	newImage("OnePlot", "32-bit black", bins, bins, 1);
	onePlotId = getImageID;
	newImage("Correlation Plots", "32-bit black", bins, bins, doCells);
	allPlotsID = getImageID;
	

	for (img = 1; img <= ojNImages(); img++){	
		if (doCells >= ojFirstObject(img)){
			ojShowImage(img);
			Stack.getDimensions(ww, hh, channels, slices, frames);
			if(channels < channelX || channels < channelX)
				exit("Stack has not enough channels");
			srcID = getImageID;
			for (obj = ojFirstObject(img); obj <= ojLastObject(img) && obj <= doCells; obj++){
				thr = ojResult("Thr", obj);
				ojShowObject(obj);
				selectImage(srcID);
				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);

				if (withPearson){
					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(allPlotsID);
					setSlice(obj);			
					run("Paste");
					setMetadata("label", "pearson=" + d2s(pearson, 3));			
					showProgress(obj/doCells);
				}
			}
			selectImage(srcID);
			close;
		}
	}
	if (withPearson){
		selectImage(allPlotsID);
		setSlice(1);	
		setBatchMode("show");
		selectImage(overlapsID);
		setSlice(1);	
		setBatchMode("show");	
		setBatchMode(false);	
		selectImage(allPlotsID);
		run("Set... ", "zoom=400");
		run("Enhance Contrast", "saturated=0.35");
		if(showSumPlot){
			run("Z Project...", "projection=[Sum Slices]");
			run("Set... ", "zoom=300");
			run("In [+]");
			run("Enhance Contrast", "saturated=0.35");
		}
		//updateDisplay;
	}
	//setBatchMode(false);	
	//selectImage(onePlotId);
	//run("Close");
	ojShowResults();
	showStatus("" + doCells + " objects correlated in " + round ((getTime - time)) + " ms");
	if (withPearson && withOverlap)
		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);
}