การค้นหาจุดสูงสุดและเซนทรอยด์อย่างรวดเร็วใน python

ฉันกำลังพยายามพัฒนาอัลกอริธึมที่รวดเร็วใน python เพื่อค้นหาพีคในรูปภาพ จากนั้นค้นหาเซนทรอยด์ของพีคเหล่านั้น ฉันได้เขียนโค้ดต่อไปนี้โดยใช้ scipy.ndimage.label และ ndimage.find_objects เพื่อค้นหาวัตถุ นี่ดูเหมือนจะเป็นจุดคอขวดในโค้ด และใช้เวลาประมาณ 7 มิลลิวินาทีในการค้นหาวัตถุ 20 ชิ้นในรูปภาพขนาด 500x500 ฉันต้องการปรับขนาดภาพให้ใหญ่ขึ้น (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()

ผลลัพธ์สำหรับขนาด=500: ป้อนคำอธิบายรูปภาพที่นี่

แก้ไข: หากจำนวนพีคมีขนาดใหญ่ (~ 100) และขนาดของรูปภาพเล็ก แสดงว่าคอขวดคือส่วนที่อยู่ตรงกลาง ดังนั้นบางทีความเร็วของส่วนนี้ก็จำเป็นต้องได้รับการปรับให้เหมาะสมเช่นกัน


person DanHickstein    schedule 01.10.2013    source แหล่งที่มา
comment
ดู github.com/tacaswell/trackpy และ github.com/nkeim/trackpy อาจเป็นที่สนใจ (ข้อจำกัดความรับผิดชอบ: หนึ่งในนั้นคือรหัสของฉันและอีกอันคือทางแยกของรหัสของฉันโดยอดีตเพื่อนร่วมห้องทดลอง)   -  person tacaswell    schedule 01.10.2013
comment
โอ้ดูน่าสนใจมาก! ฉันน่าจะตรวจสอบฟังก์ชันการระบุตัวตน.find_local_max และการระบุตัวตน.subpixel_centroid ใช่ไหม   -  person DanHickstein    schedule 02.10.2013


คำตอบ (4)


วิธีการหาจุดสูงสุดของคุณ (การกำหนดเกณฑ์ขั้นต่ำแบบง่าย) แน่นอนว่ามีความละเอียดอ่อนมากต่อการเลือกเกณฑ์: หากตั้งค่าต่ำเกินไปแล้วคุณจะ "ตรวจจับ" สิ่งต่างๆ ที่ไม่ใช่จุดสูงสุดได้ หากตั้งค่าสูงเกินไป คุณจะพลาดยอดเขาที่ถูกต้อง

มีทางเลือกอื่นที่มีประสิทธิภาพมากกว่า ซึ่งจะตรวจจับค่าสูงสุดเฉพาะจุดทั้งหมดในความเข้มของภาพ โดยไม่คำนึงถึงค่าความเข้มของแสงเหล่านั้น สิ่งที่ฉันชอบคือการใช้การขยายด้วยองค์ประกอบโครงสร้างขนาดเล็ก (5x5 หรือ 7x7) จากนั้นค้นหาพิกเซลที่ภาพต้นฉบับและเวอร์ชันที่ขยายมีค่าเท่ากัน วิธีนี้ได้ผลเพราะตามคำจำกัดความ การขยาย (x, y, E, img) = { สูงสุดของ img ภายใน E โดยมีศูนย์กลางที่พิกเซล (x, y) } และด้วยเหตุนี้ การขยาย (x, y, E, img) = img(x , y) เมื่อใดก็ตามที่ (x,y) คือตำแหน่งของค่าสูงสุดเฉพาะที่ในระดับ E

ด้วยการใช้งานตัวดำเนินการทางสัณฐานวิทยาอย่างรวดเร็ว (เช่น ตัวดำเนินการใน OpenCV) อัลกอริธึมนี้จะมีขนาดของรูปภาพเป็นเส้นตรงทั้งในพื้นที่และเวลา (บัฟเฟอร์ขนาดรูปภาพเพิ่มเติมหนึ่งบัฟเฟอร์สำหรับรูปภาพที่ถูกขยาย และอีกหนึ่งบัฟเฟอร์ส่งผ่านทั้งคู่) ในเวลาสั้นๆ ก็สามารถนำไปใช้งานออนไลน์ได้โดยไม่ต้องมีบัฟเฟอร์เพิ่มเติมและมีความซับซ้อนเพิ่มขึ้นอีกเล็กน้อย และยังคงเป็นเวลาเชิงเส้น

หากต้องการเสริมความแข็งแกร่งให้มากขึ้นเมื่อมีเกลือและพริกไทยหรือเสียงรบกวนที่คล้ายกัน ซึ่งอาจก่อให้เกิดค่าสูงสุดที่ผิดพลาดจำนวนมาก คุณสามารถใช้วิธีนี้ได้สองครั้ง โดยมีองค์ประกอบโครงสร้างที่มีขนาดแตกต่างกัน (เช่น 5x5 และ 7x7) จากนั้นคงไว้เพียงค่าคงที่ สูงสุด โดยที่ความเสถียรสามารถกำหนดได้โดยการเปลี่ยนตำแหน่งของจุดสูงสุด หรือโดยตำแหน่งที่ไม่เปลี่ยนแปลงมากกว่าหนึ่งพิกเซล เป็นต้น นอกจากนี้ คุณอาจต้องการลดค่าสูงสุดใกล้เคียงที่ต่ำเมื่อคุณมีเหตุผลที่เชื่อได้ว่าเกิดจากสัญญาณรบกวน วิธีที่มีประสิทธิภาพในการทำเช่นนี้คือการตรวจจับค่าสูงสุดเฉพาะที่ด้านบนทั้งหมดก่อน จากนั้นเรียงลำดับจากมากไปน้อยตามความสูง จากนั้นลงไปตามรายการที่เรียงลำดับแล้วเก็บไว้หากค่าในภาพไม่เปลี่ยนแปลง และหากคงไว้ ให้ตั้งค่าเป็น พิกเซลทั้งหมดเป็นศูนย์ในย่านใกล้เคียง (2d+1) x (2d+1) โดยที่ d คือระยะห่างขั้นต่ำระหว่างค่าสูงสุดใกล้เคียงที่คุณยินดียอมรับ

person Francesco Callari    schedule 02.10.2013
comment
ขอบคุณสำหรับข้อเสนอแนะ! ฉันคิดว่าอัลกอริทึมที่คุณแนะนำค่อนข้างคล้ายกับฟังก์ชัน find_local_maxima ที่แนะนำโดย @tcaswell ในโค้ด trackpy ของเขา: github.com/tacaswell/trackpy/blob/master/trackpy/ วิธีการนี้ใช้ได้ผลดีมาก และค้นหาจุดสูงสุดที่แท้จริงได้ แม้ว่าจะทับซ้อนกันบางส่วน หรือหากบางจุดนั้นทับซ้อนกัน เล็กกว่ามาก อย่างไรก็ตาม โค้ดช้ากว่ามาก โดยต้องใช้ ~200 ms สำหรับรูปภาพขนาด 500x500 ฉันยังไม่เคยใช้ OpenCV มาก่อน แต่จะลองดู คุณช่วยแนะนำฟังก์ชัน OpenCV ที่ฉันควรใช้ให้เจาะจงมากขึ้นได้ไหม - person DanHickstein; 02.10.2013
comment
การขยาย OpenCV: docs.opencv.org/modules/imgproc /doc/ (หรือ docs.opencv ที่เร็วกว่า .org/modules/ocl/doc/image_filtering.html#ocl-dilate หากคุณสามารถใช้ OCL) โดยใช้ docs.opencv.org/modules/imgproc/doc/ หากคุณรู้สึกขี้เกียจ เพื่อการเปรียบเทียบที่รวดเร็วกับการใช้งานดั้งเดิม docs.opencv.org/modules/core /doc/ - person Francesco Callari; 02.10.2013

หากคุณมีพีคหลายจุด การใช้ scipy.ndimage.center_of_mass จะเร็วกว่า คุณสามารถแทนที่โค้ดของคุณโดยเริ่มต้นด้วยคำจำกัดความ peak_slices จนถึงการพิมพ์เวลาทั้งหมด โดยมีสองบรรทัดต่อไปนี้:

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 เท่า ดังนั้นตัวเลือกที่ดีที่สุดของคุณจะขึ้นอยู่กับข้อมูลจริงของคุณ

person Jaime    schedule 01.10.2013
comment
อ่า โค้ดสะอาดกว่าและเร็วกว่าหลาย ๆ จุดสูงสุด สมบูรณ์แบบ! ขอบคุณ! แต่ฉันก็อยากได้วิธีเพิ่มความเร็วในการติดป้ายกำกับรูปภาพและกระบวนการค้นหาวัตถุด้วย - person DanHickstein; 02.10.2013
comment
ฉันสงสัยจริงๆ ว่าจะมีวิธีที่เร็วกว่าฟังก์ชัน ndimage ที่คุณสามารถเขียนโค้ดใน Python ได้... - person Jaime; 02.10.2013

อีกวิธีหนึ่งคือการหลีกเลี่ยง sum(), meshgrid() และสิ่งต่างๆ ทั้งหมด แทนที่ทุกอย่างด้วยพีชคณิตเชิงเส้นตรง

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

สามารถนำไปใช้ในมิติที่สูงขึ้นได้เช่นกัน 60% ของการเร่งความเร็วเทียบกับ centroid()

person CT Zhu    schedule 01.10.2013

การคำนวณเซนทรอยด์ต่อไปนี้เร็วกว่าทั้งสองอย่าง โดยเฉพาะสำหรับข้อมูลขนาดใหญ่:

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
ความคิดเห็นใด ๆ ว่าทำไมมันถึงเร็วขึ้น? - person Palu Macil; 28.01.2018