// MANTIS Multimodal Cell Stim/Damage processing code

// The goal here is to 
// 1) segment out each cell from the SRS images, 
// 2) obtain SRS and 2P viability (_PMT1) time series, 
// 3) and cell centroid locations
// 4) save to a file to send to matlab. 
// 5) Repeat for the rest of the experiments
// ===================================================

//srspath = File.openDialog("Select Segmentable Image Stack");
//viapath = File.openDialog("Select Functional Image Stack");
//directory = File.getParent(srspath);

//open(srspath);
//srs = getTitle();
//open(viapath);
//via = getTitle();

function multimode_seg_extract(seg, func, dir) { 
// This function segments the stack 'seg' and applies the ROIs to 
// Functional image stack 'func' and 'seg' while saving results to 
// directory 'dir'

// Segmentation done with modified AstroSeg_Virtual routine

// All data saved to be imported and analyzed in program of choice

// Written by Wilson Adams, Vanderbilt Biophotonics Center, July 2020

srs = seg; 	// Segmentable Image Stack
via = func;	// Functional Image Stack (Viability, Calcium, etc)

// Get ROIs from AvgIP of SRS stack
multimode_seg_extracter(srs);
srsavgip = getTitle();

// Get Centroids of ROIs for location transform
run("Set Measurements...", "centroid redirect=None decimal=6");
roiManager("multi measure one");
saveAs("Results", dir+srs+"_centroid");
roiManager("save", dir+srs+"_roi.zip");

// Get timeseries
run("Set Measurements...", "mean redirect=None decimal=6");

vianame = File.getNameWithoutExtension(dir+via);
selectWindow(via);
roiManager("show all");
roiManager("multi-measure one");
saveAs("Results", dir+vianame+"_timeseries");
close("Results");

srsname = File.getNameWithoutExtension(dir+srs);
selectWindow(srsname);
roiManager("show all");
roiManager("multi-measure one");
saveAs("Results", dir+srsname+"_timeseries");
close("Results");

selectWindow(srsavgip);
saveAs("tif", dir+srsname+"_avgip"+".tif");

roiManager("reset");
}

// ============================================

function multimode_seg_extracter(segstack) { 
	// given a segmentable (segstack) and functional (funcstack) 
	// image stacks,this macro will segment cells from image 

	// step 1 - get segmentable image (Avg IP from seg stack)
	selectWindow(segstack);
	run("Z Project...", "projection=[Average Intensity]");
	segAvgIP = getTitle();

	// Step 2 - run segmenting code on it to get ROIs. 
	AstroSeg_virtual(segAvgIP); // Modify function below
}




// ===========================================

// Run on AvgIP to debug
//test = getTitle();
//AstroSeg_virtual(test);

function AstroSeg_virtual(avgip) { 
// Modified from Astrocyte Segmentation

// 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=2");
run("Median...", "radius=1");
run("Subtract Background...", "rolling=30");
run("Enhance Contrast", "saturated=5.00");




// Additional Flatfield Correction if needed for cell activation
run("Enhance Local Contrast (CLAHE)", "blocksize=20 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=10-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");

close(wkimg);
close(mask);
close(grad);
close(markers);
close(wtsh);
close(wsmsk);
}