I have an intensity image, similar to the image below. I would like to very quickly, segment this image, into separate peaks. Where each peak is labelled, and segmented. A notional segmentation of 3 blobs is given below by the black lines.

Definition of a peak segment is: The contour that just crosses the highest saddle point between this peak and an adjacent peak.

This has to be a well studied problem, but I can't figure it out. I have been looking at the Matlab functions:

all of which are close.... It is not hard to write an algorithm to do this (see below) using imreconstruct to operate on the tallest peak, segment out the first peak, and repeat. The algorithm below is O(N) in number of peaks. Such an algorithm can and should be O(1) in number of peaks (and this is important for my application.)

enter image description here

The above image was rendered in matlab with the command


Here is an algorithm that performs the segmentation:

a = peaks;
peakValsX = [21, 25, 35];
peakValsY = [20, 38, 25];

N = numel(peakValsX);

clear('peakImgs', 'peakMask');
for ix = 1 : N;
    aptImg       = zeros(size(a));
    aptImg(peakValsY(ix), peakValsX(ix)) = a(peakValsY(ix),peakValsX(ix));
    peakImgs(:,:,ix) =  imreconstruct(aptImg,a);
    % All we've done is eliminated the other peaks.
    assert(all(reshape(peakImgs(:,:,ix) <= a,[],1)));

for ix = 1 : N
    peakMask(:,:,ix) = ix*(max(peakImgs(:,:,setdiff(1:N, ix)),[],3) < a);

peakMaskResult = max(peakMask,[],3);

And the resulting segmentation: enter image description here

To Reiterate. The question is, "Is there a matlab implementation of the peak segmentation algorithm (as described above) that can be executed in O(1) of number of peaks?"


I'm not sure why you need to call imreconstruct more than once. The way you normally use morphological reconstruction is described here: You basically subtract some constant (the peak threshold) from your source image, then reconstruct it. What you get is an image that is identical to your original image except at the peaks:

enter image description here

So you simply take the difference and binarize, or compare the two images pixel-wise.

  • $\begingroup$ Consider my example. I find the three peaks in the center of each of my identified regions. I run imreconstruct on that image with those three peaks, the reconstructed image will be identical to the starting image. What you describe only works if I run imreconstruct on each peak individually. OR, on a carefully chosen subset of the peaks, so instead of N times, I could probably get by with log_2 (N) times. Definitely a major improvement, but not equal to 1 time. $\endgroup$
    – John
    Sep 5 '14 at 19:59
  • $\begingroup$ @John: That's what I said: You're using it wrong. You shouldn't run imreconstruct with the peaks at all. You run it with source image - threshold. Once. (You didn't post any real image data, so I can't show actual code or results.) $\endgroup$ Sep 5 '14 at 21:12
  • $\begingroup$ your approach presumes knowledge about the heights of the peaks which I don't have. I really need to segment out the entire peak, which I don't see how to do with imreconstruct. I've rewritten the question considerably, to hopefully clarify the issues. Thanks!! $\endgroup$
    – John
    Sep 7 '14 at 3:03
  • $\begingroup$ @John: My Matlab is little rusty, so I'm not sure I understand what you're doing 100%. (you didn't post a sample image I could try this on, either.) But it looks similar to watershed segmentation. (Watershed looks for basins, not peaks, so you'll have to invert your input image). $\endgroup$ Sep 8 '14 at 20:25
  • $\begingroup$ watershed!!!! That's what it's freaking called!! Excellent. Thanks. $\endgroup$
    – John
    Sep 9 '14 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.