title: Find peak (regions) in 2D data
body: <p>I am looking to find <strong>peak regions</strong> in 2D data (if you will, grayscale images or 2D landscapes, created through a Hough transform). By <em>peak region</em> I mean a <strong>locally maximal peak</strong>, yet <strong><em>NOT a single point</em></strong> but a part of the surrounding <strong><em>contributing region</em></strong> that goes with it. I know, this is a vague definition, but maybe the word <em>mountain</em> or the images below will give you an intuition of what I mean. </p>
  
  <p>The peaks marked in red (1-4) are what I want, the ones in pink (5-6) examples for the "grey zone", where it would be okay if those smaller peaks are not found but also okay if they are.</p>
  
  <p><a href="https://i.stack.imgur.com/HdXT8.png" rel="noreferrer"><img src="https://i.stack.imgur.com/HdXT8.png" alt="an optimal result in 3D"></a></p>
  
  <p>Images contain between 1-20 peaked regions, different in height. The 2D data for above surf plot is shown below with a possible result (orange corresponds to Peak 1, green corresponds to Peak 2 a/b, ...). Single images for tests can be found in the description links:</p>
  
  <p><em>Image left</em>: <a href="https://i.stack.imgur.com/udIuy.png" rel="noreferrer">input image</a> - - - - <em>middle:</em> (okaish) <a href="https://i.stack.imgur.com/4GYrI.png" rel="noreferrer">result</a> - - - - <em>right</em>: result overlayed over image.</p>
  
  <p><a href="https://i.stack.imgur.com/vdOqO.png" rel="noreferrer"><img src="https://i.stack.imgur.com/vdOqO.png" alt="input, result and overlayed result"></a></p>
  
  <p>The result above was produced using simple thresholding (MATLAB code):</p>
  
  <pre><code>% thresh_scale = 15;                     % parameter: how many thresholding steps 
  % thresh_perc = 6;                       % parameter: threshold at which we clip
  thresh = multithresh(H,thresh_scale);    
  q_image = imquantize(H, thresh);         
  
  q_image(q_image &lt;= thresh_perc) = 0;     % regions under threshold are thrown away
  q_image(q_image &gt; thresh_perc) = 1;      % ... while all others are preserved
  q_image = imbinarize(q_image);           % binarize for further processing
  B = bwareaopen(q_image, nhood_minsize);  % Filter really small regions
  [L, L_num] = bwlabel(B); % &lt;- result     % Label connected components
  </code></pre>
  
  <p>Some values like these (15 and 6) often work fine if there are few similar peaks, but this isn't consistent if more peaks are present or they vary a lot. I mainly have two problems, that also can't be fixed by simply adjusting the parameters:</p>
  
  <ul>
  <li>higher peaks can mask lower (but clearly distinguishable) peaks. Since the threshold is relative to the highest peak, other peaks may fall below.</li>
  <li>in some cases the valley between two peaks is above the threshold, merging several peaks into one (as can be observed with Peak 2 a/b). </li>
  </ul>
  
  <p>I also don't want a huge region for a high peak, so the <em>peak region</em> should probably be defined as some percentage of the mountain. I figured instead of a global thresholding, I'd rather have a method that finds peak regions relative to their immediate environment. I've looked into mean-shift and MSER segmentation, but those seem to be suited for segmenting real images, not kind of synthetic data. </p>
  
  <p>Somehow I imagined filling a negative of the landscape with a certain amount of water would give me the regions I'm looking for: basins that fill and spread with how the surrounding regions are shaped. Like pouring water over below image and the resulting waterpools are the regions I'm looking for. </p>
  
  <p><a href="https://i.stack.imgur.com/sgPqv.png" rel="noreferrer"><img src="https://i.stack.imgur.com/sgPqv.png" alt="negative image (complement), ready to pour water into it"></a></p>
  
  <p>I thought that is what the floodfill or watershed algorithm do, but floodfill seems like something completely else and the watershed results are not at all what I had in mind, also when applying some preprocessing that I thought could help (clipped at 1/10):</p>
  
  <p><a href="https://i.stack.imgur.com/d99c0.png" rel="noreferrer"><img src="https://i.stack.imgur.com/d99c0.png" alt="watershed clipped with threshold 1/10"></a></p>
  
  <p>Or when using the same clipping threshold as with above example (clipped at 6/15):</p>
  
  <p><a href="https://i.stack.imgur.com/Qhfuv.png" rel="noreferrer"><img src="https://i.stack.imgur.com/Qhfuv.png" alt="watershed clipped with threshold 6/15"></a></p>
  
  <p>Produced with this code (MATLAB):</p>
  
  <pre><code>thresh = multithresh(H, 10);    % set to either 10 || 15 for the examples
  q_image = imquantize(H, thresh);
  mask = false(size(q_image));    % create clipping mask...
  mask(q_image &gt; 1) = true;       % ... to remove lowest 10% || lowest 6/15
                                  % show with: figure, imshow(mask);
  
  % OPTIONAL: Gaussian smoothing
  H = imgaussfilt(H, 2);  % apply before adding Inf values
  % OPTIONAL: H-minima transform
  H = imhmin(H, 10);      % parameter is threshold for suppressing shallow minima
  H = -H;                 % Complement the image
  H(~mask) = Inf;         % force "ground" pixels to Inf
  
  L = watershed(D);    
  L(~mask) = 0;                               % clip "ground" from result
  imshow(label2rgb(L,'lines',[.5 .5 .5]));    % show result
  </code></pre>
  
  <hr>
  
  <p><strong><em>My question now:</em></strong> <strong>Is there an algorithm that fills a landscape and gives me the resulting waterpools (for various amounts of water poured) to do what I've tried to achieve with above methods?</strong> Or any other suggestion is welcome. I'm implementing MATLAB (or if need be Python), but I can work with any code or pseude-code.</p>
  
  <p>To distinguish this from <a href="https://stackoverflow.com/questions/9587173/algorithm-to-find-peaks-in-2d-array">this question</a>, my maxima are not separated by zero-values. What I want is similar, yet none of the suggestions there are helpful (hill-climbing/simulated annealing will give you only one point...).</p>
  
  <p><a href="https://stackoverflow.com/questions/3684484/peak-detection-in-a-2d-array">This question</a> is also interesting, but it solves the problem with constraints (assume exactly 5 peaks of a certain size) that make the suggested approaches not useful for my case.</p>