计算径向轮廓的最有效方法

Posted

技术标签:

【中文标题】计算径向轮廓的最有效方法【英文标题】:Most efficient way to calculate radial profile 【发布时间】:2014-02-10 02:17:44 【问题描述】:

我需要优化图像处理应用程序的这一部分。 它基本上是像素与中心点的距离的总和。

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

【问题讨论】:

【参考方案1】:

看起来你可以在这里使用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 

【讨论】:

不错的解决方案,大约快 3 倍。在我接受之前我会稍等片刻,永远不知道人们能想出什么。 我在这个实现中使用了 bincount 和 histogram:github.com/keflavich/image_tools/blob/master/image_tools/…【参考方案2】:

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 库的作者。

【讨论】:

这个解决方案运行速度很快,很容易扩展到其他统计数据(最小值、最大值、总和)。对于新版本,我们应该使用import diplib as dip @ucsky 实际上,还有dip.RadialSumdip.RadialMaximumdip.RadialMinimum。更新了帖子以反映新的包名称(“pydip”已在 PyPI 中使用)。【参考方案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 之前的所有内容。)

【讨论】:

【参考方案4】:

取自我正在处理的一个 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); 这样的东西。如果您有任何反馈,我们将不胜感激。

【讨论】:

以上是关于计算径向轮廓的最有效方法的主要内容,如果未能解决你的问题,请参考以下文章

在 Python 中创建大型数组 [关闭]

从轮廓中提取最外层轮廓

用opencv python填充轮廓

数据分析 第六篇:聚类的评估(簇数确定和轮廓系数)和可视化

计算机视觉OpenCV篇 - 轮廓(寻找/绘制轮廓)

OpenCV图像的轮廓的匹配