2

I am after a MATLAB code which marks each point in a 2d array (image) if it is a local maximum over a 3x3 neighborhood.

The MATLAB code will be used to generate a c code using MATLAB Coder.

A vanilla code, which tries to avoid allocations is given by:

function [is_local_max] = isLocalMax(in_arr, local_buff, is_local_max)

% is_local_max initialized a false with same size as in_arr
% local_buff is a 3x3 array with same class as in_arr (single / double)

n_rows, n_cols = size(in_arr);

for j = 2:(n_cols - 1)
  for i = 2:(n_rows - 1)
    local_buff(:) = in_arr((i - 1):(i + 1), (j - 1):(j + 1));
    is_local_max(i, j) = in_arr(i, j) == max(local_buff, [], 'all');
  end
end

end

The generated code from the above function is indeed with no allocations. Yet it is extremely slow.

Is there a better MATLAB code which will generate better code?
Should I be aware of a specific configuration for such code?

Improve this question
6
  • You could create a 3D matrix by stacking translated versions of your array so that each slice stack(i,j,:) is your local_buff at (i,j). Then you can use max on the third dimension Commented May 9, 2024 at 10:44
  • There’s a good chance that max(in_arr((i - 1):(i + 1), (j - 1):(j + 1)), [], 'all') allows for better optimization. Otherwise try to explicitly compare each of the 8 neighbors to the central pixel. Also, instead of writing the result of the comparison to is_local_max directly, maybe it’s faster to only write when it’s true (since you said it’s zero-initialized). If all else fails, write the C code by hand… Commented May 9, 2024 at 14:09
  • @beaker That will allocate an output array (and possibly also an intermediate one?) Commented May 9, 2024 at 14:13
  • 1
    The problem with this algorithm is that it doesn’t only detect local maxima. If you have a region of equal values in the image it will be falsely marked as a local maximum. You need a much more complex algorithm to robustly recognize local maxima, unfortunately. But hopefully this is sufficient for your use case. Commented May 9, 2024 at 14:18
  • 2
    “Infinite shrink”? Do you mean ultimate erosion? Leaving only a single pixel per connected component? —— To distinguish local maxima from other plateaus you need to go around the boundary of the plateau and verify that all the neighbors not on the plateau are lower. It’s not trivial to implement, but not complicated either. I’ve got a fairly efficient implementation here: github.com/DIPlib/diplib/blob/master/src/morphology/maxima.cpp Commented May 9, 2024 at 16:52

1 Answer 1

Reset to default
3

Here is the vectorized version of computing the local maxima (based on your algorithm):

is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = ...
      in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 1:(n_cols - 2)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(2:(n_rows - 1), 1:(n_cols - 2)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 1:(n_cols - 2)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 3:(n_cols    )) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(2:(n_rows - 1), 3:(n_cols    )) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 3:(n_cols    ));

It may be made more cache friendly:

is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = ...
      in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 1:(n_cols - 2)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(2:(n_rows - 1), 1:(n_cols - 2)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 1:(n_cols - 2));
    
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 2:(n_cols - 1));
    
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 3:(n_cols    )) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(2:(n_rows - 1), 3:(n_cols    )) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 3:(n_cols    ));

You can also try this third version:

is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = ...
      in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 1:(n_cols - 2));
      
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(2:(n_rows - 1), 1:(n_cols - 2));
    
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 1:(n_cols - 2));
    
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 2:(n_cols - 1));
    
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 2:(n_cols - 1));
    
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 3:(n_cols    ));
    
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(2:(n_rows - 1), 3:(n_cols    ));
    
is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) = is_local_max(2:(n_rows - 1), 2:(n_cols - 1)) ...
    & in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(3:(n_rows    ), 3:(n_cols    ));

I cannot verify, but MATLAB should be able to generate code without memory allocation.

Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

You could also compute the max over the three sub-arrays in_arr(1:end-2,:), in_arr(2:end-1,:), in_arr(3:end,:), and in the result of that compute the max over the three sub-arrays along the other direction. That would save some work, indexing 6 sub-arrays instead of 9.
It seems a good idea. I would like someone compare their timing.
I think in_arr(2:(n_rows - 1), 2:(n_cols - 1)) >= in_arr(1:(n_rows - 2), 1:(n_cols - 2)) will generate allocations in matlab.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.