This article describes a practical approach to detect local maxima (and by extension, minima) in multidimensional arays using MATLAB's Image Processing Toolbox. The method leverages morphological dilation to identify peak locations while allowing control over minimum separation between detected extrema and handling of plateaus.

Core Concept

The algorithm uses imdilate() with a structuring element sized according to the desired minimum distance between peaks. By comparing the original data with its dilated version, points where the values are equal indicate local maxima within the neighborhood defined by the structuring element.

Function Parameters

• x: Input data array (any dimension)
• minDist: Minimum distance between peaks per dimension (scalar or vector)
• excludeEqualPoints: Boolean flag to handle flat regions (plateaus)

Handling Plateaus

When excludeEqualPoints is enabled, the function adds minimal random noise—smaller than the smallest non-zero difference in the data—to break ties in flat regions. This ensures only one representative point is selected from a plateau.

Usage Examples

1D Signal Example

a = cumsum(randn(1000,1));
peaks = localMaximum(a, [100]);
figure; plot(a); hold on; plot(peaks, a(peaks), 'ro');

2D Surface Example

[x, y] = meshgrid(-6:0.1:6, -6:0.1:6);
data = sinc(x) .* sinc(y);

maxIndices = localMaximum(data, [20 20]);
minIndices = localMaximum(-data, [20 20]);

figure;
mesh(x, y, data); hold on;
plot3(x(maxIndices), y(maxIndices), data(maxIndices), 'k*', 'MarkerSize', 10, 'LineWidth', 1.5);
plot3(x(minIndices), y(minIndices), data(minIndices), 'g*', 'MarkerSize', 10, 'LineWidth', 1.5);
legend('sinc(x)sinc(y)', 'Local Maxima', 'Local Minima', 'Location', 'Best');

Implementation Details

The function automatically adjusts minDist to match the input dimensions and enforces valid bounds. It returns linear indices by default but can output subscript coordinates when multiple outputs are requested.

function varargout = localMaximum(x, minDist, excludeEqualPoints)
    if nargin < 3
        excludeEqualPoints = false;
        if nargin < 2
            minDist = size(x) / 10;
        end
    end

    if isempty(minDist)
        minDist = size(x) / 10;
    end

    dimX = length(size(x));
    if length(minDist) ~= dimX
        minDist = repmat(minDist(1), [1, dimX]);
    end

    % Ensure valid structuring element size
    minDist = ceil(minDist);
    minDist = max([minDist(:)'; ones(1, length(minDist))]);
    minDist = min([minDist; size(x)]);

    % Handle plateaus by adding sub-resolution noise
    if excludeEqualPoints
        sortedVals = sort(x(:));
        diffs = diff(sortedVals);
        nonzeroDiffs = diffs(diffs ~= 0);
        if ~isempty(nonzeroDiffs)
            minStep = min(nonzeroDiffs);
            x = x + rand(size(x)) * (minStep / 10);
        end
    end

    % Create structuring element and perform dilation
    se = ones(minDist);
    dilated = imdilate(x, se);
    
    % Find local maxima
    peakIndices = find(x == dilated);

    % Return appropriate index format
    if nargout > 1
        [varargout{1:nargout}] = ind2sub(size(x), peakIndices);
    else
        varargout{1} = peakIndices;
    end
end

Notes

• Requires MATLAB's Image Processing Toolbox
• Pre-filtering data with a low-pass filter is recommended to reduce noise-induced false peaks
• To find local minima, apply the function to -x