23

我试图python尽可能接近图像中最明显的聚类中心,如下图所示:

图片

在我之前的问题中,我询问了如何获得二维数组的全局最大值和局部最大值,并且给出的答案非常有效。问题是我可以通过对不同 bin 大小获得的全局最大值进行平均得到的中心估计值总是比我用眼睛设置的值稍有偏差,因为我只考虑最大的bin而不是一最大的 bin (就像用眼睛做的那样)。

我尝试根据我的问题调整这个问题的答案,但事实证明我的图像太嘈杂,以至于该算法无法工作。这是我实现该答案的代码:

import numpy as np
from scipy.ndimage.filters import maximum_filter
from scipy.ndimage.morphology import generate_binary_structure, binary_erosion
import matplotlib.pyplot as pp

from os import getcwd
from os.path import join, realpath, dirname

# Save path to dir where this code exists.
mypath = realpath(join(getcwd(), dirname(__file__)))
myfile = 'data_file.dat'

x, y = np.loadtxt(join(mypath,myfile), usecols=(1, 2), unpack=True)
xmin, xmax = min(x), max(x)
ymin, ymax = min(y), max(y)

rang = [[xmin, xmax], [ymin, ymax]]
paws = []

for d_b in range(25, 110, 25):
    # Number of bins in x,y given the bin width 'd_b'
    binsxy = [int((xmax - xmin) / d_b), int((ymax - ymin) / d_b)]

    H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)
    paws.append(H)


def detect_peaks(image):
    """
    Takes an image and detect the peaks usingthe local maximum filter.
    Returns a boolean mask of the peaks (i.e. 1 when
    the pixel's value is the neighborhood maximum, 0 otherwise)
    """

    # define an 8-connected neighborhood
    neighborhood = generate_binary_structure(2,2)

    #apply the local maximum filter; all pixel of maximal value 
    #in their neighborhood are set to 1
    local_max = maximum_filter(image, footprint=neighborhood)==image
    #local_max is a mask that contains the peaks we are 
    #looking for, but also the background.
    #In order to isolate the peaks we must remove the background from the mask.

    #we create the mask of the background
    background = (image==0)

    #a little technicality: we must erode the background in order to 
    #successfully subtract it form local_max, otherwise a line will 
    #appear along the background border (artifact of the local maximum filter)
    eroded_background = binary_erosion(background, structure=neighborhood, border_value=1)

    #we obtain the final mask, containing only peaks, 
    #by removing the background from the local_max mask
    detected_peaks = local_max - eroded_background

    return detected_peaks


#applying the detection and plotting results
for i, paw in enumerate(paws):
    detected_peaks = detect_peaks(paw)
    pp.subplot(4,2,(2*i+1))
    pp.imshow(paw)
    pp.subplot(4,2,(2*i+2) )
    pp.imshow(detected_peaks)

pp.show()

这是结果(改变 bin 大小):

在此处输入图像描述

显然,我的背景太嘈杂以至于该算法无法工作,所以问题是:我怎样才能使该算法不那么敏感?如果存在替代解决方案,请告诉我。


编辑

按照 Bi Rico 的建议,我尝试在将二维数组传递给局部最大值查找器之前对其进行平滑处理,如下所示:

H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)
H1 = gaussian_filter(H, 2, mode='nearest')
paws.append(H1)

这些是 asigma为 2、4 和 8 的结果:

在此处输入图像描述

编辑 2

Amode ='constant'似乎比nearest. sigma=2对于最大的 bin 大小,它以 a 收敛到右中心:

在此处输入图像描述

那么,如何获得最后一张图像中显示的最大值的坐标?

4

5 回答 5

4

回答问题的最后一部分,图像中总是有点,您可以通过按某种顺序搜索图像的局部最大值来找到它们的坐标。如果您的数据不是点源,您可以对每个峰值应用掩码,以避免在执行未来搜索时峰值邻域成为最大值。我建议以下代码:

import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
import copy

def get_std(image):
    return np.std(image)

def get_max(image,sigma,alpha=20,size=10):
    i_out = []
    j_out = []
    image_temp = copy.deepcopy(image)
    while True:
        k = np.argmax(image_temp)
        j,i = np.unravel_index(k, image_temp.shape)
        if(image_temp[j,i] >= alpha*sigma):
            i_out.append(i)
            j_out.append(j)
            x = np.arange(i-size, i+size)
            y = np.arange(j-size, j+size)
            xv,yv = np.meshgrid(x,y)
            image_temp[yv.clip(0,image_temp.shape[0]-1),
                                   xv.clip(0,image_temp.shape[1]-1) ] = 0
            print xv
        else:
            break
    return i_out,j_out

#reading the image   
image = mpimg.imread('ggd4.jpg')
#computing the standard deviation of the image
sigma = get_std(image)
#getting the peaks
i,j = get_max(image[:,:,0],sigma, alpha=10, size=10)

#let's see the results
plt.imshow(image, origin='lower')
plt.plot(i,j,'ro', markersize=10, alpha=0.5)
plt.show()

用于测试的图像 ggd4 可以从以下位置下载:

http://www.ipac.caltech.edu/2mass/gallery/spr99/ggd4.jpg

第一部分是获取有关图像中噪声的一些信息。我通过计算整个图像的标准偏差来做到这一点(实际上最好选择一个没有信号的小矩形)。这告诉我们图像中存在多少噪声。获得峰值的想法是要求连续的最大值,这些最大值高于某个阈值(例如,噪声的 3、4、5、10 或 20 倍)。这就是函数 get_max 实际在做的事情。它执行最大值搜索,直到其中一个低于噪声施加的阈值。为了避免多次找到相同的最大值,有必要从图像中删除峰值。一般来说,面具的形状很大程度上取决于一个人想要解决的问题。对于星星的情况,使用高斯函数或类似的方法去除星星应该很好。为简单起见,我选择了一个平方函数,函数的大小(以像素为单位)是变量“大小”。我认为从这个例子中,任何人都可以通过添加更通用的东西来改进代码。

编辑:

原始图像如下所示:

在此处输入图像描述

而识别出光点后的图像是这样的:

在此处输入图像描述

于 2013-07-04T04:19:15.747 回答
4

Stack Overflow 上的 n00b 太多了,无法在此处其他地方评论 Alejandro 的回答。我会稍微改进他的代码以使用预先分配的 numpy 数组进行输出:

def get_max(image,sigma,alpha=3,size=10):
    from copy import deepcopy
    import numpy as np
    # preallocate a lot of peak storage
    k_arr = np.zeros((10000,2))
    image_temp = deepcopy(image)
    peak_ct=0
    while True:
        k = np.argmax(image_temp)
        j,i = np.unravel_index(k, image_temp.shape)
        if(image_temp[j,i] >= alpha*sigma):
            k_arr[peak_ct]=[j,i]
            # this is the part that masks already-found peaks.
            x = np.arange(i-size, i+size)
            y = np.arange(j-size, j+size)
            xv,yv = np.meshgrid(x,y)
            # the clip here handles edge cases where the peak is near the 
            #    image edge
            image_temp[yv.clip(0,image_temp.shape[0]-1),
                               xv.clip(0,image_temp.shape[1]-1) ] = 0
            peak_ct+=1
        else:
            break
    # trim the output for only what we've actually found
    return k_arr[:peak_ct]

在使用他的示例图像分析此代码和 Alejandro 的代码时,此代码快约 33%(Alejandro 的代码为 0.03 秒,我的代码为 0.02 秒。)我希望在具有大量峰值的图像上,它会更快 - 将输出附加到对于更多的峰值,列表将变得越来越慢。

于 2013-08-04T05:08:16.710 回答
2

我认为这里需要的第一步是用字段的标准偏差来表示 H 中的值:

import numpy as np
H = H / np.std(H)

现在你可以在这个 H 的值上设置一个阈值。如果假设噪声是高斯的,选择一个阈值 3,你可以非常确定(99.7%)这个像素可以与一个真正的峰值而不是噪声相关联。见这里

现在可以开始进一步的选择。我不清楚你到底想找到什么。您想要峰值的确切位置吗?或者您想要一个位于该集群中间的一组峰的位置?
无论如何,从这一点开始,所有像素值都以场的标准偏差表示,你应该能够得到你想要的。如果要查找集群,可以对 >3-sigma 网格点执行最近邻搜索,并在距离上设置阈值。即仅当它们彼此足够接近时才连接它们。如果连接了多个网格点,您可以将其定义为组/集群并计算集群的一些(sigma-weighted?)中心。希望我对 Stackoverflow 的第一个贡献对您有用!

于 2013-07-03T22:12:00.120 回答
1

我会这样做的方式:

1) 将 H 归一化在 0 和 1 之间。

2) 选择一个阈值,正如 tcaswell 所建议的那样。例如,它可能介于 0.9 和 0.99 之间

3) 使用掩码数组仅保留 H 高于阈值的 x,y 坐标:

import numpy.ma as ma
x_masked=ma.masked_array(x, mask= H < thresold)
y_masked=ma.masked_array(y, mask= H < thresold)

4) 现在您可以对蒙面坐标进行加权平均,权重类似于 (H-threshold)^2,或任何其他大于或等于 1 的幂,具体取决于您的口味/测试。

评论:1)这对于您拥有的峰类型并不可靠,因为您可能必须调整阈值。这是小问题;2)这不适用于两个峰值,如果第二个峰值高于阈值,则会给出错误的结果。

尽管如此,它总是会给你一个答案而不会崩溃(有事情的优点和缺点..)

于 2013-05-30T21:19:01.153 回答
1

我添加了这个答案,因为它是我最终使用的解决方案。这是 Bi Rico 在此处的评论(5 月 30 日 18:54)和此问题中给出的答案的组合:Find peak of 2d histogram

事实证明,使用此问题中的峰值检测算法2D 阵列中的峰值检测只会使问题复杂化。在对图像应用高斯滤波器之后,需要做的就是要求最大 bin(正如 Bi Rico 指出的那样),然后获得坐标中的最大值。

因此,我没有像上面那样使用检测峰值函数,而是在获得高斯 2D 直方图后简单地添加以下代码:

# Get 2D histogram.
H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)
# Get Gaussian filtered 2D histogram.
H1 = gaussian_filter(H, 2, mode='nearest')
# Get center of maximum in bin coordinates.
x_cent_bin, y_cent_bin = np.unravel_index(H1.argmax(), H1.shape)
# Get center in x,y coordinates.
x_cent_coor , y_cent_coord = np.average(xedges[x_cent_bin:x_cent_bin + 2]), np.average(yedges[y_cent_g:y_cent_g + 2])
于 2013-07-04T13:09:59.703 回答