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

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

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.

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.

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.

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.

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.