slsqp optimize opt python optimization numpy scipy

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



scipy.optimize python (3)

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

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


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).


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.