计算径向剖面的最有效方法

Question

计算径向剖面的最有效方法

pythonnumpyimage-processingoptimizationscipy

18

我需要优化这个图像处理应用程序的这一部分。
基本上它是将像素按其与中心点的距离进行分组，并对它们进行求和。

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

输入图像：

径向剖面：这里输入图片描述

通过提取结果图中的峰值，可以准确地找到外环的半径，这也是最终的目标。

编辑：为了进一步参考，我会发布我的最终解决方案。使用cython，我相对于被接受的答案获得了约15-20倍的加速。

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

- M4rtini

5个回答

4

在 DIPlib 中有一个与此相应的函数：dip.RadialMean。您可以类似于 OP 的 radial_profile 函数来使用它：

import diplib 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()

_{免责声明：我是DIPlib库的作者。}

- Cris Luengo

这个解决方案运行速度快，而且很容易扩展到其他统计数据（最小值、最大值、总和）。看起来对于新版本我们应该使用 import diplib as dip。 - ucsky

@ucsky 确实，还有 dip.RadialSum、dip.RadialMaximum 和 dip.RadialMinimum。已更新帖子以反映新软件包名称（在 PyPI 中 "pydip" 已被占用）。 - Cris Luengo

3

您可以使用numpy.histogram来累加出现在给定“环”（从原点开始的r值范围）中的所有像素。每个环都是直方图箱之一。您可以根据需要将环的宽度选择为不同数量的箱。（在这里，我发现3个像素宽的环效果很好，可以使绘图不会太嘈杂。）

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

顺便说一句，如果这确实需要高效处理，并且有很多相同大小的图像，则可以在调用np.histogram之前预计算所有内容。

- mattsh

1

以下内容摘自我正在撰写的一个numpy增强提案：

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

mean函数用于返回R中的唯一值，并计算数据中相应值在R中相同值的平均值。

因此，与基于直方图的解决方案不完全相同；您不必重新映射到新网格，这对于想要拟合径向剖面而不丢失信息的情况很有用。性能应比原始解决方案略好。此外，标准偏差可以以相同的效率计算。

这是我最新的草案 numpy group_by EP；作为一个非常普遍的回答，它并不是非常简洁。我希望我们都同意numpy需要类似np.group_by(keys).reduce(values)的东西；如果您有任何反馈，欢迎提出。

- Eelco Hoogendoorn

0

另一种方法是使用Photutils来处理这类任务，其中您需要获取一个方位角平均的径向剖面。

from photutils.profiles import RadialProfile

img = plt.imread('crop.tif', 0)
center = (509, 546)
edge_radii = np.arange(55)
rp = RadialProfile(data, center, edge_radii)

rp.plot(label='Radial Profile')

- Bhavesh Rajpoot

1

你的回答可以通过提供更多支持性信息来改进。请编辑以添加进一步的细节，例如引用或文档，以便他人能够确认你的回答是否正确。你可以在帮助中心找到关于如何撰写好回答的更多信息。 - Community

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Bi Rico · Accepted Answer

看起来你可以在这里使用 numpy.bincount：

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