/*  macro for  ImageJ 1.50 by Michael Cammer - please cite NYULMC Microscopy Core if data published

Macro for colocalization (Pearson's correlation) at thresholds set by the user.  Pearson's correlation was coded in 2004 for use and
distribution by Analytical Imaging Facility at Albert Einstein College of Medicine.  Wrappers for ease of use changed since for many
diffferent laboratories.

v105 20150821 for  requires a RGB image.
          The image is converted to Composite and two channels in the array "compare" are corrleated.
          Multiple regions of interest are in the ROI Manager and the correlation of each is listed in the Log window.

          There's a bug that if there are too many zeroes, NaN is calculated.  Haven't figured out why.  But add 1 globally to images,
          at least for the simulations, before running to get rid of this problem and shift the resut only minimally higher.

Warning:  spatial scale is removed from the image.

Instructions:
  User selects regions of interest.  Each one is put in ROI Manager.  
  Recommend saving ROI Manager if need to do analysis on same regions again.
  (Regions selections may be sone on multiple files before doing actual analysis.)
  If opening ROI Manager files from disk, make sure previous one closed first.
  User runs macro.
  Save log window as text file or copy-paste into Excel or other spreadsheet.



*/

var thresholds = newArray(8, 16, 32, 48, 64, 80, 128);
var median = true;        // change to false if don't want to median filter images
var compare = newArray(1 ,2 );  // The channels to compare, in this case green and blue


//**************************************************************
//Macros to implement Pearson's coefficient analysis of digital images.      
//michael.cammer@med.nyu.edu
//Pearson's section written and tested  10 March 2004.
//Additions written October 2010.
//
//**************************************************************









//=================================================
/*  This macro assumes image is and RGB TIFF file.
     It is converted to a composite stack.


*/


/*  convert RGB to composite;
    user will trace multiple regions and put in ROI manger;
   analyze each region
*/
macro "Coloc on single image multiple ROIs [q]" {
  requires("1.50");
  close("\\Others");
  run("Set Scale...", "distance=1 known=1 unit=pixels");
  title = getTitle();
  run("Split Channels");
  run("Merge Channels...", "c1=[" + title + " (red)] c2=[" + title + " (green)] c3=[" + title + " (blue)] create");
  rename(title);
  outputString = "filename "+ " \t " + "0" 
  for (i=0; i<thresholds.length; i++)
    outputString = outputString + " \t " + thresholds[i];
  printDateTime();
  print(outputString);
  for (roi=0;  roi<(roiManager("count")); roi++){
    convertIToArrays(roi);
  }
  selectWindow("Log");
}



//***********************************************************************
// TheROIs are converted to 1D arrays that can be compared directly without graphics calls.
//***********************************************************************
function convertIToArrays(roi) {
  run("Set Measurements...", "area bounding fit redirect=None decimal=4");
  roiManager("select", roi);
  run("Measure");
  BX = getResult("BX", nResults()-1);
  BY = getResult("BY", nResults()-1);
  rectWidth = getResult("Width", nResults()-1);
  rectHeight = getResult("Height", nResults()-1);
  pixels = rectWidth * rectHeight;
  a = newArray(pixels);
  b = newArray(pixels);
  
  for (chan=0; chan<2; chan++) {
    arrayindex = 0;
    Stack.setChannel(compare[chan]) ;
    for (x=BX; x<BX+rectWidth; x++)
      for (y=BY; y<BY+rectHeight; y++){
        if (selectionContains(x, y)) {       // recode by copying entire image and masking instead of this command, faster?
          if (chan==0)
            a[arrayindex] = getPixel(x, y);
          else
            b[arrayindex] = getPixel(x, y);
          arrayindex++;
          //putPixel(x,y,128);
        }  // if point inside ROI
      } // end for x and y loops

  } // for each channel 
  coef = calculateCoefficient(a, b);  
  outputString =getTitle() + " \t " + coef + " \t ";

  for (t=0; t<thresholds.length; t++) {
     for (z=0; z<a.length; z++) {
       if ( (a[z]<thresholds[t]) || (b[z] < thresholds[t])  ) { a[z] = 0;  b[z] = 0; }
    }
    coef = calculateCoefficient(a, b);  
    outputString = outputString + coef + " \t ";
  }  // for t, report at each threshold

 

  print (outputString);
  //print(f, outputString);

}  //end function convertSlicesTo1Darrays()


//*****************************************************************
// function calculateCoefficient(a, b);
// The function is passed two arrays of integers.  The arrays must be of
// equal length.  They are compared using Pearson's coeffiicent.
// Values of Zero are not considered.
//*****************************************************************
function calculateCoefficient(a, b) {
  if (a.length != b.length) exit ('arrays of unequal length');
  n = a.length;
  sXY=0;  sX=0;  sY=0;  sX2=0;  sY2=0; nn = 0;
  for (i=0; i<n; i++){
    if ((a[i] != 0) && (b[i] != 0)) {
      sXY = sXY + (a[i] * b[i]);
      sX = sX + a[i];
      sY = sY + b[i];
      sX2 = sX2 + (a[i] * a[i]);
      sY2 = sY2 + (b[i] * b[i]);
      nn++;
    }  // if
  } // for
  top = sXY - ((sX * sY) / nn);
  bottomleft = sX2 - ((sX * sX) / nn);
  bottomright = sY2 - ((sY * sY) / nn);
  bottom = pow((bottomleft*bottomright),0.5);
  Pearson = top / bottom;
   
  //print('sXY='+sXY+'  sX='+sX+'  sY='+sY+'  sX2='+sX2+'  sY2='+sY2+'  n='+n+' nn='+nn+'  top='+top+'  Pearson='+Pearson); 
  //print(Pearson);
  //selectWindow("Log");
  return Pearson;
} //calculateCoefficient



//===============================================================
// This macro demonstrates how to use the getDateAndTime() 
// function, available in ImageJ 1.34n or later.

function printDateTime() {
     MonthNames = newArray("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec");
     DayNames = newArray("Sun", "Mon","Tue","Wed","Thu","Fri","Sat");
     getDateAndTime(year, month, dayOfWeek, dayOfMonth, hour, minute, second, msec);
     TimeString ="Date: "+DayNames[dayOfWeek]+" ";
     if (dayOfMonth<10) {TimeString = TimeString+"0";}
     TimeString = TimeString+dayOfMonth+"-"+MonthNames[month]+"-"+year+"   Time: ";
     if (hour<10) {TimeString = TimeString+"0";}
     TimeString = TimeString+hour+":";
     if (minute<10) {TimeString = TimeString+"0";}
     TimeString = TimeString+minute+":";
     if (second<10) {TimeString = TimeString+"0";}
     TimeString = TimeString+second;
     print(TimeString);
  }