Calculating radial profiles

Hi Gleb

I thought a bit more about how to speed up calculating the radial profiles. At the moment, you have a triple loop, the outer loops goes through the grid of radius values, the inner two loops through the pixels of the image. We should turn this around, as I show below.

Specification

Input:

• an image img as matrix of size nx by ny
• a maximum radius rmax
• coordinates cx, cy of the center of the circles
• a kernel function h(x) that is nonzero only for abs(x) < 1
• a smoothing bandwidth s

Output:

• a vector res, of length rmax, with kernel-smoothed average intensities of the pixels at the distances 1, 2, 3, …, rmax from (cx,cy).

$$\text{res}[r]=\frac{\sum _{x=1}^{n_x}\sum _{y=1}^{n_y}h\left(\frac{d_{x,y}-r}{s}\right)\cdot \text{img}[x,y]}{\sum _{x=1}^{n_x}\sum _{y=1}^{n_y}h\left(\frac{d_{x,y}-r}{s}\right)}$$

with

$d(x,y)=\sqrt{(x-c_x{)}^2+(y-c_y{)}^2}$.

Algorithm:

First the naive form. We simply swap the loops. For this, we need to calculate numerator and denominator of res[i] separately and divide only at the end.

num = Vector of length maxr, initialized with zeroes  # numerators
den = Vector of length maxr, initialized with zeroes  # denominators
for x from 1 to nx:
   for y from 1 to *nx*
      d = distance of (x,y) tp (cx,cy)
      for r from 1 to maxr
         w = h( ( r-d ) / s )  # the kernel weight
         num[i] += w * img[x,y]
         den[i] += w
      end for
   end for
end for
res = num/den (elementwise division)

Now have a look at the inner loop: For most of the values of r, w will be zero, and we could skip all those itertaions of the loop.

In fact, we already know before entering the loop that w can only be nonzero if
abs( r - d ) < s. Therefore, the innermost loop has to run only from d-s to d+s This reduces run-time by a factor of rmax/(2*s), which may well by around 100.