# Most efficient way to calculate radial profile

Posted on

### Question :

Most efficient way to calculate radial profile

I need to optimize this part of an image processing application.
It is basically the sum of the pixels binned by their distance from the central spot.

``````def radial_profile(data, center):
y,x = np.indices((data.shape)) # first determine radii of all pixels
r = np.sqrt((x-center)**2+(y-center)**2)
ind = np.argsort(r.flat) # get sorted indices
sr = r.flat[ind] # sorted radii
sim = data.flat[ind] # image values sorted by radii
ri = sr.astype(np.int32) # integer part of radii (bin size = 1)
# determining distance between changes
deltar = ri[1:] - ri[:-1] # assume all radii represented
rind = np.where(deltar) # location of changed radius
nr = rind[1:] - rind[:-1] # number in radius bin
csim = np.cumsum(sim, dtype=np.float64) # cumulative sum to figure out sums for each radii bin
tbin = csim[rind[1:]] - csim[rind[:-1]] # sum for image values in radius bins

center, radi = (509, 546), 55

plt.show()
`````` By extracting the peaks of the resulting plot, I can accurately find the radii of the outer rings, which is the end goal here.

Edit: For further reference I’ll post my final solution of this. Using `cython` I got about a 15-20x speed up compared to the accepted answer.

``````import numpy as np
cimport numpy as np
cimport cython
from cython.parallel import prange
from libc.math cimport sqrt, ceil

DTYPE_IMG = np.uint8
ctypedef np.uint8_t DTYPE_IMG_t
DTYPE = np.int
ctypedef np.int_t DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cdef void cython_radial_profile(DTYPE_IMG_t [:, :] img_view, DTYPE_t [:] r_profile_view, int xs, int ys, int x0, int y0) nogil:

cdef int x, y, r, tmp

for x in prange(xs):
for y in range(ys):
r =<int>(sqrt((x - x0)**2 + (y - y0)**2))
tmp = img_view[x, y]
r_profile_view[r] +=  tmp

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def radial_profile(np.ndarray img, int centerX, int centerY):

cdef int xs, ys, r_max
xs, ys = img.shape, img.shape

cdef int topLeft, topRight, botLeft, botRight

topLeft = <int> ceil(sqrt(centerX**2 + centerY**2))
topRight = <int> ceil(sqrt((xs - centerX)**2 + (centerY)**2))
botLeft = <int> ceil(sqrt(centerX**2 + (ys-centerY)**2))
botRight = <int> ceil(sqrt((xs-centerX)**2 + (ys-centerY)**2))

r_max = max(topLeft, topRight, botLeft, botRight)

cdef np.ndarray[DTYPE_t, ndim=1] r_profile = np.zeros([r_max], dtype=DTYPE)
cdef DTYPE_t [:] r_profile_view = r_profile
cdef DTYPE_IMG_t [:, :] img_view = img

with nogil:
cython_radial_profile(img_view, r_profile_view, xs, ys, centerX, centerY)
return r_profile
``````

It looks like you could use numpy.bincount here:

``````import numpy as np

y, x = np.indices((data.shape))
r = np.sqrt((x - center)**2 + (y - center)**2)
r = r.astype(np.int)

tbin = np.bincount(r.ravel(), data.ravel())
nr = np.bincount(r.ravel())
``````

You can use numpy.histogram to add up all the pixels that appear in a given “ring” (range of values of r from the origin). Each ring is one of the histogram bins. You choose the number of bins depending on how wide you want the rings to be. (Here I found 3 pixel wide rings work well to make the plot not too noisy.)

``````def radial_profile(data, center):
y,x = np.indices((data.shape)) # first determine radii of all pixels
r = np.sqrt((x-center)**2+(y-center)**2)

r_max = np.max(r)

ring_brightness, radius = np.histogram(r, weights=data, bins=r_max/3)
plt.show()
``````

(By the way, if this really needs to be efficient, and there are a lot of images the same size, then everything before the call to np.histogram can be precomputed.)

There is a function in PyDIP that does just this: `dip.RadialMean`. You can use it in a similar way to OP’s `radial_profile` function:

``````import PyDIP as dip

center, radi = (509, 546), 55

``````

Disclaimer: I’m an author of PyDIP and the DIPlib library.

``````pp.plot(*group_by(np.round(R, 5).flatten()).mean(data.flatten()))