/*

   Modified Feb 2017 for Ammar Zafar to prep the image from a confocal where the intensity values were inverted.


   Macro for ImageJ 1.51h7 written 20161130 by Michael Cammer for Samir Hanna to measure the widths of 
   nanotubes as stained for f-actin and imaged using the OMX SIM microscope at the bioimaging core at
   the Rockefeller University in NYC. 

   Draw line approx 4 um long perpendicular to the nanotube with center near or on the nanotube.  
   Run macro.

   The macro:
      Gets the intensity profile.
      Fits to a Gaussian.
      Normalizes lowest value to 0 and the peak to 1.
      Finds the width at intensity 0.32 which should be between 1 and 2 standard deviations.
      Finds the width at intensity 0.01 which is more than 99.5% certainty wider than or equal to the deconvolved width of the f-actin.  
      Reports both answers.
*/

var tab = " \t";


macro "Prep Bacteria Image" {
  run("Bandpass Filter...", "filter_large=256 filter_small=0 suppress=None tolerance=5");
  run("Invert", "slice");
  selectWindow("Recorder");
  run("Invert LUT");
}




macro "Measure width [q]" {
  title = getTitle();
  path = getDirectory("image");
  run("Set Measurements...", "area redirect=None decimal=5");
  profile = getProfile();
  run("Measure");
  len = getResult("Length", nResults()-1);
  pixelsize = len / profile.length;   // each pixel in um
  xpoints = newArray(profile.length);    // will be distance in um
  for (i=0; i<profile.length; i++)
    xpoints[i] = i * pixelsize;                    // populate linearly with distance for curve fitting
  Fit.doFit("Gaussian", xpoints, profile);
  Fit.plot;
  run("Save", "save=["+path+"Gaussianfit_"+title+"]");
  close();
  if (Fit.rSquared < 0.8) exit("not robust enough:  r squared is less than 0.8");
  fitresult = newArray(profile.length);       // array of fit values
  for (i=0; i<profile.length; i++)                 // get the fitted results
    fitresult[i] = Fit.f(xpoints[i]);

//  for (i=0; i<profile.length; i++)                 // show results for debugging
//    print(xpoints[i] + tab + profile[i] + tab + fitresult[i]);

  normalized = newArray(profile.length); 
  Array.getStatistics(fitresult, min, max, mean, stdDev);
  for (i=0; i<profile.length; i++)                 // subtract the min and normalize to 1
    normalized[i] = (fitresult[i] - min) / (max - min);
  //  Array.print(normalized);  // check the math is right

  // Report width at base.-----------------------------------------------------------------------------
  threshold = 0.01;           //  Used to guarranty all photons defining width within 3 SDs of mean or > 99.5%.
  // First point on curve > threshold is considered the width at the base. 
  for (i=0; i<profile.length; i++)       // get left
    if (normalized[i] > threshold) {
      left = xpoints[i];
      i = 999999;          // exits loop; would be cleaner with repeat until, but this works
    }
  for (i=profile.length-1; i>=0; i--)   // get right
    if (normalized[i] > threshold) {
      right = xpoints[i];
      i = -1;                   // exits loop
    }
  //  print (title + "\t Width at base is \t" + (right-left) + "\t um");

  // Report width at half max.--------------------------------------------------------------------------
  threshold = 0.32;            // Used to estimate between 1 and 2 SDs from the mean or between 70% and 95% photons defining width.
  for (i=0; i<profile.length; i++)       // get left
    if (normalized[i] > threshold) {
      left1 = xpoints[i];
      left0 = xpoints[i-1];
      lefthm = (left1 + left0) / 2;
      i = 999999;          // exits loop; would be cleaner with repeat until, but this works
    }
  for (i=profile.length-1; i>=0; i--)   // get right
    if (normalized[i] > threshold) {
      right0 = xpoints[i];
      right1 = xpoints[i+1];
      righthm = (right1 + right0) / 2;        //  really, this should be weighted or fit to a line.
      i = -1;                   // exits loop
    }
  //  print(lefthm, righthm);
  print (title + "\t Width at 1-2 SDs is \t" + (righthm-lefthm) + "\t Width at base is \t" +  (right-left) + "\t um");

  selectWindow("Log");
}  // end macro