// AstroSegVirtual
//================================================================
// This macro does the bulk of the important stuff in the calcium image processing pipeline. This script will attempt to identify any number of cells in the given AvgIP image based on a marker-controlled watershed segmentation approach. This allows for optimal segmentation of cell bodies in images where cells may be overlapping, dendritic, really packed together, or any number of reason that would prohibit basic thresholding from working. 

// If it does screw up, you have the option to delete and add any ROIs before the data is extracted and saved for record. 

// Written by Wilson R Adams | Vanderbilt Biophotonics Center | Nov 2019
//================================================================

// Begin by running this macro and selecting the directory of experiments you with to process. The script should do the rest. 
roiManager("reset");

// Open AvgIP file, get image stats
avgip = getTitle();
Stack.getDimensions(width, height, channels, slices, frames);
getPixelSize(unit, pixelWidth, pixelHeight);

// Sub Background (50px kernel)
run("Duplicate...", " ");
wkimg = getTitle();
	// SET BKGD SUBTRACTION KERNEL SIZE (~2-3x feature size)
run("Subtract Background...", "rolling=50");
run("Enhance Contrast", "saturated=5.00");

// Gaussian Blur (7px kernel) 
// ==== ADJUST TO YOUR NEEDS. JUST ENOUGH TO SMOOTH OUT NOISE===
run("Gaussian Blur...", "sigma=3");

// Additional Flatfield Correction if needed for cell activation (accomodating for brightness differences between cells. 
run("Enhance Local Contrast (CLAHE)", "blocksize=199 histogram=256 maximum=10.00 mask=*None* fast_(less_accurate)");
run("Duplicate...", " ");

// Threshold image to generate Mask (Otsu(astro) - Li (BV2)). Erode a bit. Fill holes.
setAutoThreshold("Li dark");
//run("Threshold...");
run("Convert to Mask");
rename("mask")
mask = getTitle();

// Make Gradient mask image
run("Morphological Filters", "operation=Gradient element=Disk radius=2");
rename("grad")
grad = getTitle();

// Find Local Maxima (Prominence 20 - Astro | 55 BV2), add to ROI manager
selectWindow(wkimg);
run("Find Maxima...", "prominence=55 output=[Point Selection]");
roiManager("add");

// Make blank marker image. Apply points to img with 'Fill'
// You cant use ROImanager point selections. you have to have an image with points identifies as maximum seed points.
newImage("markers", "8-bit black", width, height, 1);
markers = getTitle();
roiManager("Select", 0);
roiManager("Fill");
run("Make Binary");
run("Properties...", "channels=1 slices=1 frames=1 unit=&unit pixel_width=&pixelWidth pixel_height=&pixelHeight voxel_depth=1.0000");

// Marker based watershed (input = gradient, binary markers, mask)
run("Marker-controlled Watershed", "input=grad marker=markers mask=mask binary calculate use");
wtsh = getTitle(); 
run("Fire");
roiManager("select", 0);
roiManager("delete");

// Threshold watershed map 1:inf for Analyze Particles
run("Duplicate...", " ");
setThreshold(0.1000, 1000000000000000000000000000000.0000);
run("Convert to Mask");
wsmsk = getTitle(); 

// Analyze particles --> min area = 200 into ROI. Save
// === SET SIZE AND CIRCULARITY TO SUIT YOUR NEEDS =======
run("Analyze Particles...", "size=50-Infinity pixel circularity=0.10-1.00 add");

// Check and Update ROIs (if needed, probably needed...). Save to new subdir

selectWindow(avgip); 
roiManager("show all with labels");

// Put all your measurements you want here. e.g. I use mean intensity and centroid position here.  
run("Set Measurements...", "mean centroid redirect=None decimal=6");
selectWindow(avgip);
roiManager("multi-measure one");

// Close everything (or leave them open if you want. Up to you). 
close(wkimg);
close(mask);
close(grad);
close(markers);
//close(wtsh);
//close(wsmsk);