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.
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.
@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
@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))"