我正在尝试使python尽可能接近地返回图像中最明显的聚类的中心,如下图所示:

在我的previous question中,我问如何获取二维数组的全局最大值和局部最大值,给出的答案非常有效。问题在于,我可以通过平均不同大小的bin所获得的全局最大值来获得的中心估计值总是比我用眼睛设置的中心估计值略有偏离,因为我只考虑了最大的 bin 而不是组最大垃圾箱的(就像一个人用眼睛一样)。

我尝试将answer to this question调整为我的问题,但结果表明我的图像过于嘈杂,无法使用该算法。这是我实现该答案的代码:

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的建议,我尝试先将2d数组平滑,然后再将其传递到本地最大查找器,如下所示:
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更好。它使用sigma=2收敛到右侧中心,以获取最大的bin大小:

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

最佳答案

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秒。)我希望在具有更多峰值的图像上,它甚至会更快-将输出附加到列表会越来越慢,出现更多峰。

关于python - 嘈杂的二维阵列中的峰检测,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/16842823/

10-13 09:14