/*
Files are 4 channels.  Channel 2 is green.  Channel 3 is red.  They need to be compared.
Step 1.  Spilt the channels into seprate files for JACoP.
Step 2.  Run JACoP
*/

/*  Michael Cammer   20150303
    ImageJ 1.49o30
*/


macro "JACoP on all files in directory [q]" {
  run("Close All");
  //dir1 = getDirectory("Choose Source Directory ");
  inputDir = getDirectory("Choose Source Directory ");
  //outputDir = getDirectory("Choose Destination Directory ");
  outputDir = inputDir; 
  list = getFileList(inputDir);
  setBatchMode(true);
  for (i=0; i<list.length; i++) {
     showProgress(i+1, list.length);  
     open(inputDir+list[i]);
     original = getImageID();
     print("\n\n" + getTitle());
     run("Duplicate...", "title=green duplicate channels=2");
     selectImage(original);
     run("Duplicate...", "title=red duplicate channels=3");
     //selectImage(original);
     //close();
     //exit();
     run("JACoP ", "imga=[green] imgb=[red] thra=16 thrb=16 pearson mm");       //change out threshold numbers
     print("Michael's version r= ");
     convertISlicesTo1Darrays("red", "green", 1);
     run("Close All");
  }
  setBatchMode(false);
}  // end JACOP macro

============
//=======================================================================================
// Pearson's correlation added to stacks macros on 14 June 2011
//=======================================================================================

//**************************************************************
//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 version requires two stacks.  
//**************************************************************

macro "Pearson's coefficient of ROI [F7]" {
  requires('1.32c');
  // This section gets the titles of the images to operate on and converts them to PIDs.
  if (nImages<1) exit ("No image windows are open");   
  titlelist = newArray(nImages);                                      // titles of all open windows
  for (i=1; i<=nImages; i++) {
    selectImage(i);
    titlelist[i-1] = getTitle;
    run("Put Behind [tab]");
  }
  Dialog.create("scatterplot two images")
  Dialog.addMessage("Please pick two images to be compared.\nIf two different images have the same file name, results are unpredictable.");
  Dialog.addChoice("first image", titlelist);
  reversetitlelist = Array.invert(titlelist);
  Dialog.addChoice("second image", reversetitlelist);
  Dialog.show();
  image1 = Dialog.getChoice();
  image2 = Dialog.getChoice();
  print(image1, image2);
  selectWindow(image1);              // get image names and convert them to PIDs because these are more robust for use later.
  image1n = getImageID;
  selectWindow(image2);
  image2n = getImageID;

  for (i=1; i<=nSlices; i++){
    convertISlicesTo1Darrays(image1, image2, i);
  }
}  // end macro Pearson's


//***********************************************************************
// function convertImagesTo1Darrays(image1n, image2n);
// image1 and image2 are the PIDs of the stacks to be compared.
// The images are converted to 1D arrays that can be compared directly without graphics calls.
//***********************************************************************
function convertISlicesTo1Darrays(image1, image2, slice) {
  pixels = getWidth() * getHeight();
  cf = newArray(pixels);
  pip = newArray(pixels);
    arrayindex=0;
    selectImage(image1);  
    setSlice(slice);
    for (x=0; x<getWidth(); x++)
      for (y=0; y<getHeight(); y++){
        cf[arrayindex] = getPixel(x, y);
        arrayindex++;
      } // end for x and y loops
    
    arrayindex=0;
    selectImage(image2);
    setSlice(slice);
    for (x=0; x<getWidth(); x++)
      for (y=0; y<getHeight(); y++){
        pip[arrayindex] = getPixel(x, y);
        arrayindex++;
      } // end for x and y loops
     
  coef = calculateCoefficient(cf, pip);  
}  //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