python - optimize - La manera más eficiente de calcular el perfil radial

scipy.optimize python

Question

Necesito optimizar esta parte de una aplicación de procesamiento de imágenes.
Básicamente es la suma 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:

El perfil radial:

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

Editar: para una referencia posterior publicaré mi solución final de esto. Usando cython obtuve una velocidad 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 =<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[0], img.shape[1]
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

Answer 1

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 

Answer 2

Puede usar numpy.histogram para sumar todos los píxeles que aparecen en un "anillo" dado (rango de valores de r del origen). Cada anillo es uno de los contenedores de histograma. Usted elige la cantidad de contenedores según el ancho que desee que sean los anillos. (Aquí encontré anillos de 3 píxeles de ancho que funcionan bien para que la trama 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 se puede precalcular).

Answer 3

Tomado de una propuesta de mejora numpy Estoy trabajando en:

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

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

Entonces no es exactamente lo mismo que una solución basada en un histograma; no es necesario reasignar a una nueva grilla, lo que es bueno si desea ajustar un perfil radial, sin pérdida de información. El rendimiento debería ser ligeramente mejor que su solución original. Además, las desviaciones estándar se pueden calcular con la misma eficacia.

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