Radial Distribution Function

An ImageJ macro that calculates the radial distribution function (RDF) of particle centers. Needs the Radial_Profile plugin.

The macro is based on the calculation of the autocorrelation via the Fourier transform, but it does not assume periodic boundary conditions. It increases the image size to a value large enough to avoid any influence of the 'periodic' nature of the FFT and corrects for finite-size effects (edge effects) that would occur because particles close to the edge necessarily have fewer distant neighbors than particles near the center (this correction is exact in case of very small particles only; it has limited accuracy for large particles).

For the autocorrelation of the image (all pixels, not only the centers) see the Radially Averaged Autocorrelation macro.

Performance: For a 2048*2048 image (expanded to 4096*4096 internally), computing time is about half a minute on a 2.4 GHz Core2Duo; memory requirement about 400 MB.

// ImageJ macro to calculate the Radial Distribution Function (RDF) of particle centers
//
// Use with a binary input image having black particles on white background
// (the macro does not care whether "black background" is selected in Process>Binary>Options)
//
// Needs the "Radial Profile" plugin, http://rsb.info.nih.gov/ij/plugins/radial-profile.html
//
// Limitations:
// - Particle positions are rounded to nearest full pixel.
// - Distance is in pixels, irrespective of any spatial calibration of the image
// - Particles touching the edge will be ignored; this will limit the accuracy
//   if the particles are not much smaller than the image size.
// - Do not extend the image size for avoiding edge effects; the macro takes care of this.
//
// Version history:
// 2008-Dec-09 Michael Schmid - first version
// 2009-May-14 Michael Schmid -  "select none", more specific comments
//
run("Select None");
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;
//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 positions
  run("Find Maxima...", "noise=10 output=[Single Points] light exclude");
  tempID=getImageID();
  tempTitle="temp-"+random();
  rename(tempTitle);
  run("Canvas Size...", "width="+ size+" height="+ size+" position=Center zero");
  run("FD Math...", "image1=["+tempTitle+"] operation=Correlate image2=["+tempTitle+"] result=AutoCorrelation do");
  psID=getImageID();
  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);
  makeRectangle(size/2, size/2, 1, 1);
  run("Set...", "value=0");
  run("Select None");
  circleSize=2*floor(maxRadius)+1;
  run("Specify...", "width="+circleSize+" height="+circleSize+" x="+(size/2+0.5)+" y="+(size/2+0.5)+" oval centered");
  getRawStatistics(nPixels, mean);
  run("Select None");
  run("Divide...", "value="+mean);
  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);
  rename("RDF of "+title);
  selectImage(psID);
  close();
  if (doStack) {
    Plot.getValues(x, y);
    if (slice==firstSlice) ySum = newArray(y.length);
    for (i=0; i<y.length; i++)
      ySum[i]+=y[i]/nSlices;
    close();
  }
}
if (doStack) {
  Plot.create("RDF of "+" title (stack)", "Distance (pixels)", "RDF", x, ySum);
}
setBatchMode("exit and display");