# Calculate a radial distribution function from a 2D gray scale image

Hi,

I am trying to compute a radial distribution function (or pair distribution function) from a 2D gray scale image (it could be a binary image).

The radial distribution function is a kind of a simple sum of two points products of all the pixels on the image.

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3     // Use modern global access method and strict wave access.

wave pWave

variable Xmax, Ymax
Xmax = DimSize(pWave,0) - 1
Ymax = DimSize(pWave,1) - 1

variable r_max = round(sqrt(Xmax^2 + Ymax^2))

variable reftime, seconds
reftime = StartMSTimer

variable i, j
for(i=0; i<5; i++) // Idealy, I need to put Xmax and Ymax here.
for (j=0; j<5; j++)
GetCircularAverage(pWave, i, j)
wave temp_profile
endfor
endfor

seconds = StopMSTimer(reftime)/1e6
Print seconds

End

/// Circular average loop
Function GetCircularAverage(pWave, X0, Y0)
Wave pWave
variable X0, Y0

variable Xmax, Ymax
Xmax = DimSize(pWave,0) - 1
Ymax = DimSize(pWave,1) - 1

variable r_max = round(sqrt(Xmax^2 + Ymax^2))

/// Create a vector map from the perdefined ceneter
Make/FREE/D/O/N=(Xmax+1,Ymax+1,2) pvecMap
Make/FREE/D/O/N=2 xvec, yvec
xvec = {1,0}
yvec = {0,1}
Multithread pvecMap = (p-X0)*xvec[r] + (q-Y0)*yvec[r]
MatrixOP/FREE/O pscalarMap = round( sqrt( pvecMap[][][0]*pvecMap[][][0] + pvecMap[][][1]*pvecMap[][][1]) )

/// Create a refwave to store values. added int and count number.
make/FREE/D/O/N=(r_max+1) ref_yy, yy, count

/// Start circular average
variable i, j, r
count = 0 //initialize count
/// Add pixels to get circular sum
for(i=0; i<Xmax+1; i++) //gXmax is the coordinates.
for (j=0; j<Ymax+1; j++)

r = pscalarMap[i][j]
ref_yy[r] += pWave[i][j]
count[r] += 1 //calculate the pixel number that corresponds to distance r, considering the mask.

endfor
endfor

yy = ref_yy/count
Duplicate/O yy, temp_profile

End

Currently, I compute a radial distribution from a selected pixel as the center and then move to next pixel and ...

After calculating the radial distribution functions from all the pixels, then I add all these profiles together.

It works if the image is small, but it could take hours to complete an image with 1000x1000 pixels.

Is there a better way to do this calculation.

I haven't looked at your code in detail, but it is always useful to time the different steps in your code to figure out where the bottleneck is

variable t=StartMSTimer
<do something>
Print(StopMSTimer(t))

You could probably also calculate the distances between pixels once, instead of redoing it for each pixel. Say you have a 100 x 100 pixel image. You could expand that to a 199 x 199 pixel image and calculate the distance of each pixel to the center. Later on you then simply duplicate the relevant 100 x 100 pixel subsection of the 199 x 199 image to get the distance of all pixels to your selected center.

I'm sure there is a lot of tricks you can play, but it really helps to know where the bottleneck is before you start wasting time

In addition to the above, it looks to me that a number of parameters are always identical in each call of "GetCircularAverage", like e.g. Xmax, Ymax, pvecMap and are always evaluated. This could all be precalculated once and then passed to the function.

Apart from what has been said, your algorithm will approximately scale like O(sidelength^2), so the computational time will blow up eventually.

If in a typical scenario the RDF approaches zero for large r (I guess it does for most applications), you could save a lot of computational time by only calculating the circular average up to a certain r_max. After r_max it's mostly zero anyway.

And also remove code which does nothing useful. radial_r  is calculated in GetCircularAverage but never used. Ideally you get GetCircularAverage into a state where it does only calculate stuff but not make/duplicate/redimension waves. And execute it on all cores using multithread.

Would this help.

https://www.wavemetrics.com/node/21289

Or the discussion here.