// Michael Cammer; Michael Dustin laboratory at Skirball at NYU Medical    Dec 2011

// These macros were written specifically for Catarina Sacristan's multi-field timelapse experiments 
// on the Nikon TIRF microscope run by Elements software with files saved in .nd2 format.
// The files are imported using LOCI BioFormats importer with all series selected and channels split.

var number_of_radii = 2500;   // number of pixels on perimeter of circle to be sampled.
var clipsize = 150;  // edge size of the box for clipping cells for montage
var snipwindow = -1;
var projectiondepth = 16;  // default depth for projection of variable depth
var tab = " \t";


// Should be updated each time macros are edited
function checkCurrentVersion(){
  requires('1.51e');
}



//===========================================================================
macro "Change Values (works on stacks)"{
  for (i=1; i<=nSlices(); i++){  
    setSlice(i);
    changeValues(4,99999,0); }
}




//===========================================================================
//  This macro takes all open stacks and projects them.  Depending on the filename, in this case
//  "C=#", the projection method may be different.  The projection command automatically names
//  the new windows with the type of projection followed by the original filename.
//===========================================================================
macro "Project all open stacks" {
  projectAllOpenStacks();
}  // "Project all open stacks"

function projectAllOpenStacks(){
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
    t = getTitleStripExtension();
    if (  endsWith(getTitle(), "C=2")  )   
      run("Z Project...", "start=1 stop="+nSlices+" projection=[Min Intensity]");
    else 
      run("Z Project...", "start=1 stop="+nSlices+" projection=[Max Intensity]");
    rename(t+"_proj");
  } // for i

}  // end projectAllOpenStacks()


//===============================================================================
// Extracts the nth slice from a timelapse or Z stack.
//===============================================================================
macro "Extract nth slices from stacks"{
 extract_all_nth_timepoints(2);
}

function extract_all_nth_timepoints(extract) {        // extract is an integer between 1 and nSlices
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
    t = getTitleStripExtension();
    setSlice(extract);
    title = "M" + extract + "_" + getTitle;
    run("Duplicate...", "title="+title);
    rename(t+"_s" + extract);
  } // for i
}


macro "=========================="{}

//==============================================================
// This is for when you already have color images.
//macro "box location" {
//  getCursorLoc(x, y, z, mousestate);
//  makeRectangle(x-(boxsize/2),y-(boxsize/2), boxsize, boxsize);
//}

macro "Make snipping box" {
  makeRectangle(10,10, clipsize, clipsize); }

//==============================================================
macro "snip out the box [F5]" {
  roiManager("Add");
  run("Duplicate...", "title=[snip]");
}


//=============================================================================
//  calls function clipCellAndPopOnStack(original, clipstack);
//
macro "clip all ROIs in Roi Manager, center, rotate, and put in stack for montage [F4]" {
  original = getImageID;
  t = getTitle;
  bits = bitDepth();
  if (bits == 24) bitstring="RGB";
  else bitstring = ""+bits+"-bit";
  clipsize = getNumber("Enter side of square box in pixels: ", clipsize);
  setPasteMode("Copy");
  setBatchMode(true);
  newImage("snippedstack "+ t, bitstring+" Black", clipsize, clipsize, 1);
  clipstack = getImageID;
  for (i=0;  i<roiManager("count");  i++){
    selectImage(original);
    roiManager("select", i);
    clipCellAndPopOnStack(original, clipstack);
  }  // for  i
  selectImage(clipstack);
  run("Delete Slice");
  run("Select None");
  setBatchMode(false);
  resetMinAndMax();
}  // macro "clip multiple ROIs in a stack and put in stack for montage"

//=======================================================================
// The ROI of the cell is selected and the original image is active before this function is called.
//  clipstack contains the imageID of the clipped stacks window.
//=======================================================================
function clipCellAndPopOnStack(original, clipstack) {
  requires("1.45m");
  run("Set Scale...", "distance=1 known=1 pixel=1 unit=px");
  run("Set Measurements...", "area centroid perimeter fit shape feret's skewness kurtosis stack display redirect=None decimal=2");
  run("Measure");
  xc = getResult("X", nResults-1);
  yc = getResult("Y", nResults-1);
  selectImage(original);
  makeRectangle(xc-(clipsize/2), yc-(clipsize/2), clipsize, clipsize);
  run("Copy");
  run("Select None");
  setPasteMode("Copy");
  selectImage(clipstack);
  run("Paste");
  a = getResult("Angle", nResults-1) - 90;
  run("Rotate... ", "angle="+a+" grid=1 interpolation=Bilinear ");
  run("Add Slice");
}  // clip cell  




macro "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^"{}

//=============================================================================
//   Sorts cells with bright central spots of tcr from those without.
//
macro "sort cells by ratioing center to edge intensities" {
 // make central ROI & measure
 // make annular ROI & measure
//  depending on ratio, copy to bright cSMAC or lowcSMAC

}

//=============================================================================
//  This macro measures all the openn windows using all the ROIs in the ROI Manager and saves all the results
//  in the same directory as the original images.
//=============================================================================
macro "measure all ROI Manager selections in all open images and save results" {
  for (k=1; k<=nImages; k++) {
    selectImage(k);
    measureAllROIMangagerAndSave();
  }
}


function measureAllROIMangagerAndSave() {
  path = getDirectory("image");
  t = getTitle;
  if (endsWith(t, "tif")) 
    t = replace(t, "tif", "");
  t = t + "txt";
  t = "automeasurements_"+t;
  run("Clear Results");
  run("Set Measurements...", "area mean standard centroid shape redirect=None decimal=3");
  a = newArray(roiManager("count"));
  for (i=0; i<roiManager("count"); i++) a[i]=i;
  roiManager("select", a);
  roiManager("Measure");
  selectWindow("Results");
  IJ.renameResults(t);
  saveAs("text", path+t);
  run("Close");
}  


//=============================================================================
/*
  Framework for iterating through ROI Manager list and spitting out measurments in the Log window.
  In this version, intensity and area are measured.  A band is drawn arounf the outside of the area and this
  mean intensity is treated as the background which is subtracted from the mean of the ROI.
*/
macro "measurements in ROI Manager list"{
  run("Set Measurements...", "area mean standard shape redirect=None decimal=3");
  numberCells = roiManager("count");
  print(getTitle);
  print("Area \t Mean");
  for (i=0; i<numberCells; i++) {
    roiManager("select", i);
    run("Measure");
    run("Enlarge...", "enlarge=2 pixel");
    run("Make Band...", "band=1");   // 1 um band
    run("Measure");
    mean = getResult("Mean", nResults-2) - getResult("Mean", nResults-1); 
    out = "" + getResult("Area", nResults-2) + " \t " + mean;
    print(out);
  } // for i

}

//=============================================================================
/*
    This macro needs to be edited.
*/
macro "measure all ROIs in Roi Manager only including pixel intensities above a threshold" {
  threshold = 220;  // or prompt user for the number -- or start with a background area circled and get this way
//  original = getImageID;
  print(" \t" + getTitle);
  run("Set Measurements...", "area mean redirect=None decimal=5");
  histmax = 16384;
  histmin = threshold;
  for (i=0;  i<roiManager("count");  i++){
//    selectImage(original);
    roiManager("select", i);
    getHistogram(values, counts, histmax-histmin, histmin, histmax-histmin);
    pixelcount = 0;
    integratedint = 0;
    for (h=0; h<histmax-histmin; h++){
      pixelcount = pixelcount + counts[h];
      //integratedint = integratedint + (counts[h] * values[h]);
      integratedint = integratedint + (h * counts[h]);
    }  // for h
  print(i + " \t" + pixelcount + " \t" + integratedint+ " \t" + integratedint/pixelcount);
  }  // for  i
}  // end macro "measure all ROIs in Roi Manager only including pixel intensities above a threshold"


//=========================================================
/*  Since this function was written, a command was added to ImageJ that does this:
    close(pattern)
    Thus, all calls to this function could be replaced.
*/
function closeWindowsNotInterestedIn() {
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
    if ( !(endsWith(getTitle(), "C=2")) )   // use not and the title feature you want to keep
      run("Close");
}   // end function closeWindowsNotInterestedIn()



//=====================================================================================
// Timelapse macros v115   November 2011
//  Macro for ImageJ by Michael Cammer michael@coxcammer.com
//  A roundish (or circular, but with pixel edges, not vectors) ROI is drawn on a timelapse stack such that the 
//  center of the circle is the centroid of the feature to be measured over time and that the entire object at every 
//  time is contained within the perimeter.
//  Then radii spaced by step size perimeter / n step around and are resampled in the time axis
//  and pasted into a new window.
//  After I wrote this, somebody posted a similar plugin version at http://rsbweb.nih.gov/ij/plugins/radial-reslice/
//  This macro produces radial reslices from the center out.  The original version of this reslicing went from the
//  periphery such that the alignment was along the traced outer edge and could be used for the annuluar 
//  intensity and spatial distribution in relation to the outer edge of spreading cells instead of in relation to the center
//  and some slices never even went to the very center.
//=====================================================================================
macro"Radial Sweep as a circle"{
  requires("1.45m");
  run("Set Scale...", "distance=1 known=1 pixel=1 unit=px");
  // MOVED TO GLOBAL VAR  number_of_radii = 2500;                            // number of pixels on perimeter of circle to be sampled.
  run("Set Measurements...", "  centroid redirect=None decimal=0");
  run("Measure");
  xc = getResult("X", nResults-1);
  yc = getResult("Y", nResults-1);
  selectWindow("Results");
  run("Close");
  slices = nSlices();

  step = 2 * PI / number_of_radii;
  extrasteps = 16;  // allows for overshooting to adjust for when slices are projected together

  maxradius = 5 + round(calculateMaximumRadius()/2);  // this is to find the longest radius to make a new window this width or to rescale the shorter lines this length
  original = getImageID();
  originaltitle = getTitle();
  bits = bitDepth();
  newtitle = "reslice_"+originaltitle+"";
  newslices = 1+ number_of_radii + extrasteps;
  newbits = "" + bits + "-bit Black";
  newImage(newtitle, newbits, maxradius, slices, newslices);
  stream = getImageID();

  setBatchMode(true);

  setPasteMode("Copy");
  s=1;
    for (i=0; i<2*PI+(step*extrasteps); i=i+step) {
      edgex = xc + maxradius * sin(i);
      edgey = yc + maxradius * cos(i);
      selectImage(original);
      makeLine(xc, yc, edgex, edgey);
      run("Reslice [/]...", "input=1.000 output=1.000 start=Top");
      run("Select All");
      run("Copy");
      run("Close");
      selectImage(stream);  
      setSlice(s); 
      s++;
      run("Paste");
  }  //for i
 setBatchMode(false);
 resetMinAndMax();
 //projectStackVariableWidth("Max Intensity", 4);
}  //macro("Radial Sweep")

//=========================================================
// This function is required by the macro "Diameter Sweep" to calculate the
// longest possible radius if the selection is non-circular.
// Uses Feret's diameter.
function calculateMaximumRadius (){
  run("Clear Results");
  run("Set Measurements...", "  mean feret's redirect=None decimal=3");
  run("Measure");
  diameter = getResult("Feret",(0));
  selectWindow("Results");
  run("Close");
  return (diameter);
} // function calculateMaximumRadius




//-----------------------------------------------
macro 'Multiply to fix bleaching based on ROI'{
  //checkCurrentVersion();
  original = getImageID();
  run("Set Measurements...", "  mean  redirect=None decimal=0");
  run("Plot Z-axis Profile");    run("Close"); 
  standard = getResult("Mean",(0));
  selectImage(original);
  for (i=1; i<=nSlices; i++){
    bgmean = standard / getResult("Mean",(i-1));
    run("Set Slice...", "slice="+i);
    run("Select All");
    run("Multiply...", "slice value="+bgmean);
  } // for
  run("Select None");
  selectWindow("Results");
  run("Close");
} // macro 'Multiply to fix bleaching based on ROI'


//================================================================
//  See http://microscopynotes.com/imagej/float/ for explanation.
//  originally by Kenton Arkill [underwaterginge@GMAIL.COM]
//    Floating means that each image has the same mean pixel value and the same standard deviation. 
//    I find it really useful as it means the same threshold can often be used on all the images. Here is my macro, 
//    it floats the stack to the values of the first slice (so make sure it is a good image):
//
macro "Float the stack"{
run("32-bit");
setSlice(1);
getRawStatistics(nPixels, meanref, min, max, stdref, histogram);
  for (i = 1; i <= nSlices; i++) {
       setSlice(i);
       getRawStatistics(nPixels, meani, min, max, stdi, histogram);
       run("Subtract...", "value="+meani+" slice");
       run("Divide...", "value="+stdi+" slice");
       run("Multiply...", "value="+stdref+" slice");
       run("Add...", "value="+stdref+" slice");
  }
  run("Enhance Contrast", "saturated=0.35");
  setSlice(1);
}

//===================================================================
//  This macro measures the same region of interest in each slice of a stack and subtracts the mean plus
//  specified standard deviation of the mean from each corresponding slice.  It is useful for situations where
//  there is background intensity changing over time due to something like a flickering light source and you 
//  want the stack's background to be uniform.  It may not be appropriate for situations where precise intensity
//  measurements will be made that should be uniform through the stack.  Subtracting one constant from the
//  entire stack is likely better.
//  To use, put a ROI over a region that is background in all slices.  (Choose the ROI on a maximum pixel 
//  projection if features to assure only bg pixels are included.)  Run the macro.
//
macro 'Subtract background based on ROI'{
  n = 0;  // this is the number of standard deviations to add to the mean bg to subtract
  checkCurrentVersion();
  image_to_process = getImageID();
  run("Set Measurements...", "  mean  standard  redirect=None decimal=5");
  run("Plot Z-axis Profile");    
  run("Close"); 
  selectImage(image_to_process);
  run("Select All");
  for (i=1; i<=nSlices; i++){
    bgmean=getResult("Mean",(i-1)) + n*getResult("StdDev",(i-1));
    run("Set Slice...", "slice="+i);
    run("Subtract...", "slice value="+bgmean);
  } // for loop
  run("Select None");
  run("Min...", "value=0 stack");
  resetMinAndMax();
} // Subtract background based on ROI


//================================================================
//This macro is for drawing a ROI onto a range of serial slices of a stack.
//1.  Draw the ROI.
//2.  Add it to the manager (ctrl-t).
//3.  Choose the color using Edit>Options>Colors...
//4.  Run this macro.
//
macro 'draw ROI  on range of slices in stack'{
  checkCurrentVersion();
  image_to_process = getImageID();  
  Dialog.create("draw on slices...");
  Dialog.addNumber("first slice to draw on", 1);
  Dialog.addNumber("last slice to draw on", nSlices());
  Dialog.show();
  first =  Dialog.getNumber();
  last = Dialog.getNumber();
  if (last < first) { temp=first; first=last; last=temp; }
  if (first < 1) exit("first slice must be 1 or greater");
  if (last > nSlices()) exit("last slice must be the stack size or smaller");
  n = roiManager("count");  // get the last ROI in the manager  
  selectImage(image_to_process);
  roiManager('select', n-1);
  for (i=first; i<=last; i++){
    setSlice(i);
    run("Draw", "slice");
  } // for loop
  run("Select None");
  resetMinAndMax();
} // draw ROI  on range of slices in stack


//=========================================================
macro "RGB merge adjacent slices into a new stack"{
  //checkCurrentVersion();

  setPasteMode("Copy");
  original = getImageID();
  selectImage(original);   
  slices = nSlices()-2; 
  for (i=1;i<=slices;i++) {

    selectImage(original);    
    run("Set Slice...", "slice="+i);  
    run("Select All");  
    run("Duplicate...", "title=[A]");

    selectImage(original);    
    run("Set Slice...", "slice="+(i+1));  
    run("Select All");  
    run("Duplicate...", "title=[B]");

    selectImage(original);    
    run("Set Slice...", "slice="+(i+2));  
    run("Select All");  
    run("Duplicate...", "title=[C]");

    run("RGB Merge...", "red=[A] green=[B] blue=[C]");

    if (i == 1) {
      newrgb = getImageID();
    }
    else {
      run("Select All");
      run("Copy");
      close();
      selectImage(newrgb);   
      run("Add Slice");
      run("Paste");
    }
   
  }  //for i
  selectImage(original); 
  run("Set Slice...", "slice="+nSlices); run("Delete Slice");
  run("Set Slice...", "slice="+1);  run("Delete Slice");
  run("Brightness/Contrast...");

  selectImage(newrgb);  
  run("Select None");

  exit();
} //RGB merge adjacent slices into a new stack


//==================================================================
/*
    This macro is for making a movie of a stack that loops, especially useful for animated GIFs of 3D projections
     that rock back and forth.  Downside:  doubles the size of the stack.
*/
macro "Loop stack" {
    t = getTitle();
    run("Duplicate...", "title=xyztemp duplicate");
    run("Reverse");
    run("Concatenate...", "  title="+t+"_3D image1=["+t+"] image2=[xyztemp] image3=[-- None --]");
}




macro "make test stack for Projections of variable width" {
  size = 25;
  slices = 50;
  width = 256;
  newImage("Untitled", "8-bit black", width, width, slices);
  setColor(255, 255, 255);
  for (i=1; i<slices+1; i++) {
    setSlice(i);
    x = ((i * size) % width);
    y = size * floor ((i*size)/width) + size;
    //print(x,y);
    drawString(""+i ,x , y);
    //makeRectangle(x, y, 10, 10);
    //run("Fill", "slice");
  }
  run("Projections of variable width")
}


//=================================================================
macro "paste with subtract adjacent slices"{
  checkCurrentVersion();
  original = getImageID();
  selectImage(original);    
  setPasteMode("Subtract"); 
  for (i=2;i<=nSlices;i++) {
    run("Set Slice...", "slice="+i);  
    run("Select All");  
    run("Copy");
    run("Set Slice...", "slice="+(i-1));  
    run("Paste");
  }
  selectImage(original); 
  run("Select None");
  run("Set Slice...", "slice="+nSlices); 
  run("Delete Slice");
  run("Set Slice...", "slice="+1);  
  run("Delete Slice");
  run("Brightness/Contrast...");
}  // end  "paste with subtract adjacent slices"


//==================================================================
/*  Version edited 20161002 corrects the global variable projectiondepth getting set to NaN.  It is
    now set correctly to the zdepth entered by the user to be recycled the next time the macro is run.
    A global could also be added for projection type.
*/
macro "Projections of variable width"{
  projtypes = newArray("Max Intensity", "Average Intensity", "Standard Deviation", "Min Intensity", "Median");
  Dialog.create("Projections of variable width");
  Dialog.addMessage("Each slice of the stack is replaced with a projection of the depth\ndefined.");
  Dialog.addChoice("Label", projtypes, "Max Intensity");
  Dialog.addNumber("Slices to project: ", projectiondepth);
  Dialog.show();
  projection = Dialog.getChoice();
  projectiondepth = Dialog.getNumber();
  projectStackVariableWidth(projection, projectiondepth);
}  // Projections of variable width

//==================================================================
/*  Version edited 20161002 corrects the number of slices projected.  Previous versions did zdepth+1 slices.
     Also projects maximum number available.  Previous version stopped one short of the end.
     Also, the slices out to the end are projected at progressively shorter depths.
*/
function projectStackVariableWidth(projtype, zdepth) {
  setBatchMode(true);
  setPasteMode("Copy");
  original = getImageID();
  end = nSlices();
  for (i=1; i<end; i++){
    stop = i + zdepth - 1;
    if (i<(end - zdepth + 2)) 
      run("Z Project...", "start=" + i + " stop=" + stop + " projection=["+projtype+"]");
    else
      run("Z Project...", "start=" + i + " stop=" + end + " projection=["+projtype+"]");
    run("Select All");
    run("Copy");
    run("Close");
    selectImage(original);
    run("Set Slice...", "slice=" + i);
    run("Paste");
  }
  setBatchMode(false);
  run("Select None");
} // end projectStackVariableWidth()

//X================================================X
macro "measure mean  intensity above threshold [q]" {
  // only works on 8 and 16 bit images
  for (s=1; s<=nSlices; s++){
    setSlice(s);
    histogramAboveThreshold(1.01);
  }
}

function histogramAboveThreshold(threshold) {
  getHistogram(values, a, 256);
  counts = 0;
  sum = 0;
  for (i=threshold; i<a.length; i++){
    sum = sum + (i * a[i]);
    counts = counts + a[i];
  }
  print(sum/counts);
}

//===========================================================
// This macro returns stacks resliced from the left by run("Reslice [/]...", "output=1.000 start=Left");
// to their original XYZ orientation.
macro "reorient left resliced stack" {
  run("Reslice [/]...", "output=1.000 start=Top");
  run("Rotate 90 Degrees Right");
  run("Flip Horizontally", "stack");
}

//==============================================
macro "correct depth light absorption" {
  for (i=1; i<=60; i++){
    setSlice(i);
    multby = 1+ ((nSlices-i) / nSlices);
    run("Multiply...", "value="+multby+" slice");
  }
}

//==============================================
macro "multiple range by a constant"{
  for (i=90; i<=nSlices; i++){
    setSlice(i);
    multby = 1.8;
    run("Multiply...", "value="+multby+" slice");
  }
}
//==============================================


//===========================================
macro "Recolor with ratio LUT [F6]"{
  requires ('1.38f');
  makeFavoriteAIFpseudocolorLUT();
}

//===========================================
function makeFavoriteAIFpseudocolorLUT(){
  red = newArray(0,4,8,12,16,20,24,28,32,36,40,44,48,52,56,60,64,68,72,76,80,84,88,92,96,100,104,108,112,116,120,124,128,128,128,128,128,128,128,128,128,128,128,128,128,128,128,128,128,120,112,104,96,88,80,72,64,56,48,40,32,24,16,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,7,11,15,19,23,27,31,35,39,43,47,51,55,59,63,67,71,75,79,83,87,91,95,99,103,107,111,115,119,123,127,131,135,139,143,147,151,155,159,163,167,171,175,179,183,187,191,195,199,203,207,211,215,219,223,227,231,235,239,243,247,251,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255);
  green = newArray(0,1,2,3,4,5,6,7,8,9,10,9,9,8,8,7,7,6,6,5,5,5,4,4,3,3,2,2,1,1,0,0,0,7,15,23,31,39,47,55,63,71,79,87,95,103,111,119,127,135,143,151,159,167,175,183,191,199,207,215,223,231,239,247,255,247,239,231,223,215,207,199,191,183,175,167,159,151,143,135,128,129,131,132,134,135,137,139,140,142,143,145,147,148,150,151,153,154,156,158,159,161,162,164,166,167,169,170,172,174,175,177,178,180,181,183,185,186,188,189,191,193,194,196,197,199,201,202,204,205,207,208,210,212,213,215,216,218,220,221,223,224,226,228,229,231,232,234,235,237,239,240,242,243,245,247,248,250,251,253,255,251,247,243,239,235,231,227,223,219,215,211,207,203,199,195,191,187,183,179,175,171,167,163,159,155,151,147,143,139,135,131,127,123,119,115,111,107,103,99,95,91,87,83,79,75,71,67,63,59,55,51,47,43,39,35,31,27,23,19,15,11,7,3,0,8,16,24,32,41,49,57,65,74,82,90,98,106,115,123,131,139,148,156,164,172,180,189,197,205,213,222,230,238,246,254);
  blue = newArray(0,7,15,23,31,39,47,55,63,71,79,87,95,103,111,119,127,135,143,151,159,167,175,183,191,199,207,215,223,231,239,247,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,251,247,243,239,235,231,227,223,219,215,211,207,203,199,195,191,187,183,179,175,171,167,163,159,155,151,147,143,139,135,131,128,126,124,122,120,118,116,114,112,110,108,106,104,102,100,98,96,94,92,90,88,86,84,82,80,78,76,74,72,70,68,66,64,62,60,58,56,54,52,50,48,46,44,42,40,38,36,34,32,30,28,26,24,22,20,18,16,14,12,10,8,6,4,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,16,24,32,41,49,57,65,74,82,90,98,106,115,123,131,139,148,156,164,172,180,189,197,205,213,222,230,238,246,254);
  setLut(red, green, blue);
}

//============================================================================================
//  LOCI open to 5D
//  Adjust colors
//  Run this for RGB results

macro "convert WYS Image5D to RGB stack" {
  convertImage5D_to_RGB_stack()
}  // macro


//==============================================================
//  This function converts 5D stacks as they appear ("what you see") to RGB.
//  Resets the Z size of voxels to equal XY so that kymographs are scaled 1:1.
//  
function convertImage5D_to_RGB_stack() {

  setBatchMode(true);
  original5D = getImageID();
  getDimensions(nw,nh,nch,nsl,nfr);
  getVoxelSize(width, height, depth, unit);
  depth = width;                                       // Resets the Z size of voxels to equal XY so that kymographs are scaled 1:1.
  timeinterval = Stack.getFrameInterval();                     // Returns the frame interval in T units
  Stack.getUnits(X,Y,Z,T,Value);             // Returns the x, y, z, t and value units


  // make new RGB image for result
  title = getTitle() + "RGB";
  newImage(title, "RGB Ramp", nw, nh, nfr);
  newRGB = getImageID();
  setPasteMode("Copy");

  for(i=1; i<=nfr; i++) {
    selectImage(original5D);
    run("Set Position", "x-position=1 y-position=1 channel=2 slice=1 frame="+i+" display=overlay");
    run("Image5D Stack to RGB");
    run("Select All");
    run("Copy");
    close();
    selectImage(newRGB);
    setSlice(i);
    run("Paste");

  } // for i

  setBatchMode(false);
  selectImage(newRGB);
  setVoxelSize(width, height, depth, unit); 
  Stack.setFrameInterval(timeinterval);         // Sets the frame interval in T units
  Stack.setTUnit(T);                          // Sets the T (time) unit 

} // convertImage5D_to_RGB_stack() 


//=====================================================================================
//  This macro is not ready for general use because it requires programming to work right.
//  In its current state it moves the field a fixed amount per step and this amount is 1/22.
macro "Drift Compensation" {
  setPasteMode("Copy");
  for (i=0; i<nSlices; i++) {
    setSlice(i+1);
    makeRectangle(30, 24, 463, 466);
    run("Copy");
    makeRectangle(30, 24-(i/22), 463, 466);
    run("Paste");
  }  // for i
}  //  Drift Compensation


//=======================================================================================
// 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" {
  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

// Make noisy image for Pearson's test
//********************************************************************
//macro 'make noisy image for test'{
//  run("New...", "name='noisy image' type=32-bit fill=Ramp width=64 height=64 slices=2");
//  for (z=1;z<=nSlices();z++) {run("Set Slice...", "slice="+z);
//    for (x=0;x<getWidth();x++) for (y=0;y<getHeight();y++) 
//      putPixel(x,y,random());
//  } // z
//}

//macro "test the array"{
//  testa = newArray(3);  testb = newArray(3);
//  testa[0]=1;  testa[1]=2;  testa[2]=3;  testb[0]=2;  testb[1]=5;  testb[2]=6;
//  //testa[0]=1;  testa[1]=1;  testa[2]=1;  testb[0]=1;  testb[1]=1;  testb[2]=1;
//  coef = calculateCoefficient(testa, testb);  
//}  // END OF PEARSON'S MACROS

//===================================================================
//  Applies CLAHE to a stack of images.
//===================================================================
macro "CLAHE on a stack" {
  original = getImageID;
  Dialog.create("CLAHE on a stack");
  blocksize=127;   histogramsize=256;   maximum=3;
  Dialog.addNumber("Block size", blocksize);
  Dialog.addNumber("Histogram bins", histogramsize);
  Dialog.addNumber("Contrast max (3 to 12)", maximum);  // from 3 to 12
  Dialog.show();
  blocksize=Dialog.getNumber(); histogramsize=Dialog.getNumber(); maximum=Dialog.getNumber(); 
  for (i=1; i<=nSlices; i++) {
    selectImage(original);
    setSlice(i);
    run("CLAHE ", "blocksize="+blocksize+" histogram="+histogramsize+" maximum="+maximum);
  } // for i
  resetMinAndMax();
}  // CLAHE on a stack



//===================================
macro "Save ROI Manager [F7]" {
  if ( getInfo("window.type") != "Image" )  exit("Make sure running macro from image window.");
  dir = getDirectory("image");                         // where to save results
  if (lengthOf(dir) < 1) exit("Make sure the image was opened from disk.");
  title = getTitleStripExtension();    
  newName =dir + title + ".zip";
  // add command to select none in ROIManager to mae sure all saved?
  roiManager("Save", newName);
  roiManager("reset");
}  // end Save ROI Manager [F7]



/====================================
// Use this function to strip any number of extensions // off images.
// Returns the title without the extension.  
//====================================
function getTitleStripExtension() {
  t = getTitle();
  extensions = newArray(".tif", ".tiff", ".lif", ".lsm", ".czi", ".nd2", ".ND2");    
  for(i=0; i<extensions.length; i++)
    t = replace(t, extensions[i], "");  
  return t;
}  // end getTitleStripExtension()/*  Macros for Ashira Blazer  */