//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 'main [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. **' );
  }

//----------------------------------------------------------------------------------------------------------
function measureSlices() { requires('1.31o');
//  the format of the array is RESULT[SLICE-1,ANNULUS-1]

  if (bitDepth() == 24) exit("bw image required");
  
  bgsize = 10;                                                 //how many pixels away from the traced cell edge the outer edge of the background annulus is
  bginset = 4;                                                  //how many pixels away from the traced edge of the cell the inner edge of the background annulus is                
  slices = nSlices();                                                          // the number of images in a stack to be measured
  if (slices < 1) exit ('No image chosen');
  annuli = 40;                                                                    // the maximum number of annuli that could be measured
  //annulusmean = newArray (slices, annuli);                // where the mean of each annulus for each slice is stored
  //Whoops, 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
  bg = newArray(slices);                                                   // for each slice, the calculated background for subsequent math
  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]
  
  setupMeasurementOptions();                                      //set up measurement parameters, for instance: mean, area, decimals...
  
  //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");

  
  // 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

  // Get the background from the background annulus.  The tracing steps out bgsize number of pixels.
  // The annulus is bgsize-bginset number of pixels wide.
  expandselection(bgsize);
  tempbg = newArray(slices);
  for (i=1;i<=slices;i++){
    run("Set Slice...", "slice="+i);
    run("Measure");
    tempbg[i-1] = getResult("Mean", (nResults-1));
    }
  temparea=getResult("Area", (nResults-1));
  // now get the inner edge of the the background annulus and calculate the backgrounds to be used later
  shrinkselection(bgsize-bginset);
  for (i=1;i<=slices;i++){
    run("Set Slice...", "slice="+i);
    run("Measure");
    bg[i-1] =  ((tempbg[i-1]*temparea)-(getResult("Mean", (nResults-1)) * getResult("Area",(nResults-1)))) / (temparea-getResult("Area",(nResults-1)));
    //print(bgsize+' '+bginset+'  '+bg[i-1]);    
    }
  shrinkselection(bginset);     // return to the original tracing;  a better method would be to have the original tracing saved and reloaded
  run("Clear Results"); 

  for (i=1;i<=slices;i++) cellmean[i-1] = cellmean[i-1] - bg[i-1];       // corrects the whole cell mean for background

// 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);

  innermean = newArray(slices);                 // the  mean value of the inner annulus 
  innermass = newArray(slices);                 // the  total mass of the inner annulus  (the trade off is more memory for fewer calculations)
  outermean = newArray(slices);                 // the mean value of the outer annulus 
  outermass = newArray(slices);                 // the total value of the outer annulus  (the trade off is more memory for fewer calculations)

  //a way to deal with the first outer annulus
  for (i=1;i<=slices;i++){
    run("Set Slice...", "slice="+i);
    run("Measure");
    innerarea = getResult("Area", (nResults-1));  // this is redundant after the first slice, but beats an if statement
    innermean[i-1] = getResult("Mean", (nResults-1)) - bg[i-1];     // do the background subtraction immediately
    innermass[i-1] = innermean[i-1] * innerarea;
    }

  // Do the measurements on all the annuli   
  for (j=1;j<=annuli;j++){
    shrinkselection(1);
    outerarea = innerarea;  // the previous annulus inner ring is now the outer ring
    for (i=1;i<=slices;i++){
      outermean[i-1] = innermean[i-1];        // the previous annulus inner ring is now the outer ring                       
      outermass[i-1] = innermass[i-1];        // the previous annulus inner ring is now the outer ring
      run("Set Slice...", "slice="+i);
      run("Measure");
      innermean[i-1] = getResult("Mean", (nResults-1)) - bg[i-1];     // do the background subtraction immediately
      innerarea = getResult("Area", (nResults-1));         // this is redundant after the first slice, but beats an if statement
      innermass[i-1] = innermean[i-1] * innerarea;
      annulusarea[j-1] = (outerarea - innerarea);
      //annulusmean[(i-1)][(j-1)] = (outermass - innermass) / annulusarea[j-1]) ;
      //print(j + '  ' + annulusarea[j-1] + '  ' + outermass[i-1] + '  ' + innermass[i-1] + '  ' +  ((outermass[i-1] - innermass[i-1]) / annulusarea[j-1])   );
      putNumberIntoTwoDArray(arrayID, i-1, j-1, ((outermass[i-1] - innermass[i-1]) / annulusarea[j-1]));

    }  //for i = 1 to slices

  } //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
  selectImage(arrayID);  setLocation(10,100);   run("Close");  
  run("Clear Results");   
}

//----------------------------------------------------------------------------------------------------------
//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]);
  }