//Michael Cammer  cammer@aecom.yu.edu  http://www.aecom.yu.edu/aif/
//Macros written for annulus analysis of fluorescent labeled cells.
//rewritten from three previous versions January 2004
//Reasons rewritten:
//1.  NIH-Image not suited to various image formats being used at AECOM
//2.  Other ImageJ versions too hard to edit because too little documentation.
//3.  Older versions specific for number of flour probes.
//4.  Clean up structure of precessing to a series of 2D arrays and output to clean table.
//5.  Build a library or routines.

//A major improvement would be to save the traced selection before calculating the background and
//restoring this selection after calculating the background.
//This would improvethe spatial resolution at the very edge.

macro 'Measure annuli by edge count method [q]'{
  starttime = getTime();
  if (nSlices==1) run("Add Slice");  //This is a truly disgusting way to make this work for
                      // single image analysis
  measureSlices();
  endtime = getTime();  elapsed = (endtime - starttime) / 1000;
  print('** This operation took '+elapsed+' seconds. **' );
  }
//----------------------------------------------------------------------------------------------------------
macro 'Project and Threshold [b]' { requires('1.31u');
  if (nSlices() !=2) exit('You need exactly 2 registered slices in the stack.');
  run("Z Project...", "start=1 stop=2 projection='Max Intensity'");
  run("Threshold...");
  setTool(8);
}

//----------------------------------------------------------------------------------------------------------
function measureSlices() { requires('1.31o');
//  the format of the array is RESULT[SLICE-1,ANNULUS-1]
  if (bitDepth() == 24) exit("bw image required");
  bgsize = 6;                                                                       //how many pixels away from the traced cell edge the background annulus is 
  slices = nSlices();                                                          // the number of images in a stack to be measured
  if (slices < 1) exit ('No image chosen');
  bg = newArray(slices);                                                  // the background value for each slice
  annuli = 40;                                                                     // the maximum number of annuli that could be measured
  //The macro language doesn't support 2D arrays, so here is the work around
    originalpic = getImageID();
    arrayID = makeTwoDArray('annulusmean', slices, annuli);
    selectImage(originalpic);
  annulusarea = newArray (annuli);                               // the area will be the same for all slices
  cellarea=0;                                                                       // the total cell area is the same for all slices
  cellmean = newArray(slices);                                      // for each slice the entire cell as originall traced
  // the element in the array corresponding to the slice is [slicenumber-1]
  
  run("Clear Results");
  run("Set Measurements...", "area mean  redirect=None decimal=5");
  
  //measure the whole cell as traced
  for (i=1;i<=slices;i++){
    run("Set Slice...", "slice="+i);
    run("Measure");
    cellmean[i-1] = getResult("Mean", (nResults-1));
    }
  cellarea=getResult("Area", (nResults-1));
  run("Clear Results");

//APPROX LINE 54

  // Make a new window for the mask
  // The mask window will be used overn and over to reduce processing time.
  run("New...", "name=RESIZEMASK type=8-bit fill=White width="+getWidth()+" height="+getHeight());
  run("Restore Selection");  
  setLocation(1, 1280);            // moves the window off the screen to hide it
  selectImage(originalpic);
  setLocation(1, 1280);            // moves the window off the screen to hide it

  run("Line Width...", "line=1");
  // Get the background from the background annulus.  The tracing steps out bgsize number of pixels.
  expandselection(bgsize);  
  getSelectionCoordinates(x, y);             // the edges of the ROI are put in arrays x and y
  makeSelection("freeline", x, y);             // convert to freeline so that "measure" will only measure the edge values
  for (i=1;i<=slices;i++){                            // measure the edge values for each slice
    run("Set Slice...", "slice="+i);
    run("Measure");
    bg[i-1] = getResult("Mean", (nResults-1));
    cellmean[i-1] = cellmean[i-1] - bg[i-1];             // corrects the whole cell mean for background
  }
  makeSelection("freehand", x, y);                         // restores the tracing

  shrinkselection(bgsize);     // return to the original tracing;  a better method would be to have the original tracing saved and reloaded

  // Printout the results for the whole cell
  print(getTitle());
  output='area  ';
  for (i=1;i<=slices;i++) output = output+'meanslice'+i+'  ';
  print(output);
  output = '' + cellarea + '             ';
  for (i=1;i<=slices;i++) output = output+  cellmean[i-1]  +'  ';
  print(output);

  // Do the measurements on all the cell annuli   
  for (j=1;j<=annuli;j++){  
    getSelectionCoordinates(x, y);             // the edges of the ROI are put in arrays x and y
    annulusarea[j-1] = x.length;
    makeSelection("freeline", x, y);
    for (i=1;i<=slices;i++){                           // measure the background for each slice   bg[i-1] holds the value for slice i
       run("Set Slice...", "slice="+i);
       run("Measure");
       annulusmean = getResult("Mean", (nResults-1)) - bg[i-1];
       putNumberIntoTwoDArray(arrayID, i-1, j-1, annulusmean);       
     }  // end for i 1 to slices
     makeSelection("freehand", x, y);   // restores the tracing
     shrinkselection(1);
  } // end for j 1 to annuli


  // Printout the results
  print(getTitle());
  output='annulus area ';
  for (i=1;i<=slices;i++){
    output = output+'slice'+i+'  ';
  }
  print(output);

  // note:  future versions should print out the mass of each annulus and the running total of percent of protein mass from
  // edge to current location

  for (j=0;j<annuli;j++){
    output='' + j + '                 ' + annulusarea[j] + '  ';
    for (i=0;i<slices;i++){
      number = gettNumberFromTwoDArray(arrayID, i, j);
      output = output + round(number) + '       ';
    }  //for i = 0 to slices
    print(output);
  } //for j = 0 to annuli
  print('');  // blank line to separate from next output

  // Cleanup the mask
  selectWindow("RESIZEMASK"); run("Close");  

  // Get back the hiden image window
    selectImage(originalpic);
    setLocation(10, 48);            // moves the window off the screen to hide it

  // Cleanup the 2D array and the results window
  selectImage(arrayID);  setLocation(10,100);   run("Close");  
  selectWindow("Results");  run("Close");
}

//----------------------------------------------------------------------------------------------------------
//pid is the unique pointer to the 2D array image window
//x,y is the location where the value, number, should go
function putNumberIntoTwoDArray(pid, x, y, number) {requires('1.31o');
  current = getImageID();
  selectImage(pid);
  putPixel(x,y,number);
  selectImage(current);
}

//----------------------------------------------------------------------------------------------------------
//pid is the unique pointer to the 2D array image window
//x,y is the location where the value, number, is located
function gettNumberFromTwoDArray(pid, x, y) {requires('1.31o');
  current = getImageID();
  selectImage(pid);
  number = getPixel(x,y);
  selectImage(current);
  return number;
}

//----------------------------------------------------------------------------------------------------------
function makeTwoDArray(windowname, width, height) {requires('1.31o');
  run("New...", "name="+windowname+" type=32-bit fill=Black width="+width+" height="+height+" slices=1");
  setLocation(1, 1280);     // moves the window off the screen to hide it
  return getImageID();
}

//----------------------------------------------------------------------------------------------------------
function setupMeasurementOptions(){
  run("Set Scale...", "distance=1 known=1 pixel=1 unit=");
  run("Set Measurements...", "area mean decimal=5");
  run("Grays");                                                                         // get rid of the threshold mask color so image will be measured
  run("Clear Results");
  }

//-----------------------------------------------------------------------------------
  function  expandselection(n) {
      requires("1.30i");
      if (selectionType==-1)
          exit("Selection required");
      resize(n);
  } //function  expandselection(n)

//-------------------------------------------------------------------------------------
  function shrinkselection(n) {
      if (selectionType==-1)
          exit("Selection required");
      resize(-n);
  } //function shrinkselection(n)

//-------------------------------------------------------------------------------------
  function resize(n) {
      type = selectionType();
      if (type==0) {
          getBoundingRect(x, y, width, height);
          makeRectangle(x-n, y-n, width+n*2, height+n*2);
     } else if (type==1) {
          getBoundingRect(x, y, width, height);
          makeOval(x-n, y-n, width+n*2, height+n*2);
     } else
          resizeMask(n);
  }  //function resize(n)

//-------------------------------------------------------------------------------------
function saveMask() {
      snapshot();
      id = getImageID();
      run("New...", "name=OriginalMask type=8-bit fill=White width="+getWidth()+" height="+getHeight());
      run("Restore Selection");
      setColor(0);
      fill();
      run("Select None");
      selectImage(id);
}

//-------------------------------------------------------------------------------------
  function resizeMask(n) {
      if (selectionType==-1) return;
      id = getImageID();
      //run("New...", "name=Temp type=8-bit fill=White width="+getWidth()+" height="+getHeight());
      selectWindow("RESIZEMASK"); 
      run("Restore Selection");
      setColor(0);
      fill();
      run("Select None");
      if (n>0)
         expand = true;
      else
          expand = false;
      iterations = abs(n);
      run("Options...", "iterations=1 ");  // I wrote my own iteration scheme because of a bug in 1.31o
      if (expand)
          for (a=1;a<=iterations;a++) run("Dilate");
      else
          for (a=1;a<=iterations;a++) run("Erode");
      x2=-1; y2=-1;
      getBoundingRect(xbase, ybase, width, height);
      inc = 1;
      if (width>5 && height>5)
          inc = 5;
      for (x=xbase; x<xbase+width; x+=inc) {
          for (y=ybase; y<ybase+height; y+=inc) {
              if (getPixel(x,y)==0) {
                  x2= x;
                  y2 = y;
              }
          }
      }
      if (x2==-1)
          print('Selection is too small');
     else
          doWand(x2, y2);
      //run("Close");
      selectImage(id);
      run("Restore Selection");
  }   //function resizeMask(n)

//------------------------------------------------------------------------------------------
  function printArray(a) {
      for (i=0; i<a.length; i++)
          print(i+": "+a[i]);
  }