// DAPI Nuclear Segmentation and Counting

/*
A general purpose script for segmenting nuclei from fluorescence images
Ideally DAPI or Ethidium Bromide, but potentially other nuclear shaped
objects in the right context. '

Be sure to put in the rough DIAMETER of the nuclei in your images (in pixels).
You can use the 'drawCircle' tool to estimate roughl how big the average nucleus
is in your image. 
*/

// ======================================================
// Written by Wilson R Adams | Siciliano Lab | 17 Jan 2022
// ======================================================

// USER INPUT 
nuc_diameter = 20; 

// ======================================================
// Housekeeping stuff up front:
nuc_radius = nuc_diameter/2; 
roiManager("reset");

// Using the open image
raw = getTitle(); 

Stack.getDimensions(width, height, channels, slices, frames);
//Array.show("Dims", Array.concat(channels, slices, frames));

// Max Intensity Projection if the raw data is a Z stack
if (slices>1) {
	run("Z Project...", "projection=[Max Intensity]");
}

start = getTitle(); 
Stack.getPosition(chan, sl, fr);


// Isolate Nuclei Iamge Channel (assuming it's the current channel viewed)
run("Duplicate...", " ");
nuclei = getTitle(); 

// Despeckle and smooth
run("Median...", "radius=3");
run("Subtract Background...", "rolling=&nuc_radius");

kernel = 0.2*nuc_diameter;
run("Gaussian Blur...", "sigma=&kernel");
run("Duplicate...", " ");
filt = getTitle(); 

// Thresholding of bright image features
resetMinAndMax();
run("Enhance Contrast", "saturated=0.35");
setAutoThreshold("Huang dark no-reset stack");
setOption("BlackBackground", false);
run("Convert to Mask");

// Split up binary mask
run("Distance Transform Watershed", "distances=[Borgefors (3,4)] output=[16 bits] dynamic=1 connectivity=8");
mask = getTitle(); 

setThreshold(1, 65535);
run("Convert to Mask");
bin = getTitle();

// Get ROIs from mask
run("Analyze Particles...", "size=10-Infinity clear add");

// Close out un-needed images
/*
close(filt);
close(mask); 
close(bin); 
*/

// Measure the size and brightness of each nucleus
run("Set Measurements...", "area mean redirect=None decimal=3");
selectWindow(nuclei);
roiManager("multi-measure");
roiManager("count");