Pencarian puncak dan pemusatan massa yang cepat dengan python

Saya mencoba mengembangkan algoritma cepat dengan python untuk menemukan puncak dalam suatu gambar dan kemudian menemukan pusat massa dari puncak tersebut. Saya telah menulis kode berikut menggunakan scipy.ndimage.label dan ndimage.find_objects untuk menemukan lokasi objek. Ini tampaknya menjadi hambatan dalam kode, dan dibutuhkan sekitar 7 ms untuk menemukan 20 objek dalam gambar 500x500. Saya ingin memperbesarnya menjadi gambar yang lebih besar (2000x2000), tetapi kemudian waktunya bertambah hingga hampir 100 ms. Jadi, saya ingin tahu apakah ada opsi yang lebih cepat.

Berikut adalah kode yang saya miliki sejauh ini, yang berfungsi, namun lambat. Pertama saya mensimulasikan data saya menggunakan beberapa puncak gaussian. Bagian ini lambat, namun dalam praktiknya saya akan menggunakan data nyata, jadi saya tidak terlalu peduli untuk mempercepat bagian itu. Saya ingin dapat menemukan puncaknya dengan sangat cepat.

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

Hasil untuk ukuran=500: masukkan deskripsi gambar di sini

EDIT: Jika jumlah puncaknya besar (~100) dan ukuran gambarnya kecil, maka kemacetan sebenarnya adalah bagian centroiding. Nah, mungkin kecepatan bagian ini juga perlu dioptimalkan.


person DanHickstein    schedule 01.10.2013    source sumber
comment
Lihat github.com/tacaswell/trackpy dan github.com/nkeim/trackpy Mungkin menarik. (penafian: salah satunya adalah kode saya dan yang lainnya adalah percabangan kode saya oleh mantan teman lab)   -  person tacaswell    schedule 01.10.2013
comment
Oh, sepertinya menarik sekali! Saya mungkin harus memeriksa fungsi identifikasi.find_local_max dan identifikasi.subpixel_centroid, bukan?   -  person DanHickstein    schedule 02.10.2013


Jawaban (4)


Metode Anda untuk menemukan puncak (pengambangan sederhana) tentu saja sangat sensitif terhadap pilihan ambang batas: atur terlalu rendah dan Anda akan "mendeteksi" hal-hal yang bukan puncak; atur terlalu tinggi dan Anda akan kehilangan puncak yang valid.

Ada alternatif yang lebih kuat, yang akan mendeteksi semua nilai maksimum lokal dalam intensitas gambar, berapa pun nilai intensitasnya. Pilihan saya adalah menerapkan dilatasi dengan elemen penataan kecil (5x5 atau 7x7), lalu mencari piksel di mana gambar asli dan versi yang diperbesar memiliki nilai yang sama. Ini berfungsi karena, menurut definisi, dilatasi(x, y, E, img) = { maksimal img dalam E berpusat pada piksel (x,y) }, dan oleh karena itu dilatasi(x, y, E, img) = img(x , y) kapanpun (x,y) adalah lokasi maksimum lokal pada skala E.

Dengan implementasi cepat dari operator morfologi (misalnya yang ada di OpenCV), algoritme ini memiliki ukuran gambar yang linier dalam ruang dan waktu (satu buffer berukuran gambar tambahan untuk gambar yang diperbesar, dan satu meneruskan keduanya). Dalam keadaan darurat, ini juga dapat diimplementasikan secara online tanpa buffering tambahan dan sedikit kerumitan tambahan, dan masih dalam waktu linier.

Untuk lebih memperkuatnya di hadapan garam-dan-merica atau gangguan serupa, yang dapat menimbulkan banyak maxima palsu, Anda dapat menerapkan metode ini dua kali, dengan elemen penataan dengan ukuran berbeda (misalnya, 5x5 dan 7x7), lalu hanya mempertahankan elemen stabilnya. maxima, di mana stabilitas dapat ditentukan oleh posisi maxima yang tidak berubah, atau dengan posisi yang tidak berubah lebih dari satu piksel, dll. Selain itu, Anda mungkin ingin menekan maxima rendah di dekatnya jika Anda memiliki alasan untuk meyakini bahwa maxima tersebut disebabkan oleh kebisingan. Cara yang efisien untuk melakukan hal ini adalah pertama-tama mendeteksi semua maxima lokal seperti di atas, mengurutkannya berdasarkan ketinggian, lalu turunkan daftar yang diurutkan dan menyimpannya jika nilainya dalam gambar tidak berubah dan, jika dipertahankan, setel ke nolkan semua piksel dalam lingkungan (2d+1) x (2d+1), dengan d adalah jarak minimum antara maksimum terdekat yang ingin Anda toleransi.

person Francesco Callari    schedule 02.10.2013
comment
Terima kasih untuk sarannya! Menurut saya algoritme yang Anda sarankan cukup mirip dengan fungsi find_local_maxima yang disarankan oleh @tcaswell dalam kode trackpy-nya: github.com/tacaswell/trackpy/blob/master/trackpy/ Pendekatan ini bekerja dengan sangat baik, dan memang menemukan puncak yang sebenarnya, meskipun sebagiannya tumpang tindih, atau jika beberapa puncaknya tumpang tindih. jauh lebih kecil. Namun, kodenya jauh lebih lambat, memerlukan ~200 ms untuk gambar berukuran 500x500. Saya belum pernah menggunakan OpenCV sebelumnya, tapi saya akan mencobanya. Bisakah Anda menyarankan secara lebih spesifik fungsi OpenCV mana yang harus saya gunakan? - person DanHickstein; 02.10.2013
comment
Pelebaran OpenCV: docs.opencv.org/modules/imgproc /doc/ (atau docs.opencv yang lebih cepat .org/modules/ocl/doc/image_filtering.html#ocl-dilate jika Anda dapat menggunakan OCL), menggunakan docs.opencv.org/modules/imgproc/doc/ jika Anda merasa malas. Untuk perbandingan cepat dengan penggunaan asli docs.opencv.org/modules/core /dok/ - person Francesco Callari; 02.10.2013

Jika Anda memiliki banyak puncak, lebih cepat menggunakan scipy.ndimage.center_of_mass. Anda dapat mengganti kode Anda mulai dari definisi peak_slices, hingga total waktu pencetakan, dengan dua baris berikut:

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]

Untuk num_peaks = 20 ini berjalan sekitar 3x lebih lambat dibandingkan pendekatan Anda, namun untuk num_peaks = 100 ini berjalan sekitar 10x lebih cepat. Jadi pilihan terbaik Anda akan bergantung pada data aktual Anda.

person Jaime    schedule 01.10.2013
comment
Ah, kode yang lebih bersih dan lebih cepat untuk banyak puncak, sempurna! Terima kasih! Namun saya juga sangat menginginkan cara untuk mempercepat proses pelabelan gambar dan pencarian objek. - person DanHickstein; 02.10.2013
comment
Saya benar-benar ragu ada cara yang lebih cepat daripada fungsi ndimage yang dapat Anda kode dengan Python... - person Jaime; 02.10.2013

Pendekatan lainnya adalah menghindari semua sum(), meshgrid() dan lainnya. Ganti semuanya dengan aljabar linier lurus.

>>> def centroid2(data):
    h,w=data.shape
    x=np.arange(h)
    y=np.arange(w)
    x1=np.ones((1,h))
    y1=np.ones((w,1))
    return ((np.dot(np.dot(x1, data), y))/(np.dot(np.dot(x1, data), y1)),
            (np.dot(np.dot(x, data), y1))/(np.dot(np.dot(x1, data), y1)))
#be careful, it returns two arrays

Hal ini juga dapat diperluas ke dimensi yang lebih tinggi. Kecepatan 60% dibandingkan dengan centroid()

person CT Zhu    schedule 01.10.2013

Penghitungan centroid berikut ini lebih cepat dibandingkan keduanya, terutama untuk data berukuran besar:

def centroidnp(data):
    h,w = data.shape
    x = np.arange(w)
    y = np.arange(h)
    vx = data.sum(axis=0)
    vx /= vx.sum()
    vy = data.sum(axis=1)
    vy /= vy.sum()    
    return np.dot(vx,x),np.dot(vy,y)
person CyxAndr    schedule 28.01.2018
comment
ada komentar mengapa ini lebih cepat? - person Palu Macil; 28.01.2018