# 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\limits_{x=1}^{n_x} \sum\limits_{y=1}^{n_y} h\left(\frac{d_{x,y}-r}{s}\right) \cdot \text{img}[x,y] }{\sum_\limits{x=1}^{n_x} \sum_\limits{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.