// ImageJ macro to calculate the Radially Averaged Autocorrelation Function,
// Corrected for finite size effects
// The output is normalized to a value of 1 at zero radius
// and corrected for effects of the finite image size.
// (No need to extend the image size for avoiding edge effects)
//
// Use with a binary input image.
// Needs the "Radial Profile" plugin, http://rsb.info.nih.gov/ij/plugins/radial-profile.html
//
// Version: 2008-May-14 Michael Schmid
//
requires("1.42p");
doStack=false;
if (nSlices()>1) doStack = getBoolean("Get RDF from full stack?");
setBatchMode(true);
firstSlice=getSliceNumber();
lastSlice=getSliceNumber();
if (doStack) {
  firstSlice=1;
  lastSlice=nSlices();
}
width=getWidth;
height=getHeight;
getPixelSize(unit, pixelWidth, pixelHeight);
run("Select None");
//maxRadius may be modified, should not be larger than 0.3*minOf(width, height);
maxRadius=0.3*minOf(width, height);
minFFTsize=1.3*maxOf(width, height);
title=getTitle();
size=4;
while(size<minFFTsize) size*=2;
for (slice=firstSlice; slice<=lastSlice; slice++) {
  //make autocorrelation of particle image
  tempTitle="temp-"+random();
  run("Duplicate...", "title="+tempTitle);
  tempID=getImageID();
  getRawStatistics(nPixels, mean);
  run("Canvas Size...", "width="+ size+" height="+ size+" position=Center zero");
  makeRectangle(floor((size-width)/2), floor((size-height)/2), width, height);
  run("Make Inverse");
  run("Set...", "value="+mean);
  run("Select None");
  getRawStatistics(nPixels, mean);
  run("FD Math...", "image1=["+tempTitle+"] operation=Correlate image2=["+tempTitle+"] result=AutoCorrelation do");
  psID=getImageID();
  run("Subtract...", "value="+(nPixels*mean*mean));
  selectImage(tempID);
  close();

  //make autocorrelation reference to correct finite image size effects
  newImage("frame", "8-bit White", width, height, 1);
  run("Set...", "value=255");
  tempID=getImageID();
  rename(tempTitle);
  run("Canvas Size...", "width="+ size+" height="+ size+" position=Center zero");
  run("FD Math...", "image1=["+tempTitle+"] operation=Correlate image2=["+tempTitle+"] result=AutoCorrReference do");
  refID=getImageID();
  imageCalculator("Divide", psID,refID);
  selectImage(refID);
  close();
  selectImage(tempID);
  close();

  //prepare normalized power spectrum for radial averaging
  selectImage(psID);
  circleSize=2*floor(maxRadius)+1;
  norm = getPixel(size/2, size/2);
  run("Divide...", "value="+norm);
  run("Specify...", "width="+circleSize+" height="+circleSize+" x="+(size/2+0.5)+" y="+(size/2+0.5)+" oval centered");
  run("Radial Profile", "x="+(size/2+0.5)+" y="+(size/2+0.5)+" radius="+floor(maxRadius)-1);
  Plot.getValues(x, y);
  plotID=getImageID();
  selectImage(psID);
  close();
   if (slice==firstSlice) ySum = newArray(y.length);
  for (i=0; i<y.length; i++)
    ySum[i]+=y[i]/nSlices;
  selectImage(plotID);
  close();
}
if (pixelWidth==pixelHeight) {
  for (i=0; i<x.length; i++)
    x[i] *= pixelWidth;
} else
  unit = "pixels";
if (doStack) title = title + " (stack)";
Plot.create("Autocorrelation of "+title, "Distance ("+unit+")", "Normalized Autocorrelation", x, ySum);

setBatchMode(false);