// AstroSeg.ijm

// Trying to build a better astrocyte segmentation algorithm for Ca Imaging analysis. Doing it by hand sucks. 
// Ideally, it'll be fully automated. That's hard... so some user input may be necessary

// Wilson R Adams | Vanderbilt Biophotonics Center | 1 July 2019

// Thresholding alone doesn't work well. 
// I want to try Watershed segmentation. 

stk = getTitle()
fn = File.nameWithoutExtension()
File.getName(path)
par = File.getParent(path)

// full stack stdev appears to have some risidual motion artifacts.
// stdev project over stim location appears to give a better stk to work with. 
// stdev project plus RBF bkgd subtraction makes a nice place to start. Works Well
// Trying bkgd subtract based of inital frames, then SumIP. 
stk = getTitle()

// get abg intensity projection
selectWindow(stk)
run("Z Project...", "stop=50 projection=[Average Intensity]");
bkgd_avg = getTitle()

run("Calculator Plus", "i1=&stk i2=&bkgd_avg operation=[Subtract: i2 = (i1-i2) x k1 + k2] k1=1 k2=0 create");
rename(stk+"_bkgdCorr1")
bkgd_corr1 = getTitle()
close(bkgd_avg)
run("Z Project...", "projection=[Standard Deviation]");
run("Enhance Contrast", "saturated=0.35");
stdIP = getTitle()
run("!Blueish")
run("Subtract Background...", "rolling=60 sliding disable");
run("Enhance Contrast", "saturated=0.35");

selectWindow(bkgd_corr1)
run("Z Project...", "projection=[Max Intensity]");
run("!Blueish")
maxIP = getTitle()
run("Subtract Background...", "rolling=60 sliding disable");
run("Enhance Contrast", "saturated=0.35");

// only really works for responding cells. Need something that works for nonresponding cells too. 

// Using Find Maxima on blurred MaxIP is getting very close results
// ... not sure how repeatable. but hopeful. 
// still highlights endfeet as local max... which makes sense 
// maybe a Frangi filter would help?

// I tried the Morphological Filters plugin to enhance nuclei, worked pretty well... closest I've gotten so far. +8px Gaussian Blur. find max. 

run("Morphological Filters", "operation=Opening element=Disk radius=6");
run("Gaussian Blur...", "sigma=8");
run("Find Maxima...", "prominence=40 strict output=[Segmented Particles]");

// For nonresponding stacks:
// avgIP/sumIP, big Gauss blur(50px), subtract. --> convert to 32bit
stk = getTitle()
	//filepath conventions
path = File.getName(stk)
par = File.directory
print(par); print(path)
fn = File.nameWithoutExtension()
savedir = par+"bkgd\\"
File.makeDirectory(savedir)		// Make processed images directory in parent directory
print(savedir)

	//bkgd sub
stk = getTitle()
fn = 
run("Z Project...", "stop=50 projection=[Sum Slices]");
sumIP = getTitle()
run("Duplicate...", " ");
sumBkgd = getTitle()
run("Gaussian Blur...", "sigma=50");
run("Calculator Plus", "i1=&sumIP i2=&sumBkgd operation=[Subtract: i2 = (i1-i2) x k1 + k2] k1=1 k2=0 create");

selectWindow(stk)
run("Z Project...", "stop=50 projection=[Average Intensity]");
avgBkgd = getTitle()
run("Gaussian Blur...", "sigma=50");
getStatistics(area, mean, min, max, std, histogram)
run("Calculator Plus", "i1=&stk i2=&avgBkgd operation=[Subtract: i2 = (i1-i2) x k1 + k2] k1=1 k2=&mean create");
rename("_bkgdSub)
stk_bkgd = getTitle()
close(sumIP); close(sumBkgd);
close(avgBkgd)
run("!Blueish")

selectWindow(stk_bkgd)
rename(fn+"_bkgdcorr")
stk_bkgd = getTitle()
save(savedir+stk_bkgd)


// LEVEL SET SEGMENTATION MIGHT WORK!!!! HOLY CRAP!! THANKS JEREMY!
run("Level Sets", "method=[Active Contours] use_fast_marching use_level_sets grey_value_threshold=30 distance_threshold=1 advection=2.20 propagation=1 curvature=1 grayscale=5 convergence=0.0050 region=outside");
//better. needed to apply a fixed LUT dynamic range to make this work.
run("Enhance Contrast", "saturated=20");
run("Level Sets", "method=[Active Contours] use_fast_marching use_level_sets grey_value_threshold=10 distance_threshold=1 advection=2.20 propagation=1 curvature=1 grayscale=5 convergence=0.0050 region=outside");

// Find Maxima - helpful to run on smoothed data. sensitive to noise
// Trying out different smoothing parameters. 
prom = 4000;
win = getTitle()

sig = 2;
selectWindow(win)
run("Duplicate...", getTitle()+"_dup2")
run("Gaussian Blur...", "sigma=&sig");
run("Find Maxima...", "prominence=&prom strict output=[Point Selection]");
roiManager("Add");

sig = 4
selectWindow(win)
run("Duplicate...", getTitle()+"_dup4")
run("Gaussian Blur...", "sigma=&sig");
run("Find Maxima...", "prominence=&prom strict output=[Point Selection]");
roiManager("Add");

sig = 8
selectWindow(win)
run("Duplicate...", getTitle()+"_dup8")
run("Gaussian Blur...", "sigma=&sig");
run("Find Maxima...", "prominence=&prom strict output=[Point Selection]");
roiManager("Add");

sig = 16
selectWindow(win)
run("Duplicate...", getTitle()+"_dup16")
run("Gaussian Blur...", "sigma=&sig");
run("Find Maxima...", "prominence=&prom strict output=[Point Selection]");
roiManager("Add");

sig = 32
selectWindow(win)
run("Duplicate...", getTitle()+"_dup32")
run("Gaussian Blur...", "sigma=&sig");
run("Find Maxima...", "prominence=&prom strict output=[Point Selection]");
roiManager("Add");