Databases
StackOverflow2013
Postsdbo
POFast peak-finding and centroiding in python
 

Note that there are some explanatory texts on larger screens.

plurals
  1. POFast peak-finding and centroiding in python
    primarykey
    data
    text
    <p>I am trying to develop a fast algorithm in python for finding peaks in an image and then finding the centroid of those peaks. I have written the following code using the scipy.ndimage.label and ndimage.find_objects for locating the objects. This seems to be the bottleneck in the code, and it takes about 7 ms to locate 20 objects in a 500x500 image. I would like to scale this up to larger (2000x2000) image, but then the time increases to almost 100 ms. So, I'm wondering if there is a faster option.</p> <p>Here is the code that I have so far, which works, but is slow. First I simulate my data using some gaussian peaks. This part is slow, but in practice I will be using real data, so I don't care too much about speeding that part up. I would like to be able to find the peaks very quickly.</p> <pre><code>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&lt;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() </code></pre> <p>The results for size=500: <img src="https://i.stack.imgur.com/CNONo.png" alt="enter image description here"></p> <p>EDIT: If the number of peaks is large (~100) and the size of the image is small, then the bottleneck is actually the centroiding part. So, perhaps the speed of this part also needs to be optimized.</p>
    singulars
    1. PO
      singulars
      1. PTAnswer
    1. This table or related slice is empty.
    1. PTQuestion
    1. USDanHickstein
    1. USDanHickstein
    plurals
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. PO
      singulars
      1. PTAnswer
    2. PO
      singulars
      1. PTAnswer
    3. PO
      singulars
      1. PTAnswer
    1. VO
      singulars
      1. POFast peak-finding and centroiding in python
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. POFast peak-finding and centroiding in python
      1. USunutbu
      1. VTFavorite
    3. VO
      singulars
      1. POFast peak-finding and centroiding in python
      1. This table or related slice is empty.
      1. VTUpMod
    1. COSee https://github.com/tacaswell/trackpy and https://github.com/nkeim/trackpy Might be of interest. (disclaimer: one of this is my code and the other is a fork of my code by a former lab-mate)
      singulars
      1. POFast peak-finding and centroiding in python
      1. UStacaswell
    2. COOh, looks very interesting! I should probably check out the identification.find_local_max and identification.subpixel_centroid functions, right?
      singulars
      1. POFast peak-finding and centroiding in python
      1. USDanHickstein
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload