Finding peaks in images
ImageAnalyst
peak would be defined as a region of higher gray level than the
surrounding area. For example, let's say you had topography of the
Rocky Mountains in Colorado. Well there are millions of local peaks.
Even little peaks on the sides of much larger peaks. Is there any way
to find the, say, three "largest" peaks. What algorithms are there?
I'm sure it can get quite complicated as to what exactly defines a
peak and what is too small to be a real peak and should just be
considered noise on a larger peak. But I'm open to ideas.
I've used imregionalmax() and imextendedmax() in MATLAB and they don't
really seem to do what I want, plus I don't really understand the
algorithm and they don't explain it. For example, they say:
"BW = imextendedmax(I,H) computes the extended-maxima transform, which
is the regional maxima of the H-maxima transform. H is a nonnegative
scalar.
Regional maxima are connected components of pixels with a constant
intensity value, and whose external boundary pixels all have a lower
value." Yep, clear as mud. Anyway it didn't work because it found
small dots in my largest peak and found a larger base on my smallest
peaks - very counter-intuitive.
Thanks,
ImageAnalyst
aruzinsky
Why don't you just mark pixels with luminance >= epsilon, starting
with epsilon = 255 and working downwards until you have three
segregated regions of whatever size you consider not noise?
Anyway, this paper is on a related subject but I doubt it will help:
http://www.ics.uci.edu/~irani/pubs/CameraReadyICIAP.pdf
Also, a key word is "hillocks".
Juho Pitkänen
I think that regional maximum just means that it can have several pixels
of the same value in the maximum. So largest peaks are probably more
sharp and have less pixels having same values. But Matlab is not my
tool, so just some simple ideas.
To get just locations of the largest peaks, you should get a list of the
locations and heights of the local maxima, sort it according to the
height value and select the largest one. Then you define a buffer around
the peak, perhaps making the buffer size proportional to the height of
the peak and find the next largest maximum that is not within the
buffer. And so on.
If you need the mountains as well, then one way is low pass filtering to
get rid of the little peaks on the sides, and watershed segmentation of
the negative filtered image to get segments around the local minima in
the negative image. It may be useful to select the size of the low pass
filter according to the height and filter more heavily on the high areas.
Juho
Bob Alvarez
to 2D. If you do, I'm interested :-)
http://www.mit.edu/~gari/CODE/ECGtools/ecgBag/localmax.m
function [m,i] = localmax(x, w)
% LOCALMAX Find indices and amplitudes of local maxima
% [m,i] = localmax(x, w) returns the indices and maxima defined
% over a local window of size 2w+1 given by w points on either
% side of the point being considered as a local maximum..
% inputs:
% x: vector (double) of data to find local max
% w: window
% outputs:
% m: data value at local maxima
% i: indices of local maxima
%
% by P.E.McSharry
% These routines are made available under the GNU general public
license.
% If you have not received a copy of this license, please download
from
% http://www.gnu.org/
%
% Please distribute (and modify) freely, commenting where you have
% added modifications. The author would appreciate correspondence
% regarding corrections, modifications, improvements etc.
%
% G. Clifford : ga...@ieee.org
N = length(x);
w = round(w); % change-- colon operator used below requires integer
limits
k = 2*w+1;
%-- not needed y = zeros(k,1);
m = [];
i = [];
l = 0;
for j=w+1:N-w
y = x(j-w:j+w);
[ymax,imax] = max(y);
if imax == w+1
l = l+1;
m(l) = ymax;
i(l) = j;
end
end
ImageAnalyst
>
> - Show quoted text -
-----------------------------------------------------------------------------------------------------
aruzinsky:
It's more difficult than that. Let's take a 1D example, let's say
your signal was this: [ 0 255 0 240 0 255 255 0 100 100 100 100 100
100 0]
The method you suggested would find the first 3 peaks first, but the
third peak would have less volume (510) than the 4th peak (600). So
you can't quit when you have your required number. Now let's take
another case: [0 255 30 100 100 220 100 200 100 200 100 0]
Going down from 255, you'd first have 1 peak, then you'd have 2, then
you'd have 4, and then (at less than 100) you'd be back down to only 2
peaks again, and finally (less than 30) you'd have only 1 peak. So it
can get pretty tricky, particularly when you have peaks on top of
peaks. Now I know there are lots of 1D peak detector algorithms
developed by the mass-spectrometry and chromatography community, but
I'm looking for a 2D one that makes sense.
I'm currently looking into morphological methods such as:
L. Vincent. (1993). Morphological Grayscale Reconstruction in Image
Analysis: Applications and Efficient Algorithms. IEEE Transactions on
Image Processing, vol.2 (2) 176-201
and
http://www.afscet.asso.fr/resSystemica/Crete02/Halkiotis.pdf
It seems like they have potential but are a bit tricky to understand.
Doing a watershed and then calculating volume in each zone (region)
may also work as long as I can get a decent separation of peaks and
not have it over-segmented.
Thanks,
ImageAnalyst
ImageAnalyst
>
> - Show quoted text -
------------------------------------------------------------
Juho:
Thanks for your reply. I've seen some papers that have had some
success via the watershed method and I might try that.
ImageAnalyst
ImageAnalyst
> Here is a 1D localmax function that works well. Maybe you can extend
>
> function [m,i] = localmax(x, w)
> % LOCALMAX Find indices and amplitudes of local maxima
> % [m,i] = localmax(x, w) returns the indices and maxima defined
> % over a local window of size 2w+1 given by w points on either
> % side of the point being considered as a local maximum..
> % inputs:
> % x: vector (double) of data to find local max
> % w: window
> % outputs:
> % m: data value at local maxima
> % i: indices of local maxima
> %
> % by P.E.McSharry
> % These routines are made available under the GNU general public
> license.
> % If you have not received a copy of this license, please download
> from
> %
> % Please distribute (and modify) freely, commenting where you have
> % added modifications. The author would appreciate correspondence
> % regarding corrections, modifications, improvements etc.
> %
>
> N = length(x);
> w = round(w); % change-- colon operator used below requires integer
> limits
> k = 2*w+1;
> %-- not needed y = zeros(k,1);
> m = [];
> i = [];
> l = 0;
> for j=w+1:N-w
> y = x(j-w:j+w);
> [ymax,imax] = max(y);
> if imax == w+1
> l = l+1;
> m(l) = ymax;
> i(l) = j;
> end
> end
---------------------------------------------------------------------------------------------
Bob:
Thanks for the link. I did try that code but for noisy data, it did a
terrible job. Essentially finding every other point in the whole 1D
array - since about half the values in a noisy array are local max and
half are local mins. That algorithm didn't do as good a job as the
one I'm currently using for 1D peak detection:
http://billauer.co.il/peakdet.html
There are lots and lots of 1D peak detection algorithm:
http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2631518
reviews several of them. It appears most of them come from the
chemistry community where they have a variety of instruments that
generate 1D signals (mass sepctrometry, etc.)
2D peak detection is a lot more complicated, particularly if you don't
just want to find the location of the peak's very tip top but want to
find the region that its base covers. I might take a look deeper into
watershed, such as
http://cmm.ensmp.fr/~beucher/wtshed.html
or else morphological methods like promoted by Dr. Luc Vincent
(currently at Google). L. Vincent. (1993). Morphological Grayscale
Reconstruction in Image Analysis: Applications and Efficient
Algorithms. IEEE Transactions on Image Processing, vol.2 (2) 176-201
Thanks,
ImageAnalyst
aruzinsky
> aruzinsky:
> It's more difficult than that. Let's take a 1D example, let's say
> your signal was this: [ 0 255 0 240 0 255 255 0 100 100 100 100 100
> 100 0]
> The method you suggested would find the first 3 peaks first, but the
> third peak would have less volume (510) than the 4th peak (600). So
> you can't quit when you have your required number. Now let's take
> another case: [0 255 30 100 100 220 100 200 100 200 100 0]
> Going down from 255, you'd first have 1 peak, then you'd have 2, then
> you'd have 4, and then (at less than 100) you'd be back down to only 2
> peaks again, and finally (less than 30) you'd have only 1 peak. So it
> can get pretty tricky, particularly when you have peaks on top of
> peaks. Now I know there are lots of 1D peak detector algorithms
> developed by the mass-spectrometry and chromatography community, but
> I'm looking for a 2D one that makes sense.
>
> I'm currently looking into morphological methods such as:
> L. Vincent. (1993). Morphological Grayscale Reconstruction in Image
> Analysis: Applications and Efficient Algorithms. IEEE Transactions on
> Image Processing, vol.2 (2) 176-201
> It seems like they have potential but are a bit tricky to understand.
> Doing a watershed and then calculating volume in each zone (region)
> may also work as long as I can get a decent separation of peaks and
> not have it over-segmented.
>
> Thanks,
>
> - Show quoted text -
Regarding, "Going down from 255, you'd first have 1 peak, then you'd
have 2, then you'd have 4, and then (at less than 100) you'd be back
down to only 2
peaks again, and finally (less than 30) you'd have only 1 peak.", you
didn't put a minimum size restriction on the regions as I suggested.
Putting a size restriction >=2, and working down, the first peak of
[0 255 30 100 100 220 100 200 100 200 100 0]
is
[100 100 220 100 200 100 200 100]
and the second peak is
[255 30].
Also, noise can be reduced by general methods and then region size
disregarded. Suppose we apply a median filter (which would be
appropriate for shot noise) of length three to your examples,
reflecting at borders, and disregard region size.
For
[ 0 255 0 240 0 255 255 0 100 100 100 100 100 100 0]
->
[ 0 0 240 0 240 255 255 100 100 100 100 100 100 0 ]
Now, you only have 2 peaks.
For
[0 255 30 100 100 220 100 200 100 200 100 0]
->
[0 30 100 100 100 100 200 100 200 100 100 0]
Again, you only have 2 peaks.