在嘈杂的二维数组中检测峰值

Question

我试图让python返回一个图像中最明显聚类的中心点，并尽可能接近以下图片:

在我的先前的问题中，我问了如何获取2D数组的全局最大值和局部最大值，给出的答案完美解决了这个问题。问题是，我通过对不同bin size得到的全局最大值取平均值来获得的中心估计总是稍微偏离我会凭眼睛看出的那个中心点，因为我只考虑了最大的bin而没有考虑由多个最大的bin组成的group（就像用肉眼观察一样）。

我试着将这个问题的答案应用到我的问题上，但事实证明我的图像太嘈杂了，无法使用该算法。下面是我实现该答案的代码:

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

以下是通过更改组距所得到的结果:

很明显，我的背景噪声太多了，使得该算法无法正常工作。因此，问题是:我如何使该算法变得更不敏感？如果存在另一种解决方案，请告诉我。

---

编辑

根据Bi Rico的建议，在将二维数组传递给局部最大值查找器之前，我尝试平滑化该数组，方法如下：

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

以下是当 sigma 值分别为 2、4 和 8 时的结果：

编辑 2

采用 mode ='constant' 要比 nearest 效果好得多。在最大 bin 大小为 sigma=2 时，它收敛到正确的中心位置：

那么，我该如何获取出现在最后一张图片中的最大值的坐标？

- Gabriel

---

你尝试过在应用算法之前对数据进行平滑处理吗？高斯滤波和/或中值滤波可能会有所帮助。 - Bi Rico

如何使用 np.unravel_index(array.argmax(), array.shape)？ - Bi Rico

一个简单的阈值也可能会有很大帮助。 - tacaswell

你能描述你试图检测的峰的特性吗？它总是一个单一的峰吗？你期望它是对称的，还是具有特征空间尺度？此外，背景噪声的特性是什么 - 它是否具有空间结构？ - ali_m

这些问题的答案是：我几乎无法描述它们，因为它们非常多变，不总是单个峰值，峰值不对称且没有特定的空间尺度，而背景噪声则相当随机（即：它可以平滑且非常密集，可以成团且难以与真正的“峰”区分开来，也可以非常微弱）。我知道这听起来像一个过于广泛的问题，无法通过“一刀切”的算法方法解决，这就是为什么我很难找到一个相对简单的答案。 - Gabriel

Answer 1

回答你问题的最后一部分，如果你有一张图片上的点，可以通过搜索图片中局部极大值的方式找到它们的坐标。如果你的数据不属于点源，你可以对每个峰值应用掩码，以避免在进行未来搜索时，峰值附近成为最大值。我建议使用以下代码：

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实际上正在做的事情。它执行最大值搜索，直到其中一个最大值低于噪声施加的阈值为止。为了避免多次找到同一个峰值，需要从图像中移除峰值。通常情况下，这样做的蒙版形状强烈依赖于想要解决的问题。对于星星的情况，最好使用高斯函数或类似的东西来去除星星。出于简单起见，我选择了一个正方形函数，并且函数大小（以像素为单位）是变量“size”。

我认为，从这个例子中，任何人都可以通过添加更通用的内容来改进代码。

编辑：

原始图像如下：

而识别出发光点后的图像如下：

- Alejandro

Answer 2

我认为这里需要的第一步是用场的标准差来表示H中的值：

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

现在你可以对这个H的值设置阈值。如果噪声被假设为高斯分布，选择3作为阈值，你可以相当确定（99.7%）这个像素可以与真实峰值相关联而不是噪声。请参见这里。

现在进一步的选择可以开始了。我不太清楚你到底想要找什么。你想要精确的峰值位置吗？还是你想要一个位于这个簇中心的峰值集群的位置？无论如何，从这一点开始，所有像素值都以场的标准差表示，你应该能够得到你想要的东西。如果你想要找到簇，你可以在大于3 sigma的网格点上执行最近邻搜索，并对距离设置一个阈值。即只有当它们彼此足够接近时才连接它们。如果连接了几个网格点，你可以将其定义为一个组/簇并计算一些（sigma加权？）簇的中心。希望我在Stackoverflow上的第一篇文章对你有用！

- Basileios

Answer 3

我在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秒）。我预计在具有更多峰值的图像上，它会更快 - 对输出进行列表附加将对更多的峰值变得越来越慢。

- msarahan

Answer 4

我添加这个答案是因为它是我最终使用的解决方案。它结合了Bi Rico在这里的评论（5月30日18:54）和这个问题中给出的答案：查找2D直方图峰值。

事实证明，仅使用来自此问题中的峰值检测算法 2D数组中的峰值检测只会使事情更加复杂。在对图像应用高斯滤波器后，需要做的就是询问最大的bin（如Bi Rico所指出的），然后获取坐标中的最大值。

因此，不再像上面那样使用detect-peaks函数，而是在获得高斯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])

- Gabriel

Answer 5

我的做法：

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）如果第二个峰值高于阈值，则此方法无法处理两个峰值，会给出错误结果。

尽管如此，它总是会给您一个答案而不会崩溃（具有优缺点）。

- gg349