//This macro counts the number of above threshold pixels in channel #1 
//("GPR161Channel") that colocalize with above-threshold pixels in channel 
//#2. 

//The macro will ask for the save directory, the threshold value for
//channels #1 and 2, and the user needs to identify the points of interest
//(cilia) in each image by clicking on them.  

//Along with % overlap, to help verify data integrity it will save a RGB 
//of each image with each cilium highlighted and numbered.  It will also 
//save a small RGB of each cilium in the output folder for making grids
//using a separate macro.

//For this macro you need a set of open TIFFs with composite channels, in
//16-bit or 32-bit format.  In order for the threshold values to be 
//valid each image should have been acquired with the same acquisition 
//settings.  


//1. Get the save directory, channels and GPR161 threshold values

dir = getDirectory("Save the results");

Dialog.create("GPR161 quant options");
     Dialog.addChoice("Cilia:", newArray("1", "2", "3", "4")); 
     Dialog.addChoice("GPR161:", newArray("2", "1", "3", "4"));
     Dialog.addNumber("cilia multiplier:", 4);
     Dialog.addNumber("GPR161 multiplier:", 2);
     Dialog.show();
     ciliaChannel = Dialog.getChoice();
     GPR161Channel = Dialog.getChoice();
     ciliaMultiplier = Dialog.getNumber();
     gpr161Multiplier = Dialog.getNumber();

//2. Initialize key variables, ROI manager and results table  

lower = 0
run("Set Measurements...", "area mean limit decimal=5");
run("Clear Results");

name=getTitle();
bit_depth = bitDepth();

if(bit_depth == 16){
	maxThresh = 65535;
}

if(bit_depth == 32){
	maxThresh = 4000000000;
}

Stack.setChannel(GPR161Channel);
run("Select None");

selectWindow(name);
Stack.setChannel(ciliaChannel);
run("Select None");

setBatchMode(false);

if(isOpen("GPR161_Quant")){
	selectWindow("GPR161_Quant");
	run("Close");
}

if(isOpen("Results")) {
	selectWindow("Results");
	run("Close");
}

roiManager("reset");

run("Clear Results");
		result = nResults();
            setResult("Image", result, 1); 
            setResult("Cilium", result, 1); 
            setResult("Cilia px", result, 1);
            setResult("161 px", result, 1); 
            setResult("% 161", result, 1);
            setResult("mean 161", result, 1);
            setResult("mean Arl13", result, 1);
            setResult("161/Arl13", result, 1);
            updateResults();

            IJ.renameResults("GPR161_Quant");

//3. Create a loop to iterate the macro for each open window 

Images = nImages();
for ( a=0; a<Images; a++){
	
	name=getTitle();
	run("Gaussian Blur...", "sigma=0.50 stack");
	Stack.setChannel(ciliaChannel);
	run("Red");
	Stack.setChannel(GPR161Channel);
	run("Green");
	run("Set Scale...", "distance=0 known=0 pixel=1 unit=pixel");

//4. automatically locate cilia in the image


	run("Duplicate...", "duplicate channels=["+ciliaChannel+"]-["+ciliaChannel+"]");
	run("Measure");
	overallMean = getResult("Mean", 0);
	overallLower = overallMean * ciliaMultiplier;
	setThreshold(overallLower, maxThresh);
	run("Properties...", "channels=1 slices=1 frames=1 unit=um pixel_width=1 pixel_height=1 voxel_depth=1");
	run("Set Measurements...", "centroid redirect=None decimal=5");
	run("Clear Results");
	run("Analyze Particles...", "size=100-5000 circularity=0.00-0.80 display");
	close();
	
//5. Mark the cilia and save an RGB image for auditing

	
	numPoints = nResults();
	for (q=0; q<numPoints; q++){
    	setTool("multipoint");
    	x1 = getResult("X", q);
    	y1 = getResult("Y", q);
    	setKeyDown("Shift");
    	makePoint(x1, y1);		
    	setKeyDown("none");
	}

	getSelectionCoordinates(x, y);
	
	selectWindow(name);
	run("Duplicate...", name+"_cilia");
	run("RGB Color");
	run("Restore Selection");
	saveAs("Tiff", dir+name+"_cilia");
	close();

//5a. record the cilia positions in the ROI manager

	selectWindow(name);
	n = x.length;
	for (r=0; r<n; r++){
    	makePoint(x[r], y[r]);
    	roiManager("Add");
		}

//6.  Measure each detected cilium

	selectWindow("Results"); 
    run("Clear Results");
    run("Close");		
	run("Set Measurements...", "area mean limit decimal=5");

n = roiManager("count");
for ( i=0; i<n; i++){
    nCilium=i+1;
    selectWindow(name);
    roiManager("Select", i);
    makeRectangle(x[i]-35, y[i]-35, 70, 70);
    Stack.setChannel(ciliaChannel);
    run("Duplicate...", name);
    cilium = getTitle();

    selectWindow(name);
    Stack.setChannel(GPR161Channel);
    run("Duplicate...", name);
    GPR161 = getTitle();
    
    selectWindow(cilium);
    run("Measure");
	mean1 = getResult("Mean", 0);
	ciliaLower = mean1 * ciliaMultiplier;
    setThreshold(ciliaLower, maxThresh);
    run("Measure");
    ciliaArea = getResult("Area", 1);

//6a. check whether a cilium exists at the marked position
    
    if(ciliaArea>0) {
    
    meanCilium = getResult("Mean", 1);
    run("Create Selection");
    run("Make Inverse");

//7.  Delete GPR161 data outside the cilium and measure ciliary GPR161
	
    selectWindow(GPR161);
    run("Measure");
	mean2 = getResult("Mean", 2);
	lower = mean2 * gpr161Multiplier;
    run("Restore Selection");
    run("Make Inverse");
    run("Measure");
    mean161 = getResult("Mean", 3);
    run("Make Inverse");
    run("Clear");
    run("Select None");
	setThreshold(lower, maxThresh); 
	run("Measure");
    GPR161Area = getResult("Area", 4);

    selectWindow("Results"); 
    run("Clear Results");
    run("Close");

//8.  Add results to the results table

            selectWindow("GPR161_Quant");
            IJ.renameResults("Results");
            result = nResults();
            setResult("Image", result, name); 
            setResult("Cilium", result, nCilium); 
            setResult("Cilia px", result, ciliaArea);
            setResult("161 px", result, GPR161Area); 
            setResult("% 161", result, 100*GPR161Area/ciliaArea);
            setResult("mean 161", result, mean161); 
            setResult("161/Arl13", result, mean161/meanCilium);
            setResult("Lower 161 Threshold", result, lower);
            setResult("Lower Cilia Threshold", result, ciliaLower);
            updateResults();
            IJ.renameResults("GPR161_Quant");
            selectWindow("GPR161_Quant");
			saveAs("Results", dir+"GPR161_Quant");

//9. create and save a RGB image of each cilium for making grids later               
    
    selectWindow(cilium);
    run("Make Inverse");
    rename("red");
	getMinAndMax(min, max);
	setMinAndMax(min, max*0.75);
	run("Select None");
	
    selectWindow(GPR161);
    resetThreshold();
    run("Select None");
    run("Subtract...", "value=["+mean2+"]");
    rename("green");
	gpr161Level = mean2 * 6;
	getMinAndMax(min, max);
	setMinAndMax(min, gpr161Level);

	run("Merge Channels...", "red=red green=green blue=*None* gray=*None* create");
    selectWindow("Composite");
    run("RGB Color");
    rename("grid_" + name + "_" + nCilium);
    merge=getTitle();
	saveAs("Tiff", dir+merge);
	run("Close");
	run("Clear Results");
	selectWindow("Composite");
	run("Close");
	}

//6b. if the marked position has no cilia, close the sub-windows.

	else {
		selectWindow(cilium);
		close();
		selectWindow(GPR161);
		close();
	}

//6c. Close the cilia-measurement loop
	
	}                

	selectWindow(name);
	close();

	roiManager("reset");

	}

print("Finished");

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html