# 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).

https://www.wikiwand.com/en/Radial_distribution_function

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

#pragma rtGlobals=3 // Use modern global access method and strict wave access.

//Function Get Radial distribution

Function GetRadialDistribution(pWave)

wave pWave

variable Xmax, Ymax

Xmax = DimSize(pWave,0) - 1

Ymax = DimSize(pWave,1) - 1

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

make/O/N=(r_max+1)/O radial_r = p

make/O/N=(r_max+1)/O radial_i = 0

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

radial_i += 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))

make/FREE/N=(r_max+1)/O radial_r = p

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

<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

October 19, 2020 at 10:18 pm - Permalink

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.

October 19, 2020 at 11:30 pm - Permalink

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.

October 20, 2020 at 12:23 am - Permalink

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.

October 20, 2020 at 01:34 am - Permalink

Would this help.

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

Or the discussion here.

https://www.wavemetrics.com/code-snippet/radial-profiler

EDIT: You may be after a different function than just a radial profile from a point on the image. I have to wonder indeed whether what you want is equivalent to taking the 2D Fourier transform of the 2D image and then taking a radial profile measured from the central point of the 2D Fourier transform.

October 20, 2020 at 06:15 am - Permalink

Fourier transforms look for long range ordering. I think what Maru is looking for is the local, short range ordering. It's like the difference between looking at the structure of a crystal versus something amorphous.

October 20, 2020 at 12:44 pm - Permalink

I cannot help but think that we should ask how we would do measurements on the image to obtain the specific results that are needed and then reproduce the measurement process using first principles.

I have in this regard a picture of something akin to Bragg scattering from a 2-D lattice (low energy electron diffraction) where Dirac delta functions will broaden as the fundamental surface structure goes from perfect crystalline to glassy to amorphous. Won't a radial profile from the center of the 2D FFT go from a set of Dirac delta lines to a radial distribution plot showing the abundance of local versus long-range coordinations? Or EXAFS where the measured result (as I understand) is akin to a Fourier transform of the radial distribution function for neighboring atoms from any given local state. Or neutron scattering from 2D surfaces.

October 20, 2020 at 02:24 pm - Permalink

Thank you very much for many suggestions and ideas. Using a predefined pvecMap is a nice idea and limiting r_max is crucial in calculating large images. I will definitely remove the unnecessary codes and make the calculation part multitrehad. And thank you for giving me example codes for calculating radial distribution. I will follow your advices one by one and test how fast I can make my procedure be. Sorry, but it could take me some time to rewrite the script.

I am going to analyze a 2D real space image from a confocal microscopy, not scattering patterns this time. I can perform a 2D Fourier transform on the image and analyze the correlations in the reciprocal space but I think it is more straightforward to get the radial distribution function directly from the real space images.

October 20, 2020 at 06:42 pm - Permalink