//  Macros for Inbar to measure radial plot profiles
//  Michael Cammer, May 2014


var plot1Name = "protein1";    // Alexa 488
var plot2Name = "protein2";       // OKT3 647
// The following correspond to the channel number in the hyperstack
var plot1 = 1;   //  channel 1
var plot2 = 2    //  channel 4
var maxPlotLength = 96;     // a number bigger than the largest possible cell size
var  singleProfile = newArray(maxPlotLength);
var  plot1output = newArray(maxPlotLength);   
var  plot2output = newArray(maxPlotLength); 
var  tab = "\t";
var meanCell = -1;  // used by macro "multichannel Integrated intensity through time"

//================================================================
// Radial profile plot.
macro "Radial Profile Plot multiple cells from overlays [q]" {
  requires("1.49a");              //  only run if the correct version of ImageJ is running
  openImages = nImages();

  getVoxelSize(width, height, depth, unit);     // initialize the arrays as strings
  for (i=0; i<maxPlotLength; i++) { 
    plot1output[i] = "" + (width*i) + tab + plot1Name + tab ; 
    plot2output[i] = plot2Name + tab ; 
  }   
  
  for (i=1; i<= openImages; i++) {   
    selectImage(i);
    run("Set Scale...", "distance=1 known=1 pixel=1 unit=px");   // set units to pixels
    run("Select None");          // make sure there are no ROIs on the current image
    roiManager("reset");        // clear all ROIs from the last image mesured
    run("To ROI Manager");   //  convert the overlays to ROIs sothey may be measured one by one
    roiManager("Remove Channel Info");
    roiManager("Remove Slice Info");
    roiManager("Remove Frame Info");

    // Do the measurements here.  This block of code may be replaced with whatever measurements you want.
    for (cell=0; cell<roiManager("count"); cell++) {

       // first measure by radial reslice method
       //set channel plot1
       //call measurment function with ROI #
       //move singleProfile to plot1output
      // set channel plot 2
       //call measurment function with ROI #
       //move singleProfile to plot2output

       // measure by increasing size bands method to compare to radial method
       Stack.setChannel(plot1);
       //setSlice(plot1);
       radialMeasureUsingBands(cell);
       for (k=0; k<maxPlotLength; k++) plot1output[k] = plot1output[k] + singleProfile[k] + tab;

       Stack.setChannel(plot2);
       //setSlice(plot2);
       radialMeasureUsingBands(cell);
       for (k=0; k<maxPlotLength; k++) plot2output[k] = plot2output[k] + singleProfile[k] + tab;

       //call measurment function with ROI #
       //move singleProfile to plot1output
      // set channel plot 2
       //call measurment function with ROI #
       //move singleProfile to plot2output

    } // for each cell
  }  // for i - process each open each image

  // report results
  print (getTitle());
  for (i=0; i<maxPlotLength; i++) { 
    print (plot1output[i] + plot2output[i]);
  }   
  selectWindow("Log");
}  // end Measure multiple cells from overlays






//============================================================
function radialMeasureUsingBands(ROInum) {
  roiManager("select", ROInum);
  run("Set Measurements...", "  mean centroid redirect=None decimal=1");
  run("Measure");
  xc = getResult("X", nResults-1);
  yc = getResult("Y", nResults-1);
  run("Clear Results");

  for (k=0; k<maxPlotLength; k++) {
    makeOval(xc-(k/2), yc-(k/2), k+1, k+1);
    if (k > 0) run("Make Band...", "band=1");
    run("Measure");
    singleProfile[k] = getResult("Mean", nResults-1);
  }  // end of expanding bands

} // end radialMeasureUsingBands(ROInum)


//============================================================
//  This macro requires a multichannel multi time series in a hyperstack.
//  To work properly, there need to be channels > 1, slices = 1, and frames > 1.
//  Frames is the time dimension.  To use on Z instead, the user would have to swap frames for slices.
//============================================================
macro "multichannel Integrated intensity through time" {
  Stack.getDimensions(width, height, channels, slices, frames);
  //  next two lines may be changed for different #s of channels and different proteins to be measured
  channelstomeasure = newArray(2, 3);                      // channels to measure
  channelsnames = newArray ("GP120", "OKT3");   //  names of channels to measure
  run("Set Measurements...", "  mean standard redirect=None decimal=3");
  output = newArray(frames+1);    // top row is roi #
  for (t=1; t<=output.length; t++) output[t-1] = "";

  for(roi=0; roi<roiManager("count"); roi++) {   // only use XY locations to define ROIs
    roiManager("select", roi);
    roiManager("Remove Channel Info");  
    roiManager("Remove Slice Info");
    roiManager("Remove Frame Info");
  }

  setBatchMode(true);  newImage("Untitled", "8-bit ramp", 10, 10, 1);  close();
  for (c=0; c<channelstomeasure.length; c++) {
    for(roi=0; roi<roiManager("count"); roi++) {
      for (t=1; t<=frames; t++) {
        roiManager("select", roi);
        Stack.setChannel(channelstomeasure[c]);         // this is here in case the ROI includes the wrong channel info
        Stack.setFrame(t);
           run("Enlarge...", "enlarge=5 pixel");
           run("Make Band...", "band=1");
           run("Measure");
           min = getResult("Mean", nResults()-1) + ( getResult("StdDev", nResults()-1) *2);        
        roiManager("select", roi);
        Stack.setChannel(channelstomeasure[c]);         // this is here in case the ROI includes the wrong channel info
        Stack.setFrame(t);
        
        getHistogram(values, counts,16384, 0, 16383);
        mass = caclulateIntegratedDensity(counts, min);
        if (t == 1) output[0] = output[0] + "\t " + channelsnames[c] + "_ROI#" + (roi+1) + "\t mean" ; 
        output[t] = output[t] + "\t " + mass + "\t " + meanCell ;    
      }  // for t
    }  // for roi
  }  // for c
  setBatchMode(false);
  Array.show("Integrated Intensities", output);
}  // end macro multichannel Integrated intensity through time


//===============================================================
macro "test intden" {
   for (test=0; test < 4; test++) {
        roiManager("select", test);
        run("Make Band...", "band=1");  
        run("Measure");
        min = getResult("Mean", nResults()-1) + ( getResult("StdDev", nResults()-1) ); 
        roiManager("select", test);
        getHistogram(values, counts,16384, 0, 16383);
        print(min);
        mass = caclulateIntegratedDensity(counts, min);
        print(mass);
        selectWindow("Log");
  } // test
}

function caclulateIntegratedDensity(counts, min) {
  intden = 0;
  area = 0;
  for (i=min; i<counts.length; i++) {
    intden = ((i - min) *  counts[i]) + intden;  
    if (i > min) area = area + counts[i];
  }  // for i
  meanCell = intden / area;  // this was added later so it is a global variable instead of a return
  return intden;
}