/*  This macro goes through regions of interest drawn in XY on a stack of 8 or 16 bit grayscale images and reports the mean
intensity of the top n pixels per region per slice.
The purpose is to look at change of intensity through time and allows for regions to be drawn larger than the feature of interest
so that some movement of  the sample is acceptable.  As long as two features of interest do not intersect during the duration of the
sequence, this method works well.
Variable "topn" contains the number of brightest pixels to be averaged and reported.  topn=1 is the equivalent as measuring the maximum 
intensity pixel.

Instructions:
Image must be a file already saved to disk because the data are output to a text file and needs the path.
To decide where to draw the regions, do a maximum pixel projection and darw ROIs that completely contain each feature.
Add each region to the ROI Manager.
Save the ROI manager using the macro "Save ROI Manager [F7]".
Select the original movie and run macro "Measure top n pixels mean intensity [q]".

Before running, you may want to do a median 1x1 filter onthe stack to reduce noise.


For each ROI, the intensities at each slice are measured.
There is an outer loop, for each roi, and an inner loop, for each slice.
Within these loops, a measurement method may be inserted.  

The first column reports the frame number.
If the time info were in the metadata of the file, then instead of reporting the frame number, the time could be reported.

*/



var topn = 20;
var tab = " \t ";

macro "Measure top n pixels mean intensity [q]" {
  run("Set Measurements...", "mean min stack limit display redirect=None decimal=2");
  path = getDirectory("image");
  title = getTitleStripExtension();
  s = nSlices();
  output = newArray(s);
  for (slice=0;  slice<s; slice++) 
    output[slice] = "" + (slice+1) + tab;
  for (roi=0; roi<roiManager("count"); roi++)  {
      roiManager("select", roi);
      for(slice=0; slice<s; slice++) {   // for each slice
        setSlice(slice+1);
          // begining of method to get value to report -- the single value is in variable histMax
          //run("Measure");
          //histMax = getResult("Max", nResults-1);
          histMax = returnAverage_top_n_pixels(topn);
          // end of method to get value to report -- the single value is in variable histMax        
         output[slice]= output[slice] + histMax + tab;        // there is a way to output data to a table window instead of this method
       } // for k
  }  // for each ROI

  topline = "timepoint\t";
  for (roi=0; roi<roiManager("count"); roi++)  {
    topline = topline + call("ij.plugin.frame.RoiManager.getName", roi) + tab;
  }

  outputFile = File.open(path+title+"_results.txt");
  print(outputFile, topline);
  for (slice=0; slice<output.length; slice++) {
    print (outputFile, output[slice]);
  }
  File.close(outputFile);
  // save to file in same directory 

}



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


//======================================
macro "Mean of top n pixels " {
  print (returnAverage_top_n_pixels(10));

}


function returnAverage_top_n_pixels(n) {
  if (bitDepth() > 8) 
    bins = 65536;
  else
     bins = 256;
  getHistogram(dummy0, counts, bins);        // use counts 
  numberCountedSoFar = 0;
//  currentPosition = counts.length-1;    // for full 16 bits
  currentPosition = 16384;  // used to speed up when max cannot fill all 16 bits
  sum = 0;
  do  {
     if ((numberCountedSoFar + counts[currentPosition]) > n) 
       x = n - numberCountedSoFar;
     else 
       x = counts[currentPosition];
     numberCountedSoFar = numberCountedSoFar + x;
     sum = sum + (x * currentPosition);
     currentPosition--;
  } while (numberCountedSoFar < n);
  mean = sum / n;
  return mean;
}



macro "test image" {
//  intensities = newArray(
//    1000, 1000, 1000, 1000, 1000, 1000, 2000, 2100, 2200, 2300, 900);
  intensities = newArray(
    1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000);
  Array.getStatistics(intensities, min, max, mean, stdDev);
  print(min, max, mean);
  newImage("image", "16-bit black", 256, 256, 1);
  cornerx = 10;
  cornery = 10;
  for (x=0; x<intensities.length; x++) {
    putPixel(x*3+cornerx, cornery, intensities[x]);
  }
  resetMinAndMax();
}