// 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. 

// run("Z Project...", "start=5 stop=40 projection=[Average Intensity]");
//orig = getTitle();
roiManager("reset");
avgip = getTitle();
AstroSeg_virtual(avgip, avgip) 


function AstroSeg_virtual(orig, avgip) { 
// function description

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

// Sub Background (50px kernel)
run("Duplicate...", " ");
wkimg = getTitle();
// Gaussian Blur (7px kernel)
run("Gaussian Blur...", "sigma=4");
run("Median...", "radius=8");
run("Subtract Background...", "rolling=100");
run("Enhance Contrast", "saturated=5.00");




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

// =============== USER MODIFY ===============
// Threshold image to generate Mask (Otsu(astro) - Li (BV2)). Erode a bit. Fill holes.
setAutoThreshold("Huang 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();

// =============== USER MODIFY ===============
// Find Local Maxima (Prominence 20 - Astro | 55 BV2), add to ROI manager
selectWindow(wkimg);
run("Find Maxima...", "prominence=10 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
run("Analyze Particles...", "size=300-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");

// Multimeasure on the orig stack. 
run("Set Measurements...", "mean redirect=None decimal=6")
selectWindow(orig);
roiManager("multi-measure one");

// Save output into a new subdirectory
// This will be done in the parent function

// Close everything
close(wkimg);
close(mask);
close(grad);
close(markers);
close(wtsh);
close(wsmsk);

}