//SAMA
version=0.99;
//Copyright (C) 2015 Maël Montévil and Tessie Paulose 

//    This program is free software: you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation, either version 3 of the License, or
//    (at your option) any later version.

//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.

//    You should have received a copy of the GNU General Public License
//    along with this program.  If not, see <http://www.gnu.org/licenses/>.




print("***************************");
print("Starting SAMA images  version "+version);

//creates a directory if it does not exist yet
function createdir(dir, suffix) {
  dirres=dir+suffix+"/";
  if(!File.isDirectory(dirres)){
    File.makeDirectory(dirres);
  }
  return dirres;
}


//function to bypass 3D suite weird renaming
function measuresave(dir, name) {
Ext.Manager3D_SaveResult("M",dir  + name);
    File.rename(dir+"M_"+name , dir+name);
    return true;
}

//enhances the edges of structure in order for them not to depend much on the threshold
function enhancedges(){
  name=getTitle();
  run("Duplicate...", "title=Duplicate duplicate range slices=1-"+nSlices);
  run("Variance 3D...", "x=2 y=2 z=2");
  run("Gamma...", "value=0.50 stack");
  maxx=0;
  minn=10000;
  for (j=1; j<nSlices+1; j++) {
    setSlice(j);
    getRawStatistics(nPixels, mean, min, max, std, histogram);
    minn=minOf(minn,min);
    maxx=maxOf(maxx,max);
  }
  setMinAndMax(minn, maxx);
  run("8-bit");
  run("Multiply...", "value=0.7 stack");
  imageCalculator("Add create 32-bit stack", name,"Duplicate");
  close("Duplicate");
  
  maxx=0;
  minn=10000;
  for (j=1; j<nSlices+1; j++) {
    setSlice(j);
    getRawStatistics(nPixels, mean, min, max, std, histogram);
    minn=minOf(minn,min);
    maxx=maxOf(maxx,max);
  }
  setMinAndMax(minn, maxx);
  run("8-bit");
  close(name);
  rename(name);
}


//function to correct the loss of luminosity in the depper part of confocal stacks
function attenuation(){
  run("Enhance Contrast...", "saturated=0.001 normalize process_all use");
  run("Median 3D...", "x=1 y=1 z=1");
  oldm=1;
  for (j=2; j<nSlices; j++) {
    maxx=0;
    for (k=-1; k<2; k++) {
      setSlice(j+k);
      getRawStatistics(nPixels, mean, min, max, std, histogram);
      maxx=maxOf(maxx,max);
    }
    newm=255/maxx;
    run("Multiply...", "value="+newm+" ");
  }
  run("Median 3D...", "x=1 y=1 z=1");
  run("Gaussian Blur 3D...", "x=0 y=0 z=1");
  run("Minimum 3D...", "x=1 y=1 z=1");
}

// performs the preprocessing for the basic morphometrics	
function preprocessBasic(mode){
  name=getTitle();
  if(mode=="Enhance boundaries"){enhancedges();}
  if(mode=="Split neighbors"){
    run("Gradient (3D)");
    imageCalculator("Subtract create stack", name,"Rebinned");
    close("Rebinned");
  }
  run("Subtract Background...", "rolling=50 sliding stack");
  attenuation();
  attenuation();
}

//main function to peform basic morphometrics
function mainstructures(dirsource, dirmap,dirdata,thres,chull,mode) {
  list = getFileList(dirsource);
  setBatchMode(true);
  if(chull){
    run("3D Manager Options", "volume surface compactness fit_ellipse 3d_moments convex_hull integrated_density mean_grey_value std_dev_grey_value feret minimum_grey_value maximum_grey_value centroid_(pix) centroid_(unit) distance_to_surface centre_of_mass_(pix) centre_of_mass_(unit) bounding_box radial_distance surface_contact distance_between_centers=10 distance_max_contact=1.80");
  }else{
    run("3D Manager Options", "volume surface compactness fit_ellipse 3d_moments integrated_density mean_grey_value std_dev_grey_value feret minimum_grey_value maximum_grey_value centroid_(pix) centroid_(unit) distance_to_surface centre_of_mass_(pix) centre_of_mass_(unit) bounding_box radial_distance surface_contact distance_between_centers=10 distance_max_contact=1.80");
  }  
  run("3D OC Options", "volume surface nb_of_obj._voxels nb_of_surf._voxels integrated_density std_dev_gray_value centroid mean_distance_to_surface std_dev_distance_to_surface median_distance_to_surface centre_of_mass bounding_box dots_size=5 font_size=10 store_results_within_a_table_named_after_the_image_(macro_friendly) redirect_to=none");
  for (i=0; i<list.length; i++) {
    showProgress(i+1, list.length);
    open(dirsource+list[i]);
    name=getTitle();
    preprocessBasic(mode);
    //run("Gradient (3D)");
    //imageCalculator("Subtract create stack", name,"Rebinned");
    //close("Rebinned");
    //saveAs("TIFF", dir2+"Result of_"+list[i]);
    setThreshold(thres, 255);
    setOption("BlackBackground", true);
    run("Convert to Mask", "method=Default background=Dark black");
    run("Options...", "iterations=1 count=1 black edm=Overwrite do=Nothing");
    run("Fill Holes", "stack");
    run("3D Objects Counter", "threshold="+thres+" min.=200 max.=55495580  objects");
    saveAs("TIFF", dirmap+"Object map of_"+list[i]);
    run("3D Manager");
    Ext.Manager3D_AddImage();
    Ext.Manager3D_Measure();
    
    //close();
    close();
   // selectWindow("3D Measure");
   // saveAs("Results", dirdata  + name+"Structure.xls");   
   // selectWindow(name+"Structure.xls");
   // run("Close");
   measuresave(dirdata, name+"Structure.xls");
    Ext.Manager3D_CloseResult("A");
    Ext.Manager3D_Close();
    run("Close All");
  }
  setBatchMode(false);
}

//main function for lumen preprocessing 
function prelumen(dirsource,dirprelumen) {
  list = getFileList(dirsource);
  setBatchMode(true);
  for (i=0; i<list.length; i++) {
    showProgress(i+1, list.length);
    open(dirsource+list[i]);
    name=getTitle();
    attenuation();
    //attenuation();
    //run("Gaussian Blur 3D...", "x=1.5 y=1.5 z=0.5");
    run("Subtract Background...", "rolling=50 sliding stack");
    run("Duplicate...", "title=Duplicate duplicate range slices=1-"+nSlices);
    run("Gaussian Blur 3D...", "x=2 y=2 z=1");
    run("Maximum 3D...", "x=10 y=10 z=1");
    run("Gaussian Blur 3D...", "x=2 y=2 z=1");
    run("Add...", "value=20 stack");	
    imageCalculator("Divide create 32-bit stack", name,"Duplicate");
    close("Duplicate");
    close(name);
    maxx=0;
    minn=10000;
    for (j=1; j<nSlices+1; j++) {
      setSlice(j);
      getRawStatistics(nPixels, mean, min, max, std, histogram);
      minn=minOf(minn,min);
      maxx=maxOf(maxx,max);
    }
    setMinAndMax(minn, maxx);
    run("8-bit");
    saveAs("TIFF", dirprelumen+list[i]);
    run("Close All");
  }
  setBatchMode(false);
}
//function for the user to select the thresholds for the lumens
function lumenmanual(dirprelumen,dirdata) {
  list = getFileList(dirprelumen);
  for (i=0; i<list.length; i++) {
    showProgress(i+1, list.length);
    open(dirprelumen+list[i]);
    name=getTitle();
    if (isOpen("Threshold")) 
      selectWindow("Threshold"); 
    else 
      run("Threshold...");
    setAutoThreshold("Default dark stack");
    waitForUser("Adjust lower limit in Threshold Window");
    getThreshold(lower, upper);
    selectWindow(list[i]); 
    run("Close");
    setResult("name", nResults, name);
    setResult("threshold", nResults-1, lower);
  }
  updateResults();
  selectWindow("Results"); 
  saveAs("Results", dirdata  + "lumenthres.csv");
  run("Close"); 
}


//final processing of lumens in tier 3
function lumenfinal(dirprelumen,dirdata,basic,dirmap,chull) {
  list = getFileList(dirprelumen);
  open(dirdata+"lumenthres.csv");
  thres= newArray(nResults);
  for (i=0; i<nResults; i++) {
    thres[i]=getResult("threshold",i);
  }
  selectWindow("Results"); 
  run("Close"); 
  if(!basic){
    if(chull){
      run("3D Manager Options", "volume surface compactness fit_ellipse 3d_moments convex_hull integrated_density mean_grey_value std_dev_grey_value feret minimum_grey_value maximum_grey_value centroid_(pix) centroid_(unit) distance_to_surface centre_of_mass_(pix) centre_of_mass_(unit) bounding_box radial_distance surface_contact distance_between_centers=10 distance_max_contact=1.80");
    }else{
      run("3D Manager Options", "volume surface compactness fit_ellipse 3d_moments integrated_density mean_grey_value std_dev_grey_value feret minimum_grey_value maximum_grey_value centroid_(pix) centroid_(unit) distance_to_surface centre_of_mass_(pix) centre_of_mass_(unit) bounding_box radial_distance surface_contact distance_between_centers=10 distance_max_contact=1.80");
    }
    
  }
  setBatchMode(true);
  for (i=0; i<list.length; i++) {
    showProgress(i+1, list.length);
    open(dirprelumen+list[i]);
    name=getTitle();
    run("Duplicate...", "title=Duplicate0 duplicate range slices=1-"+nSlices);
    setThreshold(thres[i], 255);
    setOption("BlackBackground", true);
    run("Convert to Mask", "method=Default background=Dark black");
    
    run("Duplicate...", "title=Duplicate duplicate range slices=1-"+nSlices);
    run("Options...", "iterations=1 count=1 black edm=Overwrite do=Nothing");
    run("Fill Holes", "stack");
    imageCalculator("Subtract stack", "Duplicate","Duplicate0");
    run("Gaussian Blur 3D...", "x=1 y=1 z=1");
    setSlice(1);
    setPixel(1, 1, 255);
    run("3D Objects Counter", "threshold=40 min.=0 max.=51961984 objects");
    run("3D Manager");
    Ext.Manager3D_AddImage();
    Ext.Manager3D_Measure();
    measuresave(dirmap,"temp.csv");
    //selectWindow("3D Measure");
    //IJ.renameResults("3D Measure","Results");
    open(dirmap+"temp.csv");
    File.delete(dirmap+"temp.csv");
    //IJ.renameResults("temp.csv","Results");
    open(dirmap+"Object map of_"+list[i]);
    run("Maximum 3D...", "x=1 y=1 z=1");
    for (j=0; j<nResults; j++) {
      setSlice(getResult("CZ (pix)", j)+1);
      setResult("Structure", j, getPixel(getResult("CX (pix)", j), getResult("CY (pix)", j)));
    }
    selectWindow("Results");
    saveAs("Results", dirdata  + name+"lumen.xls");
    run("Close");
    Ext.Manager3D_CloseResult("A");
    Ext.Manager3D_Close();
    run("Close All");
  } 
  setBatchMode(false);
}

//function for skeleton analysis
function skeletons(dirmap,dirtree,dirdata) {
  list = getFileList(dirmap);
  setBatchMode(true);
  for (i=0; i<list.length; i++) {
    showProgress(i+1, list.length);
    open(dirmap+list[i]);
    
    name=getTitle();
    
    run("Duplicate...", "title=Duplicate duplicate range slices=1-"+nSlices);
    setAutoThreshold("Default dark");
    //run("Threshold...");
    setThreshold(1, 1000000);
    setOption("BlackBackground", true);
    run("Convert to Mask", "method=Default background=Dark black");
    //run("Maximum 3D...", "x=1 y=1 z=1");
    // run("Gaussian Blur 3D...", "x=1 y=1 z=1");
    setThreshold(25, 255);
    setOption("BlackBackground", true);
    run("Convert to Mask", "method=Default background=Dark black");
    run("Options...", "iterations=1 count=1 black edm=Overwrite do=Nothing");
    run("3D Fill Holes");
    run("Skeletonize (2D/3D)");
    saveAs("TIFF", dirtree+"Tree of_"+list[i]);
    run("Analyze Skeleton (2D/3D)", "prune=[shortest branch] ");
    selectWindow("Results");
    saveAs("Results", dirdata  + name+"br2.xls");
    selectWindow("Results");
    IJ.deleteRows(0, nResults);
    run("Close");
  //workaround imagej bug (results do not appear sometimes)  
  while( !isOpen("Branch information")){

 	 run("Analyze Skeleton (2D/3D)", "prune=[shortest branch] show");
 selectWindow("Results");
 IJ.deleteRows(0, nResults);
 run("Close");
 	}
    IJ.renameResults("Branch information","Results");
    selectWindow(name);
    run("Maximum 3D...", "x=1 y=1 z=1");
    getVoxelSize(vwidth, vheight, vdepth, unit);
    for (j=0; j<nResults; j++) {
      setSlice(getResult("V1 z", j)/vdepth+1);
      setResult("Structure", j, getPixel(getResult("V1 x", j)/vwidth, getResult("V1 y", j)/vheight));
    }
    selectWindow("Results");
    saveAs("Results", dirdata  + name+"br1.xls");
    run("Close");
    run("Close");
    run("Close");
    run("Close All");
    while( isOpen("Results")){run("Close");}
  }
  setBatchMode(false);
}

//function for threshold analysis
function boni(dirmap,dirdata) {
  list = getFileList(dirmap);
  setBatchMode(true);
  for (i=0; i<list.length; i++) {
    showProgress(i+1, list.length);
    open(dirmap+list[i]);
    
    name=getTitle();
    run("Duplicate...", "title=Duplicate duplicate range slices=1-"+nSlices);
    setAutoThreshold("Default dark");
    //run("Threshold...");
    setThreshold(1, 1000000);
    setOption("BlackBackground", true);
    run("Convert to Mask", "method=Default background=Dark black");
    run("Collect Garbage");
    run("Particle Analyser", "surface_area thickness ellipsoids min=0.000 max=Infinity surface_resampling=2 surface=Gradient split=0.000 volume_resampling=1 labelling=Mapped slices=2");
    run("Collect Garbage");
    selectWindow(name);
    getVoxelSize(vwidth, vheight, vdepth, unit);
    bfr=0;
    for (j=0; j<nResults; j++) {
      setSlice(floor(getResult("z Cent ("+unit+")", j)/vdepth)+1);
      str=getPixel(1+floor(getResult("x Cent ("+unit+")", j)/vwidth), 1+floor(getResult("y Cent ("+unit+")", j)/vheight));
      if (str==0)
      {str=bfr+1;}
      setResult("Structure", j, str);
      bfr=str;
    }
    selectWindow("Results");
    saveAs("Results", dirdata  + name+"boni.xls");
    run("Close");
    run("Close");
    run("Close All");
  }
  setBatchMode(false);
}


//interface part

basic = true;
thres = 30;
check = true;
chull = false;
prelu = true;
skel = true;
bonibonej = false;
luman = true;
luan = true;
items=newArray("None","Enhance boundaries","Split neighbors");

while (check){
  check = false;
  Dialog.create("SAMA-Images");
  
  Dialog.addMessage("TIER 1, automatic (this might take a while)       ");
  Dialog.addCheckbox("Basic morphometrics", basic);
  Dialog.addChoice("Additional filtering", items,"Enhance boundaries")
  Dialog.addNumber("       Threshold for basic morphometrics:", thres);
  Dialog.addCheckbox("Check threshold for basic morphometrics:", check);
  Dialog.addCheckbox("Compute convex hull (optional)", chull);
  Dialog.addCheckbox("Lumen measurement - pretreatment", prelu);
  Dialog.addCheckbox("Branching measurement", skel);
  Dialog.addCheckbox("Thickness measurement (see documentation)", bonibonej);
  
  
  Dialog.addMessage("TIER 2, manual");
  Dialog.addCheckbox("Lumen measurement - input thresholds", luman);
  
  Dialog.addMessage("TIER 3, automatic");
  Dialog.addCheckbox("Lumen measurement - final step ", luan);
  
  
  Dialog.show();
  basic = Dialog.getCheckbox();
  mode = Dialog.getChoice();
  thres = Dialog.getNumber();
  check = Dialog.getCheckbox();
  chull = Dialog.getCheckbox();
  prelu = Dialog.getCheckbox();
  skel = Dialog.getCheckbox();
  bonibonej = Dialog.getCheckbox();
  luman = Dialog.getCheckbox();
  luan = Dialog.getCheckbox();
  
  if(check){
    adress=File.openDialog("Choose original file to test");
    setBatchMode(true);
    open(adress);
    preprocessBasic(mode);
    setBatchMode(false);
    if (isOpen("Threshold")) 
      selectWindow("Threshold"); 
    else 
      run("Threshold...");
    setThreshold(thres, 255);
    waitForUser("Adjust lower limit in Threshold Window");
    getThreshold(thres, upper);
    close();
  }
  
}


//create required directories (if needed)

dirsource = getDirectory("Choose Source Directory");
dir=File.getParent(dirsource)+"/";
dirmap=createdir(dir,"SAMA-Object map");
dirdata=createdir(dir,"SAMA-Image data");
dirtree=createdir(dir,"SAMA-Tree analysis");
dirprelumen=createdir(dir,"SAMA-Prelumen images");


//actual processing

if(basic){
  print("Measuring structures...");
  mainstructures(dirsource, dirmap,dirdata,thres,chull,mode);
  
}
if(skel){
  print("Measuring skeletons...");
  skeletons(dirmap,dirtree,dirdata);
  
}
if(bonibonej){
  print("Measuring thickness from Bonej...");
  boni(dirmap,dirdata);
}

if(prelu){
  print("Prepocessing lumens...");
  prelumen(dirsource,dirprelumen);	  
}
if(luman){
  print("Manual choice of threshold...");
  lumenmanual(dirprelumen,dirdata);		
}
if(luan){
  print("Measuring lumens...");
  lumenfinal(dirprelumen,dirdata,basic,dirmap,chull);		    
}
print("Sama images completed, please run SAMA-Analyze");
print("***************************");