# O(1) peak segmentation algorithm in matlab?

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.) The above image was rendered in matlab with the command

imagesc(peaks)


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)));
end

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

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


And the resulting segmentation: 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?"

## 1 Answer

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: So you simply take the difference and binarize, or compare the two images pixel-wise.

• 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. – John Sep 5 '14 at 19:59
• @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.) – Niki Estner Sep 5 '14 at 21:12
• 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!! – John Sep 7 '14 at 3:03
• @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). – Niki Estner Sep 8 '14 at 20:25
• watershed!!!! That's what it's freaking called!! Excellent. Thanks. – John Sep 9 '14 at 16:03