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.
For me I had to change the single conversion to "im2single" on line 113 for the sub-pixel peak finding.
I use fastpeakfind with the followig parameters:
thres = 15; %value for binarisation, my image is uint coded and has maximum values of about 200
filt = (fspecial('gaussian', 200, 15)); %my peaks are nearly gaussiian with an amplitude of 200 and a sigma of about 15
edg=15; %I dont know exactly, but I thingk this is the distance to the edgo of the image zo exclude from search
res=1; %1: Local maximum 20%faster than 2: weighted centroid with subpixel interpolation
peaks=FastPeakFind(imagetest, thres, filt,edg,res); %imagetest is my image
This code works well and many thanks to Nathan.
I'm trying to use your function to locate a large set of peaks in some of my very messy data, but just to begin with, I can't get it running. I get the error
??? Undefined function or method 'fspecial' for input arguments of type 'char'.
Error in ==> FastPeakFind at 81
filt = (fspecial('gaussian', 7,1)); %if needed modify the filter according to the expected peaks sizes"
When looking up fspecial, I simply don't have this function installed, and I'm not sure where to get it from.
??? Undefined function or variable 'fspecial'.
I'm on MATLAB version R2011a, but I do have the Image Processing Toolbox, at least it seems so, since
>> license('test', 'image_toolbox')
gives 1 in return. On https://se.mathworks.com/help/images/ref/fspecial.html it also says that fspecial was "Introduced before R2006a".
My inputimage variable is a 600x600 double with values between 0 and 118325, I don't know if that changes anything? I'm confused as to why the error above was complaining about input arguments of type 'char' as well.
Thank you so much Natan, i got my answer,
What i mean by peak area is that if the data is flat at the peak (i mean for example when you have a peak surface instaed of just a peak point ),
Hi Mehrann, I dont understand what you mean by "peak area" vs "peak point". For the purpose of this code, a peak is defined by a point spread function (or PSF, see here https://en.wikipedia.org/wiki/Point_spread_function ) + some noise. The filter to use should be as close as possible to that theoretical PSF in order to optimize the peak detection. The reason is that the filter is then used to "blur" the data (using conv2). The output is are "blurred" peaks that have a unique local maximum value in each of them, and then the rest of the code is used to find that local maximum pixel in each peak. If that doesn't answer your question maybe send me an example so I'll be able to better understand what you ask.
Very nice job Natan, thank you
I just have a question? The 'filt' input parameter is a filter (for example a gaussian filter) which you have to define somewhere else and then use it in your code, right? And what it does is first implementing the filter and then find the peaks, right?
Because as i understood, your code does not work if there is peak area (instead of peak point) and you do not use filt, however if you define a simple filt then you will get at least one pixel per each peak area, am i correct?
excellent work, cleanly written and therefore easily adaptable to special needs
Hi Wim and thanks for the feedback! I'll fix the first point in the next update for sure. However, having a "flat" peak means for me that the either detector was saturated
( for example,the camera was over-exposed) , or that it is not really a peak but some distribution, because by definition a peak has a well defined maxima, or that the filter was not set properly. The latter point related to the 2d-convolution with some filter (the default is a 7x7 Gaussian, with sigma=1). if the data is still flat at the peak maybe a different filter matrix should be used. Please send me an example and we can discuss this further (my email is in the code header somewhere).
Nice piece of work, Natan!
However, I think there are 2 problems with your code:
1. if one inputs a floating point matrix with peaks < 1, the data gets scaled. You should scale the threshold then too
2. it cannot handle flat peaks: peaks that have more than 1 pixel of the same value in x or y direction will not be recognised
I am testing your code now with gray background but I get strange results. Discussion here http://stackoverflow.com/q/40119075/54964
This is an excellent piece of work! I think the sensitivity and specificity of the method should be included for each method. I get quite low sensitivities with the basic `function(I)`, analysing energy maps with rapid peaks. It would be great to extend it with statistical power!
This is a correction of my former wrong rate.
Excellent work, it would be nice to also have GPU support.
Hi Faruk, looking into your image1.tif file I see the issue of not recognizing the fainter peaks. The main reason for that, as I suspected, is that a flat threshold is not suitable for this image. In order to make it work you need to apply a local threshold, or do something that is called "common mode correction" of your image. Some cameras produce imaging data that should not directly be used in analysis and need corrections. Most popular corrections are: dark rate (pedestal) subtraction, bad pixel masking, common mode correction, gain correction, etc. Without going into too many details (you can search more info about the above), you can do the following "computational" trick, taking the local mean of the image you want to analyze. Denoting your image d:
ws=100; % window size
mean_d=imfilter(d,fspecial('average',ws),'replicate'); % this a local mean operation
Then apply the code "as usual" on dd, and set a flat threshold that will capture all the peaks.
Thanks Natan! I am sending you an image.
Hi Faruk, if I had to guess, your image may need a different type of threshold as it may be a non uniform background, another reason might be related to the PSF size of your peaks vs the filtering used, but without looking at your image it's hard for me to help you, so can you please send me an example of an image that does that and we can discuss this further? (my email is commented in the code)
I have a issue though. Some peaks look like are not being selected. If I lower down the threshold value, it gets lots more spots, zooming in indicates those are not peaks.
I was wondering if you have any suggestion to overcome this issue. Thanks.
Works well in different data sets and can overcome noisy data.
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.
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 ?
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
Great job! Very useful code!
If you includes a cluster size threshold to your code, it would be more than great!!! Anyway is great!
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.
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?
@Dan, I've added a sub-pixel resolution centroid feature to the code...
This is great! Makes imregionalmax seem really slow and unreliable. Would be nice to have a subpixel centroid as well.
please add the option to get a 2d bw mask of the peaks, in addition to the awkward coordinate vector ...
Thanks. it works well.
@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.
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
What is wrong?
Great algorithm and pretty fast. Thanks !!!!
@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.
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))"
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.
Added sub-pixel resolution using weighted centroids as a feature. Improved documentation.
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.
Added automatic uint8 support for improved performance. Added an output option of a binary matrix of peak positions besides the regular coordinate vector.
corrected typos from previous release that prevented the function from running for no arguments (demo mode)
Improved performance of saving to file and additional small improvements.
Added scaling correction for the case pixel values are all between 0 and 1. Improved file documentation.
minor editing, added functionality - when no input is used, the function generates random peaks data and plot a figure.
Code runs faster by casting to appropriate numeric classes for medfilt2 and conv2
bug fixed, code now handles images of arbitrary size.