python寻峰算法_python中的快速寻峰与质心

我正在尝试用python开发一个快速算法来查找图像中的峰值，然后找到这些峰值的质心。我已经用scipy.ndimage.label和ndimage.find_对象编写了以下代码来定位对象。这似乎是代码中的瓶颈，在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

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

尺寸=500的结果：

编辑：如果峰值的数量很大(~100)，图像的大小很小，那么瓶颈实际上就是中心部分。所以，也许这部分的速度也需要优化。