La forma más eficiente de calcular el perfil radial.

Necesito optimizar esta parte de una aplicación de procesamiento de imágenes.
Es básicamente la sum de los píxeles agrupados por su distancia desde el punto central.

def radial_profile(data, center): y,x = np.indices((data.shape)) # first determine radii of all pixels r = np.sqrt((x-center[0])**2+(y-center[1])**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)[0] # 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 radialprofile = tbin/nr # the answer return radialprofile img = plt.imread('crop.tif', 0) # center, radi = find_centroid(img) center, radi = (509, 546), 55 rad = radial_profile(img, center) plt.plot(rad[radi:]) plt.show() 

Imagen de entrada: introduzca la descripción de la imagen aquí

El perfil radial: introduzca la descripción de la imagen aquí

Al extraer los picos de la ttwig resultante, puedo encontrar con precisión los radios de los anillos externos, que es el objective final aquí.

Edit: para más referencia voy a publicar mi solución final de esto. Usando cython obtuve una aceleración de 15-20x en comparación con la respuesta aceptada.

 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 =(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[0], img.shape[1] cdef int topLeft, topRight, botLeft, botRight topLeft =  ceil(sqrt(centerX**2 + centerY**2)) topRight =  ceil(sqrt((xs - centerX)**2 + (centerY)**2)) botLeft =  ceil(sqrt(centerX**2 + (ys-centerY)**2)) botRight =  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 

Parece que podrías usar numpy.bincount aquí:

 import numpy as np def radial_profile(data, center): y, x = np.indices((data.shape)) r = np.sqrt((x - center[0])**2 + (y - center[1])**2) r = r.astype(np.int) tbin = np.bincount(r.ravel(), data.ravel()) nr = np.bincount(r.ravel()) radialprofile = tbin / nr return radialprofile 

Puede usar numpy.histogram para sumr todos los píxeles que aparecen en un “anillo” dado (rango de valores de r desde el origen). Cada anillo es uno de los contenedores de histogtwig. Eliges la cantidad de bandejas dependiendo de qué tan ancho quieres que sean los anillos. (Aquí encontré que los anillos de 3 píxeles de ancho funcionan bien para hacer que la ttwig no sea demasiado ruidosa).

 def radial_profile(data, center): y,x = np.indices((data.shape)) # first determine radii of all pixels r = np.sqrt((x-center[0])**2+(y-center[1])**2) # radius of the image. r_max = np.max(r) ring_brightness, radius = np.histogram(r, weights=data, bins=r_max/3) plt.plot(radius[1:], ring_brightness) plt.show() 

(Por cierto, si esto realmente necesita ser eficiente, y hay muchas imágenes del mismo tamaño, entonces todo antes de la llamada a np.histogram puede ser precomputado).

Tomado de una propuesta de mejora numpy en la que estoy trabajando:

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

La llamada a la media devuelve los valores únicos en R, y la media de los valores correspondientes en datos sobre valores idénticos en R.

Así que no es lo mismo que una solución basada en histogtwig; no tiene que volver a asignar una nueva cuadrícula, lo cual es bueno si desea ajustar un perfil radial, sin pérdida de información. El rendimiento debe ser ligeramente mejor que su solución original. Además, las desviaciones estándar se pueden calcular con la misma eficiencia.

Aquí está mi último borrador númpy group_by EP ; No es una respuesta muy concisa como tal, sino muy general. Espero que todos estemos de acuerdo en que numpy necesita algo como np.group_by (keys) .reduce (values); Si tiene algún comentario, sería bienvenido.

Hay una función en PyDIP que hace esto: dip.RadialMean . Puede usarlo de manera similar a la función radial_profile de OP:

 import PyDIP as dip img = dip.ImageReadTIFF('crop.tif') # center, radi = find_centroid(img) center, radi = (509, 546), 55 rad = dip.RadialMean(img, binSize=1, center=center) rad[radi:].Show() 

Descargo de responsabilidad: soy un autor de PyDIP y la biblioteca DIPlib.