// Ratiometric Fluroescence Processing

// This script takes a hyperstack of timeseries images containing multiple different simultaneously applied image channels, and spits out ratiometric fluorescence and indivisual fluorescence data. 

// One fluorescence channel will be used to segment the ROIs to extract fluroescence timeseries across multiple image channels. Background subtracted fluorescence channels will be ratioed to provide a 32-bit depth ratiometric fluorescence stack based on laurdan/di-4-ANNEPS general polarization. This is calculated by ratioing the difference and sum images of the two fluroescence channels in question. 

// Since the raw fluroescence timeseries extraction does not take into account backround subtraction done computationally, an additional background ROI will need to be obtained for such calculations in each channel. 

// For di-4 experiments, ChanA(1) and ChanD(3) need to be compared for GP. ChanC (2) is the fSHG signal that kindof worked? not sure if it'll be useful for anything though. low signal. difficult to see. 

// Written by Wilson Adams | Vanderbilt Biophotonics Center | 21 Oct 2020

// =====================================================================================
// Run this section for single experiment
hypstk = getTitle();
di4ratioprocess(hypstk);


// =============================================
// Run this section for batch processing
/*
	dir = getDirectory("Navigate to the folder to be processed");
	print("STARTING file extraction: "+dir);
	
	filelist = getFileList(dir);
	Array.show(filelist);
	
	// Loop through each folder in the parent directory
	
	//for (i = 0; i < 1; i++) {			// for debugging one folder at a time
	for (i = 0; i < filelist.length; i++) {
	    if (endsWith(filelist[i], ".tif")) {  
	    	open(filelist[i]);
	    	hypstk = getTitle();
			di4ratioprocess(hypstk);
			print(hypstk+" COMPLETE"); 
	    }
	}
	print("Processing of " + dir + " FINISHED");
*/
// ==============================================
// ==============================================

function di4ratioprocess(hypstk) { 
	// Base function for ratiometric fluorescence processing. Feed into batch processing
	roiManager("reset");

	run("Duplicate...", "duplicate channels=1");
	fluorA = getTitle();
	
	selectWindow(hypstk);
	run("Duplicate...", "duplicate channels=3");
	fluorB = getTitle();

	// AvgIP to get image that will become mask
	selectWindow(hypstk);
	run("Z Project...", "projection=[Average Intensity]");
	premask = getTitle();
	run("Duplicate...", " ");
	setAutoThreshold("Li dark");
	setOption("BlackBackground", false);
	run("Convert to Mask");
	mask = getTitle();

	// bulk selection
	run("Create Selection");
	roiManager("Add");

// CELL SELECTION - Seeded watershed of inverted blurred avgip

	// Segmentation exclusion area from cell border mask
	selectWindow(premask);
	run("Duplicate...", " ");
	cellexc = getTitle();
	roiManager("Show None");
	run("Gaussian Blur...", "sigma=2");		// to optimize maxima/min finding
	setAutoThreshold("Li dark");
	//run("Threshold...");
	run("Convert to Mask");
	run("Close-");
	run("Fill Holes");		// to make a solid mask

	// Generate seed points for watershed from inverted avgIP
	selectWindow(premask);
	run("Duplicate...", " ");
	map = getTitle();
	run("Gaussian Blur...", "sigma=1");
	run("Find Maxima...", "prominence=9 light output=[Point Selection]");
	roiManager("Add");
	getDimensions(width, height, channels, slices, frames);
	newImage("ptseed", "8-bit", width, height, 1);
	run("Invert");
	roiManager("select", 1); // select the points selection
	roiManager("fill");
	roiManager("delete");	// get rid of point seed ROI, not needed anymore
	ptseed = getTitle();
	//run("Duplicate...", " ");
	//distmap = getTitle(); 
	//run("Distance Map");
	
	
	// Run Watershed
	run("Marker-controlled Watershed", "input=&cellexc marker=&ptseed mask=&cellexc binary calculate use");

	wtshd = getTitle(); 
	setThreshold(0.1500, 1000000000000000000000000000000.0000);
	run("Convert to Mask");
	run("Analyze Particles...", "add");

	// Get edge ROIs - wtshd AND mask pixels within wtshd ROIs
	//imageCalculator("AND create", "input binary full mask","AVG_Merged-1");
	//wtshdedge = getTitle(); 

	// Save ROIs. 
	selectWindow(hypstk);
	dir = getInfo("image.directory");	
	resdir = dir+"/Results/";
	hypstk1 = File.getNameWithoutExtension(hypstk);
	
	roiManager("show none");
	roiManager("save", resdir+hypstk1+"_roi.zip");
	
	// Run on hyperstack
	selectWindow(hypstk);
	roiManager("Multi Measure");
	saveAs("Results", resdir+hypstk+"_results.csv");

	// Get GP Image Stack
	selectWindow(fluorA);
	imageCalculator("Subtract create 32-bit stack", fluorA, fluorB);
	diff = getTitle();
	imageCalculator("Add create 32-bit stack", fluorB, fluorA);
	sum = getTitle();
	imageCalculator("Divide create 32-bit stack", diff, sum);
	gp = getTitle(); 

	// Get GP ROI data
	selectWindow(gp);
	saveAs("Tiff", resdir+hypstk1+"_GP.tiff");
	roiManager("Multi Measure");
	saveAs("Results", resdir+hypstk+"_GPresults.csv");
	
	// Close everything and move on
	close("*");
	close("Results");
	run("Close All");
	roiManager("reset");
	
}
	// Should i be analyzing on a per-cell basis? or should i stick with bulk experiment analysis?