#14797 Add RigContourMapPeakFinder to find local maxima in contour maps

Union-find sweep from highest to lowest value computes the prominence of each peak. Peaks are filtered on prominence and edge contact, ranked by prominence or value, and thinned by minimum distance. Flat shoulders and plateaus are rejected, and prominence is measured against the minimum of each connected data region.

ApplicationLibCode/ReservoirDataModel/ContourMap/CMakeLists_files.cmake

@@ -4,6 +4,7 @@ set(SOURCE_GROUP_SOURCE_FILES
     ${CMAKE_CURRENT_LIST_DIR}/RigContourPolygonsTools.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RigContourMapCalculator.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RigContourMapProjection.cpp
+    ${CMAKE_CURRENT_LIST_DIR}/RigContourMapPeakFinder.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RigContourMapTopography.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RigContourMapTrianglesGenerator.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RigEclipseContourMapProjection.cpp

ApplicationLibCode/ReservoirDataModel/ContourMap/RigContourMapPeakFinder.cpp

@@ -0,0 +1,239 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2026- Equinor ASA
+//
+//  ResInsight is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
+//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.
+//
+//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#include "RigContourMapPeakFinder.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+namespace RigContourMapPeakFinder
+{
+namespace internal
+{
+    //--------------------------------------------------------------------------------------------------
+    /// Returns true if the cell at (i, j), or any of its 8 neighbors, is outside the grid or undefined.
+    //--------------------------------------------------------------------------------------------------
+    bool touchesEdgeOrUndefined( const std::vector<double>& z, int nx, int ny, int i, int j )
+    {
+        for ( int dj = -1; dj <= 1; ++dj )
+        {
+            for ( int di = -1; di <= 1; ++di )
+            {
+                const int ni = i + di;
+                const int nj = j + dj;
+                if ( ni < 0 || nj < 0 || ni >= nx || nj >= ny ) return true;
+                if ( !std::isfinite( z[ni + nj * nx] ) ) return true;
+            }
+        }
+        return false;
+    }
+
+    //--------------------------------------------------------------------------------------------------
+    /// Prominence per cell (NaN for non-peak cells, >= 0 for peaks). For each peak, also records whether
+    /// its region survived the sweep (i.e. it is the highest peak of a connected defined area) and the
+    /// number of cells in that surviving region.
+    ///
+    /// A zero prominence peak that was absorbed at a saddle is never a real peak: it is a flat shoulder
+    /// or plateau cell split off by the tie-break rule. A surviving zero prominence peak is a completely
+    /// flat connected area, which only counts as a peak if it is a single isolated cell.
+    //--------------------------------------------------------------------------------------------------
+    struct ProminenceResult
+    {
+        std::vector<double> prominence;
+        std::vector<char>   isSurvivingRoot;
+        std::vector<int>    regionSize;
+    };
+
+    ProminenceResult computeProminenceAndRegionSize( const std::vector<double>& z, int nx, int ny )
+    {
+        const int    M   = nx * ny;
+        const double nan = std::numeric_limits<double>::quiet_NaN();
+
+        ProminenceResult result;
+        result.prominence.assign( M, nan );
+        result.isSurvivingRoot.assign( M, 0 );
+        result.regionSize.assign( M, 0 );
+
+        std::vector<int> order;
+        order.reserve( M );
+
+        for ( int k = 0; k < M; ++k )
+        {
+            if ( std::isfinite( z[k] ) ) order.push_back( k );
+        }
+
+        // Descending by value, ties broken by index -> plateaus yield exactly one peak
+        std::sort( order.begin(), order.end(), [&]( int a, int b ) { return z[a] > z[b] || ( z[a] == z[b] && a < b ); } );
+
+        std::vector<int>    parent( M, -1 ); // -1 = not yet processed
+        std::vector<int>    peakOf( M, -1 ); // valid for roots: index of the region's highest cell
+        std::vector<int>    liveRegionSize( M, 0 ); // valid for roots: number of cells merged into the region so far
+        std::vector<double> regionMin( M, 0.0 ); // valid for roots: lowest value merged into the region so far
+
+        auto findRoot = [&]( int a )
+        {
+            while ( parent[a] != a )
+            {
+                parent[a] = parent[parent[a]]; // path halving
+                a         = parent[a];
+            }
+            return a;
+        };
+
+        auto higherPeak = [&]( int pa, int pb ) { return z[pa] > z[pb] || ( z[pa] == z[pb] && pa < pb ); };
+
+        std::vector<int> roots;
+        roots.reserve( 8 );
+
+        for ( int c : order )
+        {
+            parent[c]         = c;
+            peakOf[c]         = c;
+            liveRegionSize[c] = 1;
+            regionMin[c]      = z[c];
+
+            const int ci = c % nx;
+            const int cj = c / nx;
+
+            roots.clear();
+            for ( int dj = -1; dj <= 1; ++dj )
+            {
+                for ( int di = -1; di <= 1; ++di )
+                {
+                    if ( di == 0 && dj == 0 ) continue;
+                    const int ni = ci + di;
+                    const int nj = cj + dj;
+                    if ( ni < 0 || nj < 0 || ni >= nx || nj >= ny ) continue;
+
+                    const int n = ni + nj * nx;
+                    if ( parent[n] == -1 ) continue;
+
+                    const int r = findRoot( n );
+                    if ( std::find( roots.begin(), roots.end(), r ) == roots.end() ) roots.push_back( r );
+                }
+            }
+
+            if ( roots.empty() ) continue; // c is a new peak
+
+            // Dominant region = the one with the highest peak
+            int dominant = roots.front();
+            for ( int r : roots )
+            {
+                if ( higherPeak( peakOf[r], peakOf[dominant] ) ) dominant = r;
+            }
+
+            // c is the saddle for all other regions: they die here
+            for ( int r : roots )
+            {
+                if ( r == dominant ) continue;
+                result.prominence[peakOf[r]] = z[peakOf[r]] - z[c];
+                liveRegionSize[dominant] += liveRegionSize[r];
+                parent[r] = dominant;
+            }
+            liveRegionSize[dominant] += liveRegionSize[c];
+            regionMin[dominant] = std::min( regionMin[dominant], z[c] ); // cells are visited in descending order
+            parent[c]           = dominant;
+        }
+
+        // Surviving regions (one per connected defined area): prominence relative to the region's own minimum
+        for ( int k : order )
+        {
+            if ( parent[k] == k )
+            {
+                result.prominence[peakOf[k]]      = z[peakOf[k]] - regionMin[k];
+                result.isSurvivingRoot[peakOf[k]] = 1;
+                result.regionSize[peakOf[k]]      = liveRegionSize[k];
+            }
+        }
+
+        return result;
+    }
+} // namespace internal
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<double> computeProminence( const std::vector<double>& z, int nx, int ny )
+{
+    return internal::computeProminenceAndRegionSize( z, nx, ny ).prominence;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<Peak>
+    findPeaks( const std::vector<double>& z, int nx, int ny, double dx, double dy, const Settings& settings, double originX, double originY )
+{
+    const internal::ProminenceResult prominenceResult = internal::computeProminenceAndRegionSize( z, nx, ny );
+    const std::vector<double>&       prominence       = prominenceResult.prominence;
+    const std::vector<char>&         isSurvivingRoot  = prominenceResult.isSurvivingRoot;
+    const std::vector<int>&          regionSize       = prominenceResult.regionSize;
+
+    const int maxPeaks = std::max( 0, settings.maxPeaks );
+    if ( maxPeaks == 0 ) return {};
+
+    std::vector<Peak> candidates;
+    for ( int k = 0; k < nx * ny; ++k )
+    {
+        if ( std::isnan( prominence[k] ) || prominence[k] < settings.minProminence ) continue;
+
+        // Zero prominence: flat shoulders and plateaus are not peaks. Only a single isolated cell (no other
+        // data to compare against) is kept as a legitimate, if trivial, peak.
+        if ( prominence[k] <= 0.0 && !( isSurvivingRoot[k] && regionSize[k] == 1 ) ) continue;
+
+        const int i = k % nx;
+        const int j = k / nx;
+        if ( settings.excludeEdgePeaks && internal::touchesEdgeOrUndefined( z, nx, ny, i, j ) ) continue;
+
+        candidates.push_back( { i, j, originX + i * dx, originY + j * dy, z[k], prominence[k] } );
+    }
+
+    std::sort( candidates.begin(),
+               candidates.end(),
+               [&]( const Peak& a, const Peak& b )
+               {
+                   const double ka = settings.rankByProminence ? a.prominence : a.z;
+                   const double kb = settings.rankByProminence ? b.prominence : b.z;
+                   return ka > kb;
+               } );
+
+    // Greedy non-maximum suppression on distance
+    const double      minDist2 = settings.minDistance * settings.minDistance;
+    std::vector<Peak> result;
+    result.reserve( maxPeaks );
+
+    for ( const Peak& c : candidates )
+    {
+        if ( static_cast<int>( result.size() ) >= maxPeaks ) break;
+
+        const bool farEnough = std::none_of( result.begin(),
+                                             result.end(),
+                                             [&]( const Peak& t )
+                                             {
+                                                 const double ddx = c.x - t.x;
+                                                 const double ddy = c.y - t.y;
+                                                 return ddx * ddx + ddy * ddy < minDist2;
+                                             } );
+        if ( farEnough ) result.push_back( c );
+    }
+
+    return result;
+}
+
+} // namespace RigContourMapPeakFinder

ApplicationLibCode/ReservoirDataModel/ContourMap/RigContourMapPeakFinder.h

@@ -0,0 +1,65 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2026- Equinor ASA
+//
+//  ResInsight is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
+//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.
+//
+//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#pragma once
+
+#include <vector>
+
+//==================================================================================================
+/// Find the N most significant peaks in a gridded surface / contour map, with a minimum distance
+/// between accepted peaks.
+///
+/// Grid layout: row-major, index = i + j * nx. Undefined cells may be NaN or +/-inf.
+/// Distances are computed in world units via cell sizes dx, dy (rotation does not affect distance).
+//==================================================================================================
+namespace RigContourMapPeakFinder
+{
+struct Peak
+{
+    int    i          = 0;
+    int    j          = 0;
+    double x          = 0.0;
+    double y          = 0.0;
+    double z          = 0.0;
+    double prominence = 0.0;
+};
+
+struct Settings
+{
+    int    maxPeaks         = 10;
+    double minDistance      = 0.0; // world units
+    double minProminence    = 0.0; // z units, filters noise / shoulders
+    bool   rankByProminence = true; // false: rank by z (after prominence filter)
+    bool   excludeEdgePeaks = false; // drop peaks touching grid border or undefined cells
+};
+
+// Returns prominence per cell: NaN for non-peak cells, >= 0 for peaks.
+// Union-find sweep from highest to lowest value (8-connectivity). O(M log M).
+std::vector<double> computeProminence( const std::vector<double>& z, int nx, int ny );
+
+// Find the top N peaks, ranked and filtered according to settings.
+std::vector<Peak> findPeaks( const std::vector<double>& z,
+                             int                        nx,
+                             int                        ny,
+                             double                     dx,
+                             double                     dy,
+                             const Settings&            settings,
+                             double                     originX = 0.0,
+                             double                     originY = 0.0 );
+
+} // namespace RigContourMapPeakFinder

ApplicationLibCode/UnitTests/CMakeLists.txt

@@ -51,6 +51,7 @@ set(SOURCE_UNITTEST_FILES
     ${CMAKE_CURRENT_LIST_DIR}/WellPathAsciiFileReader-Test.cpp
     ${CMAKE_CURRENT_LIST_DIR}/opm-flowdiagnostics-Test.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RigTofAccumulatedPhaseFractionsCalculator-Test.cpp
+    ${CMAKE_CURRENT_LIST_DIR}/RigContourMapPeakFinder-Test.cpp
     ${CMAKE_CURRENT_LIST_DIR}/HDF5FileReader-Test.cpp
     ${CMAKE_CURRENT_LIST_DIR}/HDF5FileWriter-Test.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RigCellGeometryTools-Test.cpp

ApplicationLibCode/UnitTests/RigContourMapPeakFinder-Test.cpp

@@ -0,0 +1,323 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2026- Equinor ASA
+//
+//  ResInsight is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
+//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.
+//
+//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#include "gtest/gtest.h"
+
+#include "ContourMap/RigContourMapPeakFinder.h"
+
+#include <cmath>
+#include <limits>
+
+namespace
+{
+double nanValue()
+{
+    return std::numeric_limits<double>::quiet_NaN();
+}
+} // namespace
+
+//--------------------------------------------------------------------------------------------------
+/// Single, isolated peak in an otherwise flat grid.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, SinglePeak )
+{
+    const int nx = 5;
+    const int ny = 5;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[2 + 2 * nx] = 10.0; // center cell
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks = 5;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 1u, peaks.size() );
+    EXPECT_EQ( 2, peaks[0].i );
+    EXPECT_EQ( 2, peaks[0].j );
+    EXPECT_DOUBLE_EQ( 10.0, peaks[0].z );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Two separated peaks, both should be reported when far enough apart and maxPeaks allows it.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, TwoSeparatedPeaks )
+{
+    const int nx = 10;
+    const int ny = 1;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[1] = 5.0;
+    z[8] = 8.0;
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks    = 10;
+    settings.minDistance = 1.0;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 2u, peaks.size() );
+    // Ranked by prominence/z descending: highest peak first.
+    EXPECT_EQ( 8, peaks[0].i );
+    EXPECT_EQ( 1, peaks[1].i );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Two nearby peaks: with a large minimum distance only the more prominent one should survive.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, MinimumDistanceSuppressesNearbyPeak )
+{
+    const int nx = 10;
+    const int ny = 1;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[4] = 5.0;
+    z[5] = 8.0;
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks    = 10;
+    settings.minDistance = 5.0;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 1u, peaks.size() );
+    EXPECT_EQ( 5, peaks[0].i );
+    EXPECT_DOUBLE_EQ( 8.0, peaks[0].z );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// maxPeaks limits the number of returned peaks to the most significant ones.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, MaxPeaksLimitsResultCount )
+{
+    const int nx = 10;
+    const int ny = 1;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[1] = 3.0;
+    z[4] = 5.0;
+    z[7] = 8.0;
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks    = 2;
+    settings.minDistance = 0.0;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 2u, peaks.size() );
+    EXPECT_EQ( 7, peaks[0].i );
+    EXPECT_EQ( 4, peaks[1].i );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// minProminence filters out low-significance shoulders/noise.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, MinimumProminenceFiltersNoise )
+{
+    const int nx = 10;
+    const int ny = 1;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[4] = 0.5; // low prominence "noise"
+    z[8] = 8.0;
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks      = 10;
+    settings.minProminence = 1.0;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 1u, peaks.size() );
+    EXPECT_EQ( 8, peaks[0].i );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// excludeEdgePeaks drops peaks touching the grid border.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, ExcludeEdgePeaks )
+{
+    const int nx = 5;
+    const int ny = 5;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[0]          = 8.0; // corner cell, touches the edge
+    z[2 + 2 * nx] = 5.0; // interior cell
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks         = 10;
+    settings.excludeEdgePeaks = true;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 1u, peaks.size() );
+    EXPECT_EQ( 2, peaks[0].i );
+    EXPECT_EQ( 2, peaks[0].j );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Undefined (NaN) cells are ignored and must not be returned as peaks.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, UndefinedCellsAreIgnored )
+{
+    const int nx = 5;
+    const int ny = 5;
+
+    std::vector<double> z( nx * ny, nanValue() );
+    z[2 + 2 * nx] = 5.0;
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks = 10;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 1u, peaks.size() );
+    EXPECT_EQ( 2, peaks[0].i );
+    EXPECT_EQ( 2, peaks[0].j );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// World-space x/y coordinates should reflect origin and cell sizes.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, WorldCoordinatesUseOriginAndCellSize )
+{
+    const int nx = 5;
+    const int ny = 5;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[3 + 2 * nx] = 5.0;
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks = 10;
+
+    const double originX = 100.0;
+    const double originY = 200.0;
+    const double dx      = 10.0;
+    const double dy      = 20.0;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, dx, dy, settings, originX, originY );
+
+    ASSERT_EQ( 1u, peaks.size() );
+    EXPECT_DOUBLE_EQ( originX + 3 * dx, peaks[0].x );
+    EXPECT_DOUBLE_EQ( originY + 2 * dy, peaks[0].y );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// rankByProminence = false ranks candidates by raw z value instead of prominence. The data has a high
+/// peak with a high saddle next to the global maximum, and a lower but isolated (more prominent) peak,
+/// so the two ranking modes select different second peaks.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, RankByZValue )
+{
+    const int nx = 10;
+    const int ny = 1;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[1] = 9.0; // high z, low prominence (0.5) due to the high saddle at i=2
+    z[2] = 8.5;
+    z[3] = 10.0; // global maximum
+    z[7] = 5.0; // lower z, but isolated so higher prominence (5)
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks = 2;
+
+    settings.rankByProminence = false;
+    auto peaksByZ             = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 2u, peaksByZ.size() );
+    EXPECT_EQ( 3, peaksByZ[0].i );
+    EXPECT_EQ( 1, peaksByZ[1].i );
+    EXPECT_DOUBLE_EQ( 0.5, peaksByZ[1].prominence );
+
+    settings.rankByProminence = true;
+    auto peaksByProminence    = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 2u, peaksByProminence.size() );
+    EXPECT_EQ( 3, peaksByProminence[0].i );
+    EXPECT_EQ( 7, peaksByProminence[1].i );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Flat steps (shoulders) on a slope and flat background cells have zero prominence and must not be
+/// reported as peaks.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, FlatStepIsNotAPeak )
+{
+    const int nx = 7;
+    const int ny = 1;
+
+    std::vector<double> z = { 0.0, 0.0, 3.0, 3.0, 5.0, 0.0, 0.0 };
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks = 10;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    ASSERT_EQ( 1u, peaks.size() );
+    EXPECT_EQ( 4, peaks[0].i );
+    EXPECT_DOUBLE_EQ( 5.0, peaks[0].prominence );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Each separate data region measures prominence against its own minimum, not the global minimum.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, ProminenceIsPerDataRegion )
+{
+    const int nx = 5;
+    const int ny = 1;
+
+    std::vector<double> z = { 100.0, 101.0, nanValue(), 5.0, 6.0 };
+
+    auto prominence = RigContourMapPeakFinder::computeProminence( z, nx, ny );
+
+    EXPECT_DOUBLE_EQ( 1.0, prominence[1] );
+    EXPECT_DOUBLE_EQ( 1.0, prominence[4] );
+    EXPECT_TRUE( std::isnan( prominence[0] ) );
+    EXPECT_TRUE( std::isnan( prominence[3] ) );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// A negative maxPeaks must not crash and returns no peaks.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, NegativeMaxPeaksReturnsNoPeaks )
+{
+    const int nx = 5;
+    const int ny = 5;
+
+    std::vector<double> z( nx * ny, 0.0 );
+    z[2 + 2 * nx] = 10.0;
+
+    RigContourMapPeakFinder::Settings settings;
+    settings.maxPeaks = -1;
+
+    auto peaks = RigContourMapPeakFinder::findPeaks( z, nx, ny, 1.0, 1.0, settings );
+
+    EXPECT_TRUE( peaks.empty() );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// An empty grid should not crash and must return no peaks.
+//--------------------------------------------------------------------------------------------------
+TEST( RigContourMapPeakFinderTest, EmptyGridReturnsNoPeaks )
+{
+    std::vector<double> z;
+
+    RigContourMapPeakFinder::Settings settings;
+    auto                              peaks = RigContourMapPeakFinder::findPeaks( z, 0, 0, 1.0, 1.0, settings );
+
+    EXPECT_TRUE( peaks.empty() );
+}