// 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 clipsize = 150;  // edge size of the box for clipping cells for montage
var snipwindow = -1;
var projectiondepth = 15;

function checkCurrentVersion(){
  requires('1.46a');
}

//==========================================================================
// User needs to change path and name is the last line of code.  This could be made more user friendly by
// prompting user to choose directory and click the "Select" button.
// User should change this to delete # of characters from beginning of filename.  It is used here to get rid of
// redundant filenames.  Set to zero or remark out to do nothing.
//==========================================================================
macro "Rename Open WIndows and Save All" {
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
    run("Canvas Size...", "width=600 height=600 position=Center zero");
    // save
   // close()
  }
}


//===============================================================================
macro "pad windows and save" {
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
    a = getTitle();
    a = replace(a, ".nd2", "");    // delete the .nd2 extension
    // The next line of code has to be edited manually for each file!
    a = substring(a, 6);    // <===============
    rename(a);
    saveAs("Tiff", "D:\\Catarina Sacristán\\20111214 LCMV Cl13 d7 P14 CD8 GFP 12_14_11 CSR\\"+a+".tif");
  }


//===============================================================================
// Use:  Open .nd2 file which is a multi-channel multi-field time lapse.  This will select from these fields either
// projections or single timepoints.
//  The purpose of this macro is to extract all the nth 
//===============================================================================
macro "Extract or otherwise process all open .nd2 stacks [q]" {
   // function convertProjectedImagesToStacks();    // use this or extract as first command
   extract_all_nth_timepoints(2);
   close_all_nonM_windows();
   convertProjectedImagesToStacks();

   // additional processing!!!
     n = nImages();
     for (i=1; i<=n; i++) {
       selectImage(i);
       run("Bandpass Filter...", "filter_large=160 filter_small=0 suppress=None tolerance=5 process");
     }

}


//===========================================================================
//  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.
//===========================================================================
function projectAllOpenStacks(){
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
    // run("Bandpass Filter...", "filter_large=200 filter_small=1 suppress=None tolerance=0");
    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]");
  } // for i
}


//===============================================================================
// Extracts from a timelapse or Z stack the nth slice only.
//===============================================================================
function extract_all_nth_timepoints(extract) {        // extract is an integer between 1 and nSlices
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
    setSlice(extract);
    title = "M" + extract + "_" + getTitle;
    run("Duplicate...", "title="+title);
    rename(title);
  } // for i
}

macro "Extract All First Timepoints"{
 extract_all_nth_timepoints(1);
}

macro "Extract All 2nd Timepoints"{
 extract_all_nth_timepoints(2);
}

macro "Extract All 5th Timepoints"{
 extract_all_nth_timepoints(5);
}

//===============================================================================
// Extracted or projected files begin with the letter "M".  As long as the original files do not have titles that
// begin with the letter "M", all original stacks are removed and only the extracted files remain.
//===============================================================================
function close_all_nonM_windows() {
  for (i=nImages(); i>0; i--) {
    selectImage(i);
    if (!(  startsWith(getTitle(), "M")  ) )
      close();
  } // for i
} // end function close_all_non-M_windows(); 


function convertProjectedImagesToStacks() {
  run("Images to Stack", "name=Stack title=C=0 use");
  rename("GFP");
  run("Images to Stack", "name=Stack title=C=1 use");
  rename("PDL1");
  run("Images to Stack", "name=Stack title=C=2 use");
  rename("reflection");
  run("Images to Stack", "name=Stack title=C=3 use");
  rename("ICAM");
} // end function convertProjectedImagesToStacks();


//============================================================================
// This macro was written for Hila who has a series of nd2 files of fixed stained material.
//  Each file consists of n channels.  They are all imported with channels split and this macro
//  copies them into stacks.
macro "sort open files of split channels into stacks using name C=#" {
  t = replace(getTitle(), " ", "_");
  t = substring(t, 0, (indexOf(t, "nd2")));
  t = replace(t, "\\.", "");
  run("Images to Stack", "name="+t+"_C=0 title=C=0 use");
  run("Images to Stack", "name="+t+"_C=1 title=C=1 use");
  run("Images to Stack", "name="+t+"_C=2 title=C=2 use");
  run("Images to Stack", "name="+t+"_C=3 title=C=3 use");
  run("Images to Stack", "name="+t+"_C=4 title=C=4 use");
  run("Images to Stack", "name="+t+"_C=5 title=C=5 use");
  run("Images to Stack", "name="+t+"_C=6 title=C=6 use");
  run("Images to Stack", "name="+t+"_C=7 title=C=7 use");
  run("Images to Stack", "name="+t+"_C=8 title=C=8 use");
  run("Images to Stack", "name="+t+"_C=9 title=C=9 use");
  run("Images to Stack", "name="+t+"_C=10 title=C=10 use");
}
//========================

//============================================================================
// 
macro "RENAME C=#" {
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
    t = replace(getTitle, " ", "_");
    t = replace(t, "\\.nd2", "");
    t = replace(t, "\\.tif", "");
    t = replace(t, "C=0", "ICAM red");
    t = replace(t, "C=1", "gp120 green");
    t = replace(t, "C=2", "anti-CD3 blue");
    t = replace(t, "C=3", "reflection");
    t = replace(t, "C=4", "trans");
    rename(t);
  } // i
}
//========================
//============================================================================
// 
macro "MAKE ALL RGB" {
  n = nImages();
  for (i=1; i<=n; i++) {
    selectImage(i);
   run("RGB Color");
  } // 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]");
}


//=============================================================================
//  Modified 20121001 for composite color hyperstacks.
//  calls function clipCellAndPopOnStack(original, clipstack);
//
macro "clip all ROIs in Roi Manager, center, and put in stack for montage [F4]" {
  if (! Stack.isHyperstack) exit('Expecting a hyperstack');
  original = getImageID;
  Stack.getDimensions(width, height, channels, slices, frames);
  Stack.getDisplayMode(displaymode);
  Stack.getActiveChannels(activechannel);
  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");
  clipstackID=newArray(channels);
  clipstackTitle=newArray(channels);
  for (c=1; c<=channels; c++){
    setBatchMode(true);
    newImage("snippedstack "+ t + " C=" + c, bitstring+" Black", clipsize, clipsize, 1);
    clipstackTitle[c-1] = getTitle();
    clipstackID[c-1] = getImageID;
    for (i=0;  i<roiManager("count");  i++){
      selectImage(original);
      roiManager("select", i);
      clipCellAndPopOnStack(original, clipstackID[c-1], c);
    }  // for  i
    selectImage(clipstackID[c-1]);
    run("Delete Slice");
    run("Select None");
    setBatchMode(false);
    resetMinAndMax();
  } //for c channels

 // This block of code converts the individual snipped stacks back to a composite color stack.
  frames = nSlices();
  catstring = " title=[tempcat] ";
  for(i=0;i<clipstackTitle.length; i++){
    catstring = catstring + "image"+(i+1)+"=[" + clipstackTitle[i] + "] ";
  }  // done constructing a string
  catstring = catstring +  "image" + (i+1) +"=[-- None --]";
  print(catstring);
  run("Concatenate...", catstring);
  run("Stack to Hyperstack...", "order=xytcz channels="+channels+" slices=1 frames="+frames+" display=Color");
  rename("snippedstack "+ t );
  Stack.setDisplayMode(displaymode);
  Stack.setActiveChannels(activechannel);

   // This block of code transfers the LUTs and scaling from the original images to the snippedstack
   resultImage = getImageID();
   for (c=1; c<=channels; c++){
      selectImage(original);
      Stack.setChannel(c);
      getLut(reds, greens, blues);
      getMinAndMax(min, max);
      selectImage(resultImage);
      Stack.setChannel(c);
      setLut(reds, greens, blues);
      setMinAndMax(min, max);
   } // for c  to update LUTs

}  // macro "clip multiple ROIs in a stack and put in stack for montage"

//=======================================================================
// Modified 20121001 for composite color hyperstacks.
// 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, chan) {
  requires("1.45m");
  run("Set Scale...", "distance=1 known=1 pixel=1 unit=px");
  run("Set Measurements...", "  centroid area redirect=None decimal=0");
  run("Measure");
  xc = getResult("X", nResults-1);
  yc = getResult("Y", nResults-1);
  selectWindow("Results");
  run("Close");
  selectImage(original);
  Stack.setChannel(chan);
  makeRectangle(xc-(clipsize/2), yc-(clipsize/2), clipsize, clipsize);
  run("Copy");
  run("Select None");
  setPasteMode("Copy");
  selectImage(clipstack);
  run("Paste");
  run("Add Slice");
}  // clip cell  

//================================================================================================
macro "rotate image multiple times and average" {
// get minimum value
// rotate any angle
// at end replace all zeroes with minimum value
  setBatchMode(true);
  rotation = 1;
  original = getImageID;
  run("Select None");
  for (i=0; i<(360/rotation)-1; i++){
    setSlice(1);
    run("Duplicate...", "title=[]");
    run("Rotate... ", "angle="+(rotation*(i+1))+" grid=1 interpolation=Bilinear");
    run("Select All");
    run("Copy");
    close();
    selectImage(original);
    run("Add Slice");
    run("Paste");
    deletemestack = getImageID();
  }  // for i
  setBatchMode(false);
  run("Z Project...", "start=1 stop=360 projection=[Average Intensity]");
  resultimage = getImageID();
  selectImage(deletemestack); 
  close();
  selectImage(resultimage);
  //makeLine(0, 0, getWidth(), getHeight(), 1);
  makeLine(getWidth()/2,0,getWidth()/2,getHeight(),1);
  //run("Plot Profile");
  p = getProfile();
  min = 999999999;
  for(i=0; i<p.length/2;i++)
    if (p[i] < min) min=p[i];
  for(i=0; i<p.length;i++)
    //p[i] = p[i] - min;
    print(p[i]-min);
}

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


//=============================================================================
//  This macro requires two files.  One file is a hyperstack of images where there are multiple channels (2 to n) and
//  many frames.  The second file is ROIs which are opened into the ROI Manager.
//  The macro selects each ROI from the manager and measures it in each channel.  It makes a second pass to measure an annulus
//  of the ROI immediately exterior to the ROI.  
//  At the end, the results are exported to a tab delimitted table.
macro "Simple measurements of multiple channels in hyperstack"{
  if (! Stack.isHyperstack) exit('Expecting a hyperstack');
  original = getImageID;
  Stack.getDimensions(width, height, channels, slices, frames);
  Stack.getActiveChannels(activechannel);
  t = getTitle;
  numberCells = roiManager("count");
  if (numberCells < 1) exit("Needs ROI Manager with list of ROIs.");

  run("Clear Results");
  run("Set Measurements...", "area mean standard shape integrated redirect=None decimal=2");
  for (m=0; m<numberCells; m++) {
    roiManager("select", m);
    for (c=1; c<=channels; c++){
      Stack.setChannel(c);
      Stack.setActiveChannels(c);
      run("Measure");
    } // for c
  } // for m measurements

  //Now remeasure the background intensities
  for (m=0; m<numberCells; m++) {
    roiManager("select", m);
    run("Enlarge...", "enlarge=2 pixel");
    run("Make Band...", "band=1");   // 1 um band
    for (c=1; c<=channels; c++){
      Stack.setChannel(c);
      Stack.setActiveChannels(c);
      run("Measure");
    } // for c
  } // for m measurements

  // report the results
  print(t);
  // print the column headers
  line = "cell# \t";
  for (c=1; c<=channels; c++) { 
    line = line +"ch"+c+"_rawmean \t"+ "ch"+c+"_mean \t"+ "ch"+c+"_area \t"+ "ch"+c+"_stdev \t"+ "ch"+c+"_circ \t"+ "ch"+c+"_intden \t";
  } // for c  printing column headers
  print(line);
  //Output the results
  bgOffset = numberCells * channels;
  for (i=0; i<bgOffset; i=i+channels) {
    line = ""+(1+i/channels)+" \t ";
    for (c=0; c<channels; c++) {  // the channel offset
      rawmean = getResult("Mean", i+c);
      mean = getResult("Mean", i+c) - getResult("Mean", bgOffset+i+c); 
      area = getResult("Area", i+c);
      stdev = getResult("StdDev", i+c);
      circ = getResult("Circ.", i+c);
      intden = getResult("IntDen", i+c);
      line = line + rawmean+" \t "+mean+" \t "+area+" \t "+stdev+" \t "+circ+" \t "+intden+" \t " ;
    } // for c
    print(line);
  } // for i  reporting results
} // end "Simple measurements of multiple channels"

//=====================================================================
//  This macro looks for regions of interest in the ROI Manager and measures each one on an
//  image.  Because ROIs in the manager are stored with the slice number, this macro will operate on a stack
//  of images.  The ROIs have to be closed polygons.  As of this version, the macro does not error check
//  to make sure each ROI is the correct type.  A future version should check the ROI type and skip it if it is not
//  a closed shape,
//  The function radialSweepCircle(ROInumber) is called to generate a sweeping radial kymograph.  This
//  macro could easily be modified to make other radial measurments such as line profile plots such as center-out
//  version of the perimeter-in measurement method generated for the Condeelis lab in 2004-2005 as posted at
//  http://cammer.net/aecom/radialsweep/01/  user: aecom  pwd: newpassowrd
//=============================================================================
macro "Radial Sweep multiple ROIs in a stack for intensity profiles" {
   if (! Stack.isHyperstack) exit('Expecting a hyperstack');
   original = getImageID;
   Stack.getDimensions(width, height, channels, slices, frames);
   Stack.getActiveChannels(activechannel);
   t = getTitle;
   numberCells = roiManager("count");
   if (numberCells < 1) exit("Needs ROI Manager with list of ROIs.");
   resultstring = newArray(256);

  for (channel=1; channel<=channels; channel++){
    dataIndex = 0;
    //setBatchMode(true);
    newImage("data", "32-bit Black", 256, 512, 1);   // This is the data window where the arrays of line profiles are stored
    datawindow = getImageID;
    for (i=0;  i<numberCells;  i++){
      selectImage(original);
      roiManager("select", i);
      Stack.setChannel(channel);
      radialSweepCircleLineProfiles(i);
    }  // for  i

    selectImage(datawindow);
    makeRectangle(0, 0, 256, 79, numberCells);
    run("Crop");

    run("Select All");
    //run("Plot Profile");
    run("Clear Results");
    profile = getProfile();
    for (i=0; i<profile.length; i++)
      resultstring[i] = "" + resultstring[i] + profile[i] +"\t";
    selectImage(datawindow);
    close();
  } // for channel
  //setBatchMode(false);
  print(t);
  for (i=0; i<profile.length; i++)
    print(resultstring[i]);

}  // macro "Radial Sweep multiple ROIs in a stack for intensity profiles"

//=======================================================================
function radialSweepCircleLineProfiles(ROInumber) {
  run("Set Scale...", "distance=1 known=1 pixel=1 unit=px");
  number_of_radii = 360;                            // number of pixels on perimeter of circle to be sampled.
  run("Set Measurements...", "  centroid area redirect=None decimal=0");
  run("Measure");
  xc = getResult("X", nResults-1);
  yc = getResult("Y", nResults-1);
  area = getResult("Area", nResults-1);
  selectWindow("Results");
  run("Close");
  step = 2 * PI / number_of_radii;
  // make a new window; could be broken out into a function
  maxradius = 20 + round(calculateMaximumRadius()/2);  // this is to find the longest radius to make a new window; a few addl pixels are added to length for bg
  original = getImageID();
  originaltitle = getTitle();
  bits = bitDepth();
  newtitle = "" + ROInumber+1 + " radialprofile_"+area+" "+originaltitle+"";
  newY = number_of_radii;
  newbits = "" + bits + "-bit Black";
  newImage(newtitle, newbits, maxradius, newY, 1);
  profilewindow = getImageID();
  setBatchMode(true);               // speeds up the radial resampling
  setPasteMode("Copy");
  Y=0;
    for (i=0; i<2*PI+(step); i=i+step) {
      edgex = xc + maxradius * sin(i);
      edgey = yc + maxradius * cos(i);
      selectImage(original);
      makeLine(xc, yc, edgex, edgey);
      profile = getProfile();
      selectImage(profilewindow);  
      for(X=0; X<profile.length; X++){
        putPixel(X,Y,profile[X]);
      }  // for X
      Y++;
    }  //for i
 setBatchMode(false);              //  end of radial resampling
 // report the profile from the center out
 selectImage(profilewindow); 
 run("Select All");
 resetMinAndMax();
 profile = getProfile();
 close();
 selectImage(datawindow);
 for(i=0; i<profile.length; i++){
   putPixel(i,dataIndex,profile[i]);
 }
 dataIndex++;
}  //macro Radial sweep for line profiles

//========================================================
macro "reslice line ROIs from ROIManager"{
  original = getImageID;
  t = getTitle;
  numberCells = roiManager("count");
  if (numberCells < 1) exit("Needs ROI Manager with list of ROIs.");
  run("Clear Results");
  for (m=0; m<numberCells; m++) {
    selectImage(original);
    roiManager("Select", m);
    run("Reslice [/]...", "output=0.000 start=Top avoid");
    run("Properties...", "channels=1 slices=1 frames=1 unit=µm pixel_width=0.3321064 pixel_height=5 voxel_depth=0.3321064 frame=[0 sec] origin=0,0");
    rename("line#"+m+" "+t);
  } // for m measurements
} // reslice line ROIs from ROIManager


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



//=============================================================================
//   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 open 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");
}  


//=============================================================================
//  
//
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 + (h * counts[h]);
    }  // for h
  print(i + " \t" + pixelcount + " \t" + integratedint+ " \t" + integratedint/pixelcount);
  }  // for  i
}  // "measure all ROIs in Roi Manager only including pixel intensities above a threshold"




//===================================================
// Criteria for counting the cells are:
// cell area greater than 10 um^2
macro "Count cells by center brighter than edges" {
  countSmallCells = 0;
  getPixelSize(unit, pixelWidth, pixelHeight);
  roiManager("deselect");
  run("Set Measurements...", "area mean standard centroid redirect=None decimal=3");
  print("Reporting measurements for");
  print(getDirectory("image") +" "+getTitle);
  original = getImageID;
  print(roiManager("count")+"  total number of cells counted by reflection, GFP, or user.");

  ICAMrois = newArray(roiManager("count"));
  for (i=0; i<roiManager("count"); i++) ICAMrois[i]=-1;
  i=0;
  while (i<roiManager("count")) {
    roiManager("select", i);
    run("Measure");
    xc = getResult("X", nResults-1);   yc = getResult("Y", nResults-1);
    selectionArea = getResult("Area", nResults-1); 
    cellMean=getResult("Mean", nResults-1); 
    if (selectionArea > 10) {  // if cell footprint is bigger than 10 um^2
      run("Enlarge...", "enlarge=1 pixel");
      run("Make Band...", "band=1");
      run("Measure");
      bandMean=getResult("Mean", nResults-1); 
      //print(cellMean / bandMean);
      if (cellMean / bandMean > 1.10) {  // This is the amount over bg the cell needs to be to be considered positive
        ICAMrois[i] = i;    //run("Draw");
      } // if cellmean
    } // if selectionarea
    else {
      countSmallCells++;
    } // else
    i++;
  } // while
  print(countSmallCells+" cells too small (less than 10^2).");
  //print(roiCount+" cells considered ICAM positive.");
  positiveCount=0;
  for (i=0; i<roiManager("count"); i++) {
    if ( ICAMrois[i] > 0) {
      roiManager("select", i);
      positiveCount++;
    } // if
  }  // for
  print (positiveCount+"  cells considered positive.");
  print ("");
  selectWindow("Log");
}  // end "Count cells by center brighter than edges"


//=========================================================
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()

//=========================================================
macro "Process movies for timelapse analysis"{
  open the full .nd2 file

  //eliminate all stacks except the channels used for cell extraction
  closeWindowsNotInterestedIn();

for each reflection stack:
  cover up the black edges
  //fft 255-0 twice
  run("Bandpass Filter...", "filter_large=200 filter_small=0 suppress=None tolerance=0");
  maybe smooth?
  linear contrast if necessary
  threshold huang
  make binary
  save for Imaris or later ImageJ

In Imaris, have tracks made then manually eliminate touching cells 
This involves opening each movie.

}  // end Process movies for timelapse analysis





//=====================================================================================
// 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");
  number_of_radii = 360;                            // 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'


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


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


//==================================================================
macro "Remove Outliers and AVG Projection of variable width"{
  checkCurrentVersion();

  projectiondepth = 5;
  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, "Average Intensity");
  Dialog.addNumber("Slices to project: ", projectiondepth);
  Dialog.show();
  projection = Dialog.getChoice();
  projectiondepth = projection;
  zdepth = Dialog.getNumber();
  run("Remove Outliers...", "radius=1 threshold=0 which=Bright stack");
  run("Remove Outliers...", "radius=1 threshold=0 which=Bright stack");
  projectStackVariableWidth(projection, zdepth);
}  // Projections of variable width




//==================================================================
macro "Projections of variable width"{
  checkCurrentVersion();
  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 = projection;
  zdepth = Dialog.getNumber();
  projectStackVariableWidth(projection, zdepth);
}  // Projections of variable width

//==================================================================
function projectStackVariableWidth(projtype, zdepth) {
  setBatchMode(true);
  setPasteMode("Copy");
  original = getImageID();
  end = nSlices();
  for (i=1;i<(end-zdepth+2);i++){
    stop = i + zdepth;
    run("Z Project...", "start="+i+" stop="+stop+" 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 [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

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