我正在尝试用 python 开发一种快速算法，用于查找图像中的峰值，然后找到这些峰值的质心。我使用 scipy.ndimage.label 和 ndimage.find_objects 编写了以下代码来定位对象。这似乎是代码中的瓶颈，在 500x500 图像中定位 20 个对象大约需要 7 毫秒。我想将其放大到更大的 (2000x2000) 图像，但时间会增加到几乎 100 毫秒。所以，我想知道是否有更快的选择。
这是我到目前为止的代码，它可以工作，但速度很慢。首先，我使用一些高斯峰值模拟我的数据。这部分很慢，但实际上我将使用真实数据，所以我不太关心加快这部分的速度。我希望能够很快找到峰值。
import time
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage
import matplotlib.patches
plt.figure(figsize=(10,10))
ax1 = plt.subplot(221)
ax2 = plt.subplot(222)
ax3 = plt.subplot(223)
ax4 = plt.subplot(224)
size = 500 #width and height of image in pixels
peak_height = 100 # define the height of the peaks
num_peaks = 20
noise_level = 50
threshold = 60
np.random.seed(3)
#set up a simple, blank image (Z)
x = np.linspace(0,size,size)
y = np.linspace(0,size,size)
X,Y = np.meshgrid(x,y)
Z = X*0
#now add some peaks
def gaussian(X,Y,xo,yo,amp=100,sigmax=4,sigmay=4):
return amp*np.exp(-(X-xo)**2/(2*sigmax**2) - (Y-yo)**2/(2*sigmay**2))
for xo,yo in size*np.random.rand(num_peaks,2):
widthx = 5 + np.random.randn(1)
widthy = 5 + np.random.randn(1)
Z += gaussian(X,Y,xo,yo,amp=peak_height,sigmax=widthx,sigmay=widthy)
#of course, add some noise:
Z = Z + scipy.ndimage.gaussian_filter(0.5*noise_level*np.random.rand(size,size),sigma=5)
Z = Z + scipy.ndimage.gaussian_filter(0.5*noise_level*np.random.rand(size,size),sigma=1)
t = time.time() #Start timing the peak-finding algorithm
#Set everything below the threshold to zero:
Z_thresh = np.copy(Z)
Z_thresh[Z_thresh<threshold] = 0
print 'Time after thresholding: %.5f seconds'%(time.time()-t)
#now find the objects
labeled_image, number_of_objects = scipy.ndimage.label(Z_thresh)
print 'Time after labeling: %.5f seconds'%(time.time()-t)
peak_slices = scipy.ndimage.find_objects(labeled_image)
print 'Time after finding objects: %.5f seconds'%(time.time()-t)
def centroid(data):
h,w = np.shape(data)
x = np.arange(0,w)
y = np.arange(0,h)
X,Y = np.meshgrid(x,y)
cx = np.sum(X*data)/np.sum(data)
cy = np.sum(Y*data)/np.sum(data)
return cx,cy
centroids = []
for peak_slice in peak_slices:
dy,dx = peak_slice
x,y = dx.start, dy.start
cx,cy = centroid(Z_thresh[peak_slice])
centroids.append((x+cx,y+cy))
print 'Total time: %.5f seconds\n'%(time.time()-t)
###########################################
#Now make the plots:
for ax in (ax1,ax2,ax3,ax4): ax.clear()
ax1.set_title('Original image')
ax1.imshow(Z,origin='lower')
ax2.set_title('Thresholded image')
ax2.imshow(Z_thresh,origin='lower')
ax3.set_title('Labeled image')
ax3.imshow(labeled_image,origin='lower') #display the color-coded regions
for peak_slice in peak_slices: #Draw some rectangles around the objects
dy,dx = peak_slice
xy = (dx.start, dy.start)
width = (dx.stop - dx.start + 1)
height = (dy.stop - dy.start + 1)
rect = matplotlib.patches.Rectangle(xy,width,height,fc='none',ec='red')
ax3.add_patch(rect,)
ax4.set_title('Centroids on original image')
ax4.imshow(Z,origin='lower')
for x,y in centroids:
ax4.plot(x,y,'kx',ms=10)
ax4.set_xlim(0,size)
ax4.set_ylim(0,size)
plt.tight_layout
plt.show()
The results for size=500:
编辑：如果峰值数量很大（~100）并且图像尺寸很小，那么瓶颈实际上是质心部分。所以，也许这部分的速度也需要优化。
当然，您寻找峰值的方法（简单阈值处理）对阈值的选择非常敏感：将阈值设置得太低，您将“检测到”不是峰值的东西；将其设置得太高，您将错过有效的峰值。
有更强大的替代方案，可以检测图像强度中的所有局部最大值，无论其强度值如何。我的首选方法是使用小型（5x5 或 7x7）结构元素进行扩张，然后找到原始图像及其扩张版本具有相同值的像素。这是有效的，因为根据定义， dilation(x, y, E, img) = { E 内以像素 (x,y) 为中心的最大 img }，因此 dilation(x, y, E, img) = img(x , y) 只要 (x,y) 是 E 尺度上的局部最大值的位置。
通过快速实现形态学算子（例如 OpenCV 中的算子），该算法在空间和时间上的图像大小都是线性的（一个额外的图像大小的缓冲区用于扩张图像，并且对两者都进行一次传递）。在紧要关头，它也可以在线实现，无需额外的缓冲区和一点额外的复杂性，而且仍然是线性时间。
为了在存在椒盐或类似噪声（可能会引入许多虚假最大值）的情况下进一步增强其鲁棒性，您可以使用不同大小的结构元素（例如 5x5 和 7x7）应用该方法两次，然后仅保留稳定的最大值，其中稳定性可以通过最大值的位置不变或位置变化不超过一个像素等来定义。此外，当您有理由相信附近的最大值是由噪声引起时，您可能希望抑制它们。一种有效的方法是首先检测上面的所有局部最大值，按高度降序对它们进行排序，然后沿着排序列表向下查找，如果它们在图像中的值没有改变，则保留它们，如果保留，则设置为将 (2d+1) x (2d+1) 邻域中的所有像素归零，其中 d 是您愿意容忍的附近最大值之间的最小距离。