10

画像内のピークを見つけて、それらのピークの重心を見つけるための高速アルゴリズムを 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()

size=500 の結果: ここに画像の説明を入力

編集: ピークの数が多く (~100)、画像のサイズが小さい場合、ボトルネックは実際には重心部分です。したがって、おそらくこの部分の速度も最適化する必要があります。

4

4 に答える 4

9

ピークを見つけるための方法 (単純なしきい値処理) は、もちろん、しきい値の選択に非常に敏感です。設定が低すぎると、ピークではないものを「検出」します。設定が高すぎると、有効なピークを見逃すことになります。

強度値に関係なく、画像強度のすべての極大値を検出する、より堅牢な代替手段があります。私の好みは、小さな (5x5 または 7x7) 構造要素で膨張を適用し、元の画像とその膨張バージョンが同じ値を持つピクセルを見つけることです。これは、定義により、 dilation(x, y, E, img) = { ピクセル (x,y) を中心とする E 内の img の最大値 }、したがって dilation(x, y, E, img) = img(x , y) (x,y) が E のスケールでの極大値の位置であるときはいつでも。

形態学的演算子 (OpenCV のものなど) の高速な実装により、このアルゴリズムは、空間と時間の両方で画像のサイズに線形になります (拡張画像用に 1 つの余分な画像サイズのバッファーと、両方に 1 つのパス)。ピンチの場合は、追加のバッファーや少し余分な複雑さなしにオンラインで実装することもできますが、それでも線形時間です。

多くの偽の最大値を導入する可能性のあるごま塩または同様のノイズの存在下でそれをさらに堅牢化するには、異なるサイズの構造化要素 (たとえば、5x5 と 7x7) を使用してメソッドを 2 回適用し、安定したものだけを保持します。ここで、安定性は、最大値の位置が変化しないこと、または位置が 1 ピクセル以上変化しないことなどによって定義できます。さらに、ノイズが原因であると信じる理由がある場合は、近くの低い最大値を抑制したい場合があります。これを行う効率的な方法は、最初に上記のようにすべての極大値を検出し、それらを高さの降順に並べ替えてから、並べ替えられたリストを下に移動し、画像内の値が変更されていない場合はそれらを保持し、保持されている場合はに設定することです。それらの (2d+1) x (2d+1) 近傍のすべてのピクセルをゼロにします。ここで、d は許容できる近くの最大値間の最小距離です。

于 2013-10-02T12:49:40.390 に答える
5

ピークが多い場合は を使った方が早いscipy.ndimage.center_of_massです。の定義からpeak_slices合計時間の出力までのコードを次の 2 行に置き換えることができます。

centroids = scipy.ndimage.center_of_mass(Z_thresh, labeled_image,
                                         np.arange(1, number_of_objects + 1))
centroids = [(j, i) for i, j in centroids]

これは、アプローチよりもnum_peaks = 20約 3 倍遅くnum_peaks = 100実行されますが、約 10倍速く実行されます。したがって、最適なオプションは実際のデータによって異なります。

于 2013-10-01T19:03:30.240 に答える