// Astro Multipulse Automated PRocessing
//================================================================
// This macro is the first preprocessing step in a series of macros to somewhat automagically extract calcium time series data from calcium imaging experiments performed on the STORM scope. 

// This macro will open up each experiment, and calculate Avg and Std Z projections from each experiment, then repeats this action for each experiment in the directory (assuming it is a stack, and not a static image). 

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

// Select folder directory to process
dir = getDirectory("Navigate to the folder to be processed");
filelist = getFileList(dir);
Array.show(filelist);

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

// Generate file list

// Loop through file list and process each image
for (i = 0; i < lengthOf(filelist); i++) {
    if (endsWith(filelist[i], ".nd2")) { 
    	
    	// Generate file name string
        file = dir + File.separator + filelist[i];
        
        // Open as a VIRTUAL STACK (Faster, memory saving)
        run("Bio-Formats Importer", "open=&file color_mode=Default rois_import=[ROI manager] split_channels view=Hyperstack stack_order=XYCZT use_virtual_stack");
        orig = getTitle();
        run("Out [-]");
		filename = File.nameWithoutExtension();		
		print(filename);
        // Go to frame 65
        if (nSlices == 1) {
        	break;
        }
	    setSlice(65);
		run("Enhance Contrast", "saturated=0.35");
		
	    // Get Avg and StdIP
	    run("Z Project...", "start=55 stop=70 projection=[Average Intensity]");
	    run("Out [-]");
	    rename("avgip_"+filename);
	    saveAs("tiff", zpd+"\\"+getTitle());
	    avgip = getTitle();

		selectWindow(orig);
	    run("Z Project...", "start=55 stop=70 projection=[Standard Deviation]");
	    run("Out [-]");
	    rename("stdip_"+filename);
	    saveAs("tiff", zpd+"\\"+getTitle());
	    stdip = getTitle();
	    selectWindow(stdip);
	    close(stdip);     

		AstroSeg_virtual(orig, avgip);
		
		selectWindow(orig); 
		restitle = File.nameWithoutExtension();
		saveAs("Results", resdir+'\\'+restitle+"_results.csv");
		close(orig);
		close(avgip);
		roiManager("Save", resdir+'\\'+restitle+"_roi.zip");
		close("Results");
		roiManager("reset")
    } 
}


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

// Gaussian Blur (7px kernel)
run("Gaussian Blur...", "sigma=7");

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

}