Code covered by the BSD License  

Highlights from
Fast 2D peak finder

5.0
5.0 | 9 ratings Rate this file 236 Downloads (last 30 days) File Size: 4.14 KB File ID: #37388 Version: 1.12
image thumbnail

Fast 2D peak finder

by

Natan (view profile)

 

03 Jul 2012 (Updated )

Find local maxima \ peak positions in noisy 2D arrays

| Watch this File

File Information
Description

A simple and fast 2D peak finder. The aim was to be faster than more sophisticated techniques yet good enough to find peaks in noisy data. The code analyzes noisy 2D images and find peaks using robust local maxima finder (1 pixel resolution) or by weighted centroids (sub-pixel resolution). The code is designed to be as fast as possible, so I kept it pretty basic. It best works when using uint16 \ uint8 images, and assumes that peaks are relatively sparse.

The code requires Matlab's Image Processing Toolbox, and can be used inside parfor for faster processing times.

Required Products Image Processing Toolbox
MATLAB release MATLAB 8.2 (R2013b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (20)
07 May 2015 Gael Leveque  
05 May 2015 Natan

Natan (view profile)

Hi Gael, I'd like to help you, however, without the images that you're talking about and the exact way you used the code, I can't see how I can do that. Please email me (my email in commented in the code itself) with more info and I'll try to help you solve the problem if I can.

Comment only
05 May 2015 Gael Leveque

Great job!
I have some issues for a few pictures though. When the noisy peak I want to detect is too sharp, if I use the first mode too many peaks are found (I can't increase too much the matrix filter value or the real peak is flattened) and if I use the second mode, the function returns roughly the middle of the picture...
Any idea of how I could get trough this ?

Comment only
23 Jan 2015 Natan

Natan (view profile)

Bruno, you can edit the code to include the cluster size threshold very easily. Just write: bw = bwareaopen(logical(d), thresC);
after line 116. where tresC is the cluster threshold you want. Then edit the line that reads "[x y]=find(..." to [x y]=find(bw(edg:sd(1)...

Read more about bwareaopen here: http://www.mathworks.com/help/images/ref/bwareaopen.html

Comment only
21 Jan 2015 Brunno Machado de Campos

Great job! Very useful code!
If you includes a cluster size threshold to your code, it would be more than great!!! Anyway is great!

19 Aug 2014 Yuanwei  
10 May 2014 Natan

Natan (view profile)

The default filter matrix is a 7x7 Gaussian, and the smoothing is done via 2D convolution of the image with that filter. If you don't understand why this smooths the data you can refer to this simple demo (http://blogs.mathworks.com/videos/2012/04/17/using-convolution-to-smooth-data-with-a-moving-average-in-matlab/). The size of the filter was selected for typical PSF's, and what it does is to take the weighted average for each pixel with it's 7x7 pixels surrounding it (so it's a moving average), where the weights are normally distributed (hence a Gaussian). If you have more questions email me.

Comment only
09 May 2014 Jochen

Jochen (view profile)

Thanks a lot for that brilliant function!
But I don't understand how the filter matrix is supposed to look like to smooth the image. Can you give me hint?
Thanks!

18 Apr 2014 Andrey

Andrey (view profile)

 
11 Oct 2013 Natan

Natan (view profile)

@Dan, I've added a sub-pixel resolution centroid feature to the code...

Comment only
09 Oct 2013 Dan Hickstein

This is great! Makes imregionalmax seem really slow and unreliable. Would be nice to have a subpixel centroid as well.

08 May 2013 Roman Voronov

please add the option to get a 2d bw mask of the peaks, in addition to the awkward coordinate vector ...

02 Apr 2013 Sujay

Sujay (view profile)

Thanks. it works well.

01 Apr 2013 Natan

Natan (view profile)

@Sujay, as the code documentation say in line 8: "The 2D data raw image - assumes a Double\Single-precision ... or unit16 array. Please note that the code casts the raw image to uint16, if your image is 8-bit depth, I recommend to change all the places that uint16 is used to uint8 for faster run times." So, either use FastPeakFind(uint16(image),...) or change all uint16 to uint8 in the code. In future updates I'll include this feature automatically.

Comment only
27 Mar 2013 Sujay

Sujay (view profile)

I have a 512x512 .tif image that I can read in matlab as, A=imread(''). I am trying to use this program as,
FastPeakFind(A). It is showing the following error:
??? Error using ==> times
Integers can only be combined with integers of the same
class, or scalar doubles.

Error in ==> FastPeakFind at 81
d=d.*uint16(d>threshold);

What is wrong?

Comment only
23 Oct 2012 Slawomir

Great algorithm and pretty fast. Thanks !!!!

Comment only
21 Aug 2012 Natan

Natan (view profile)

@Carl Witthoft, Thanks for the comments, I've updated the information regarding the image processing toolbox requirement the scaling in case of a 0:1 range of image values. I've also already tested some more compact local maxima conditions, such as:
if all ( reshape( d(x(j),y(j)) >= d(x(j)-1:x(j)+1,y(j)-1:y(j)+1),9,1))

if d(x(j),y(j)) == max(max(d((x(j)-1):(x(j)+1),(y(j)-1):(y(j)+1))))

if d(x(j),y(j)) == max( reshape( d(x(j),y(j)) >= d(x(j)-1:x(j)+1,y(j)-1:y(j)+1),9,1))

All these were a factor of 2-3 slower than the simple 8-fold if condition. It takes Matlab more time to do the max or reshape functions than to go through the simple conditions...

Additionally, Matlab's image processing toolbox has the function imregionalmax, but I've found it is 6-7 times slower than this code.

Comment only
13 Aug 2012 Carl Witthoft

First commment lost.. try again. 1) Please specify that ImageProcessingToolbox is required. 2) For speed reasons, it may help to replace the 8-fold logical test with something like "d(j,j) > max(max(d(j-1:j+1,j-1:j+1))"

Comment only
13 Aug 2012 Carl Witthoft

One more thought: to allow for input files with float ranges (e.g. 0:1.0), replace the uint16 cast with "d = uint16( d.*2^16./(max(max(d))" or similar scaling function.

Comment only
25 Jul 2012 Song

Song (view profile)

 
Updates
09 Jul 2012 1.1

bug fixed, code now handles images of arbitrary size.

16 Jul 2012 1.2

Code runs faster by casting to appropriate numeric classes for medfilt2 and conv2

10 Aug 2012 1.3

minor editing, added functionality - when no input is used, the function generates random peaks data and plot a figure.

21 Aug 2012 1.4

Added scaling correction for the case pixel values are all between 0 and 1. Improved file documentation.

05 Feb 2013 1.6

Improved performance of saving to file and additional small improvements.

11 Feb 2013 1.7

corrected typos from previous release that prevented the function from running for no arguments (demo mode)

05 Mar 2013 1.8

bug fixes

09 May 2013 1.9

Added automatic uint8 support for improved performance. Added an output option of a binary matrix of peak positions besides the regular coordinate vector.

06 Jun 2013 1.10

the case binary matrix output when the image is all zeros returned an error (bug now fixed). Many thanks to Roman Voronov for spotting this.

11 Oct 2013 1.12

Added sub-pixel resolution using weighted centroids as a feature. Improved documentation.

Contact us