// Optimized hiPSC Co Culture in vitro segmentation - Bardy Lab - SAHMRI

// For Bardy lab, load the first raw image of your timeseries stack with this command:
// ================================================
run("Image Sequence...");

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

// For Wilson, Load images with this command:
// ================================================
// run("TIFF Virtual Stack...");
// =================================================


orig = getTitle();
filename = File.nameWithoutExtension();	
dir = getDirectory("image")
filelist = getFileList(dir);
Array.show(filelist);

zpd = dir+"\\zproj";
File.makeDirectory(zpd);
crp = dir+"\\crop"; 
File.makeDirectory(crp);
resdir = dir+"\\results"; 
File.makeDirectory(resdir);

// Get rid of heading
//makeRectangle(0, 8, 2048, 2040);
//run("Duplicate...", "duplicate");


// Denoise
imgc = getTitle();

run("Enhance Contrast", "saturated=0.5");
//run("Despeckle", "stack");



// Z-projections to get ranges
run("Z Project...", "start=5 stop=60 projection=[Average Intensity]");
makeRectangle(0, 9, 2048, 2040);
run("Crop");
avgip = getTitle();	// ** SAVE AVGIP **
Stack.getDimensions(width, height, channels, slices, frames);
run("Subtract Background...", "rolling=50");
//run("Gaussian Blur...", "sigma=2");
run("Enhance Contrast", "saturated=1");
run("Green Fire Blue ");
run("Enhance Local Contrast (CLAHE)", "blocksize=85 histogram=256 maximum=20.00 mask=*None* fast_(less_accurate)");


// Threshold the difference range and convert to ROIs
selectWindow(avgip)
run("Morphological Filters", "operation=Opening element=Disk radius=5");
rename("filt");
filt = getTitle();

// Segment out optical fiber (if you need to... otherwise comment out). 
waitForUser("Use the Rectangle ROI to cut out the fiber from the FOV")

run("Find Maxima...", "prominence=15 output=[Point Selection]");
roiManager("add");

// Make Marker image for watershed
newImage("markers", "8-bit black", width, height, 1);
markers = getTitle();
roiManager("Select", 0);
roiManager("Fill");
run("Make Binary");

// Get Mask for watershed Input
selectWindow(filt)
run("Duplicate...", " ");
setAutoThreshold("Huang dark")
run("Threshold...")
run("Convert to Mask");
rename("mask");
mask = getTitle()  // ** SAVE MASK **

// Get gradient image for watershed input
run("Morphological Filters", "operation=Gradient element=Disk radius=2");
rename("grad");
grad = getTitle();

// Run watershed Segmentation
run("Marker-controlled Watershed", "input=grad marker=markers mask=mask binary calculate use");
wtsh = getTitle(); 
run("Fire");
roiManager("reset")

run("Duplicate...", " ");
setThreshold(0.1000, 1000000000000000000000000000000.0000);
run("Convert to Mask");
wsmsk = getTitle();

// Circle cell data
run("Analyze Particles...", "size=250-Infinity pixel circularity=0.10-1.00 add");
close(wsmsk);

// Apply ROIs to Original image
selectWindow(orig)

roiManager("Show None");
roiManager("Show All");
roiManager("multi-measure one")	// **SAVE RESULTS **

resname = "results_"+orig;
print(resname)
restitle = filename;
saveAs("Results", resdir+'\\'+restitle+"_results.csv");

roiManager("Save", resdir+'\\'+restitle+"_roi.zip");
close("Results");
roiManager("reset")

selectWindow(avgip)
rename("avgip_"+filename);
roiManager("Show None");
roiManager("Show All");
saveAs("tiff", zpd+"\\"+getTitle());
avgip = getTitle();


// CLOSE stuff
close(markers);
close(grad);
close(wtsh);
close(mask);
close(filt)
//close(orig);
//close(avgip);





// Interactive data viewing


// Make a pretty video
// AvgIP --> Create Background --> Subtract from denoised stack
//