diff --git a/autotest/cpp/CMakeLists.txt b/autotest/cpp/CMakeLists.txt
index 493f1d59f9bd..d8420177a9ff 100644
--- a/autotest/cpp/CMakeLists.txt
+++ b/autotest/cpp/CMakeLists.txt
@@ -77,6 +77,7 @@ add_executable(
   test_gdal_aaigrid.cpp
   test_gdal_dted.cpp
   test_gdal_gtiff.cpp
+  test_gdal_minmax_element.cpp
   test_gdal_pixelfn.cpp
   test_gdal_typetraits.cpp
   test_ogr.cpp
diff --git a/autotest/cpp/test_gdal_minmax_element.cpp b/autotest/cpp/test_gdal_minmax_element.cpp
new file mode 100644
index 000000000000..ec92329797d2
--- /dev/null
+++ b/autotest/cpp/test_gdal_minmax_element.cpp
@@ -0,0 +1,614 @@
+///////////////////////////////////////////////////////////////////////////////
+//
+// Project:  C++ Test Suite for GDAL/OGR
+// Purpose:  Test gdal_minmax_element.hpp
+// Author:   Even Rouault <even.rouault at spatialys.com>
+//
+///////////////////////////////////////////////////////////////////////////////
+// Copyright (c) 2023, Even Rouault <even.rouault at spatialys.com>
+/*
+ * SPDX-License-Identifier: MIT
+ ****************************************************************************/
+
+#include "gdal_unit_test.h"
+
+#include "gdal_minmax_element.hpp"
+
+#include "gtest_include.h"
+
+#include <limits>
+
+namespace
+{
+
+struct test_gdal_minmax_element : public ::testing::Test
+{
+};
+
+TEST_F(test_gdal_minmax_element, uint8)
+{
+    using T = uint8_t;
+    constexpr GDALDataType eDT = GDT_Byte;
+    T min_v = 1;
+    T max_v = 3;
+    {
+        T nodata = 0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T nodata = 0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v{static_cast<T>((min_v + max_v) / 2), max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v(129, static_cast<T>((min_v + max_v) / 2));
+        v[5] = min_v;
+        v[31] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, int8)
+{
+    using T = int8_t;
+    T min_v = -1;
+    T max_v = 3;
+    constexpr GDALDataType eDT = GDT_Int8;
+    {
+        T nodata = 0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T nodata = 0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v{static_cast<T>((min_v + max_v) / 2), max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v(129, static_cast<T>((min_v + max_v) / 2));
+        v[5] = min_v;
+        v[31] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, uint16)
+{
+    using T = uint16_t;
+    constexpr GDALDataType eDT = GDT_UInt16;
+    T min_v = 1000;
+    T max_v = 2000;
+    {
+        T nodata = 0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T nodata = 0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v{static_cast<T>((min_v + max_v) / 2), max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v(129, static_cast<T>((min_v + max_v) / 2));
+        v[5] = min_v;
+        v[31] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, int16)
+{
+    using T = int16_t;
+    constexpr GDALDataType eDT = GDT_Int16;
+    T min_v = -1000;
+    T max_v = 2000;
+    {
+        T nodata = 0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T nodata = 0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v{static_cast<T>((min_v + max_v) / 2), max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v(129, static_cast<T>((min_v + max_v) / 2));
+        v[5] = min_v;
+        v[31] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, uint32)
+{
+    using T = uint32_t;
+    constexpr GDALDataType eDT = GDT_UInt32;
+    T min_v = 10000000;
+    T max_v = 20000000;
+    {
+        T nodata = 0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T nodata = 0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v{static_cast<T>((min_v + max_v) / 2), max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v(129, static_cast<T>((min_v + max_v) / 2));
+        v[5] = min_v;
+        v[31] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, int32)
+{
+    using T = int32_t;
+    constexpr GDALDataType eDT = GDT_Int32;
+    T min_v = -10000000;
+    T max_v = 20000000;
+    {
+        T nodata = 0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T nodata = 0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v{static_cast<T>((min_v + max_v) / 2), max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v(129, static_cast<T>((min_v + max_v) / 2));
+        v[5] = min_v;
+        v[31] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, uint64)
+{
+    using T = uint64_t;
+    constexpr GDALDataType eDT = GDT_UInt64;
+    T min_v = 100000000000000;
+    T max_v = 200000000000000;
+    {
+        T nodata = 0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T nodata = 0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v{static_cast<T>((min_v + max_v) / 2), max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v(129, static_cast<T>((min_v + max_v) / 2));
+        v[5] = min_v;
+        v[31] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, int64)
+{
+    using T = int64_t;
+    constexpr GDALDataType eDT = GDT_Int64;
+    T min_v = -100000000000000;
+    T max_v = 200000000000000;
+    {
+        T nodata = 0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T nodata = 0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v{static_cast<T>((min_v + max_v) / 2), max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        std::vector<T> v(129, static_cast<T>((min_v + max_v) / 2));
+        v[5] = min_v;
+        v[31] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, float32)
+{
+    using T = float;
+    constexpr GDALDataType eDT = GDT_Float32;
+    {
+        T min_v = 1.0f;
+        T max_v = 1.5f;
+        T nodata = 2.0f;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T min_v = 1.0f;
+        T max_v = 1.5f;
+        T nodata = 2.0f;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0f;
+        T max_v = 1.5f;
+        T nodata = 2.0f;
+        std::vector<T> v{std::numeric_limits<T>::quiet_NaN(),
+                         std::numeric_limits<T>::quiet_NaN(), nodata, max_v,
+                         min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0f;
+        T max_v = 1.5f;
+        std::vector<T> v{std::numeric_limits<T>::quiet_NaN(),
+                         std::numeric_limits<T>::quiet_NaN(),
+                         max_v,
+                         std::numeric_limits<T>::quiet_NaN(),
+                         min_v,
+                         std::numeric_limits<T>::quiet_NaN()};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0f;
+        T max_v = 1.5f;
+        std::vector<T> v{max_v, std::numeric_limits<T>::quiet_NaN(), min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0f;
+        T max_v = 1.5f;
+        std::vector<T> v(33, std::numeric_limits<T>::quiet_NaN());
+        v[5] = min_v;
+        v[15] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0f;
+        T max_v = 1.5f;
+        std::vector<T> v(33, 1.2f);
+        v[5] = min_v;
+        v[15] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0f;
+        T max_v = 1.5f;
+        std::vector<T> v(33, std::numeric_limits<T>::quiet_NaN());
+        v[v.size() - 2] = min_v;
+        v.back() = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+TEST_F(test_gdal_minmax_element, float64)
+{
+    using T = double;
+    constexpr GDALDataType eDT = GDT_Float64;
+    {
+        T min_v = 1.0;
+        T max_v = 1.5;
+        T nodata = 2.0;
+        std::vector<T> v{max_v, nodata, min_v};
+        {
+            auto idx_min =
+                gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            auto idx_max =
+                gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+        {
+            auto [idx_min, idx_max] =
+                gdal::minmax_element(v.data(), v.size(), eDT, true, nodata);
+            EXPECT_EQ(v[idx_min], min_v);
+            EXPECT_EQ(v[idx_max], max_v);
+        }
+    }
+    {
+        T min_v = 1.0;
+        T max_v = 1.5;
+        T nodata = 2.0;
+        std::vector<T> v{nodata, max_v, min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0;
+        T max_v = 1.5;
+        T nodata = 2.0;
+        std::vector<T> v{std::numeric_limits<T>::quiet_NaN(),
+                         std::numeric_limits<T>::quiet_NaN(), nodata, max_v,
+                         min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, true, nodata);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0;
+        T max_v = 1.5;
+        std::vector<T> v{max_v, std::numeric_limits<T>::quiet_NaN(), min_v};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0;
+        T max_v = 1.5;
+        std::vector<T> v{std::numeric_limits<T>::quiet_NaN(),
+                         std::numeric_limits<T>::quiet_NaN(),
+                         max_v,
+                         std::numeric_limits<T>::quiet_NaN(),
+                         min_v,
+                         std::numeric_limits<T>::quiet_NaN()};
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+    {
+        T min_v = 1.0;
+        T max_v = 1.5;
+        std::vector<T> v(33, std::numeric_limits<T>::quiet_NaN());
+        v[5] = min_v;
+        v[15] = max_v;
+        auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_min], min_v);
+        auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0);
+        EXPECT_EQ(v[idx_max], max_v);
+    }
+}
+
+}  // namespace
diff --git a/gcore/CMakeLists.txt b/gcore/CMakeLists.txt
index 57a114bc3205..08fe7e646b7b 100644
--- a/gcore/CMakeLists.txt
+++ b/gcore/CMakeLists.txt
@@ -209,6 +209,7 @@ target_public_header(
   gdalsubdatasetinfo.h
   gdal_typetraits.h
   gdal_adbc.h
+  gdal_minmax_element.hpp
 )
 
 set(GDAL_DATA_FILES
diff --git a/gcore/gdal_minmax_element.hpp b/gcore/gdal_minmax_element.hpp
new file mode 100644
index 000000000000..490f809185dc
--- /dev/null
+++ b/gcore/gdal_minmax_element.hpp
@@ -0,0 +1,981 @@
+/******************************************************************************
+ * Project:  GDAL Core
+ * Purpose:  Utility functions to find minimum and maximum values in a buffer
+ * Author:   Even Rouault, <even dot rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2024, Even Rouault <even dot rouault at spatialys.com>
+ *
+ * SPDX-License-Identifier: MIT
+ ****************************************************************************/
+
+#ifndef GDAL_MINMAX_ELEMENT_INCLUDED
+#define GDAL_MINMAX_ELEMENT_INCLUDED
+
+// NOTE: This header requires C++17
+
+// This file may be vendored by other applications than GDAL
+// WARNING: if modifying this file, please also update the upstream GDAL version
+// at https://github.com/OSGeo/gdal/blob/master/gcore/gdal_minmax_element.hpp
+
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <utility>
+
+#include "gdal.h"
+
+#ifdef GDAL_COMPILATION
+#define GDAL_MINMAXELT_NS gdal
+#elif !defined(GDAL_MINMAXELT_NS)
+#error "Please define the GDAL_MINMAXELT_NS macro to define the namespace"
+#endif
+
+#if defined(__x86_64) || defined(_M_X64)
+// SSE2 header
+#include <emmintrin.h>
+#endif
+
+#include "gdal_priv_templates.hpp"
+#if GDAL_VERSION < GDAL_COMPUTE_VERSION(3, 10, 0)
+// For vendoring in other applications
+namespace GDAL_MINMAXELT_NS
+{
+template <class T> inline bool GDALIsValueExactAs(double dfValue)
+{
+    return GDALIsValueInRange<T>(dfValue) &&
+           static_cast<double>(static_cast<T>(dfValue)) == dfValue;
+}
+
+template <> inline bool GDALIsValueExactAs<float>(double dfValue)
+{
+    return std::isnan(dfValue) ||
+           (GDALIsValueInRange<float>(dfValue) &&
+            static_cast<double>(static_cast<float>(dfValue)) == dfValue);
+}
+
+template <> inline bool GDALIsValueExactAs<double>(double)
+{
+    return true;
+}
+}  // namespace GDAL_MINMAXELT_NS
+#endif
+
+namespace GDAL_MINMAXELT_NS
+{
+namespace detail
+{
+/************************************************************************/
+/*                         extremum_element()                           */
+/************************************************************************/
+
+template <class T, bool IS_MAX>
+size_t extremum_element(const T *v, size_t size, T noDataValue)
+{
+    static_assert(
+        !(std::is_same<T, float>::value || std::is_same<T, double>::value));
+    if (size == 0)
+        return 0;
+    size_t idx_of_extremum = 0;
+    T extremum = v[0];
+    bool extremum_is_nodata = extremum == noDataValue;
+    size_t i = 1;
+    for (; i < size; ++i)
+    {
+        if (v[i] != noDataValue &&
+            (((IS_MAX && v[i] > extremum) || (!IS_MAX && v[i] < extremum)) ||
+             extremum_is_nodata))
+        {
+            extremum = v[i];
+            idx_of_extremum = i;
+            extremum_is_nodata = false;
+        }
+    }
+    return idx_of_extremum;
+}
+
+/************************************************************************/
+/*                         extremum_element()                           */
+/************************************************************************/
+
+template <class T, bool IS_MAX> size_t extremum_element(const T *v, size_t size)
+{
+    static_assert(
+        !(std::is_same<T, float>::value || std::is_same<T, double>::value));
+    if (size == 0)
+        return 0;
+    size_t idx_of_extremum = 0;
+    T extremum = v[0];
+    size_t i = 1;
+    for (; i < size; ++i)
+    {
+        if ((IS_MAX && v[i] > extremum) || (!IS_MAX && v[i] < extremum))
+        {
+            extremum = v[i];
+            idx_of_extremum = i;
+        }
+    }
+    return idx_of_extremum;
+}
+
+#if defined(__x86_64) || defined(_M_X64)
+
+/************************************************************************/
+/*                    extremum_element_with_nan()                       */
+/************************************************************************/
+
+static inline int8_t Shift8(uint8_t x)
+{
+    return static_cast<int8_t>(x + std::numeric_limits<int8_t>::min());
+}
+
+static inline int16_t Shift16(uint16_t x)
+{
+    return static_cast<int16_t>(x + std::numeric_limits<int16_t>::min());
+}
+
+CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW
+static inline int32_t Shift32(uint32_t x)
+{
+    x += static_cast<uint32_t>(std::numeric_limits<int32_t>::min());
+    int32_t ret;
+    memcpy(&ret, &x, sizeof(x));
+    return ret;
+}
+
+// Using SSE2
+template <class T, bool IS_MAX>
+inline size_t extremum_element_with_nan(const T *v, size_t size)
+{
+    static_assert(
+        std::is_same<T, uint8_t>::value || std::is_same<T, int8_t>::value ||
+        std::is_same<T, uint16_t>::value || std::is_same<T, int16_t>::value ||
+        std::is_same<T, uint32_t>::value || std::is_same<T, int32_t>::value ||
+        std::is_same<T, float>::value || std::is_same<T, double>::value);
+    if (size == 0)
+        return 0;
+    size_t idx_of_extremum = 0;
+    T extremum = v[0];
+    [[maybe_unused]] bool extremum_is_nan = std::isnan(extremum);
+    size_t i = 1;
+    constexpr auto set1 = [](T x)
+    {
+        if constexpr (std::is_same<T, uint8_t>::value)
+            return _mm_set1_epi8(Shift8(x));
+        else if constexpr (std::is_same<T, int8_t>::value)
+            return _mm_set1_epi8(x);
+        else if constexpr (std::is_same<T, uint16_t>::value)
+            return _mm_set1_epi16(Shift16(x));
+        else if constexpr (std::is_same<T, int16_t>::value)
+            return _mm_set1_epi16(x);
+        else if constexpr (std::is_same<T, uint32_t>::value)
+            return _mm_set1_epi32(Shift32(x));
+        else if constexpr (std::is_same<T, int32_t>::value)
+            return _mm_set1_epi32(x);
+        else if constexpr (std::is_same<T, float>::value)
+            return _mm_set1_ps(x);
+        else
+            return _mm_set1_pd(x);
+    };
+
+    constexpr size_t VALS_PER_REG = sizeof(set1(extremum)) / sizeof(extremum);
+    constexpr int LOOP_UNROLLING = 4;
+    // If changing the value, then we need to adjust the number of sse_valX
+    // loading in the loop.
+    static_assert(LOOP_UNROLLING == 4);
+    constexpr size_t VALS_PER_ITER = VALS_PER_REG * LOOP_UNROLLING;
+
+    for (; i < VALS_PER_ITER && i < size; ++i)
+    {
+        if ((IS_MAX && v[i] > extremum) || (!IS_MAX && v[i] < extremum))
+        {
+            extremum = v[i];
+            idx_of_extremum = i;
+            extremum_is_nan = false;
+        }
+        else if constexpr (std::is_same<T, float>::value ||
+                           std::is_same<T, double>::value)
+        {
+            if (extremum_is_nan && !std::isnan(v[i]))
+            {
+                extremum = v[i];
+                idx_of_extremum = i;
+                extremum_is_nan = false;
+            }
+        }
+    }
+
+    auto sse_extremum = set1(extremum);
+
+    const auto loadv = [](const T *x)
+    {
+        if constexpr (std::is_same<T, float>::value)
+            return _mm_loadu_ps(x);
+        else if constexpr (std::is_same<T, double>::value)
+            return _mm_loadu_pd(x);
+        else
+            return _mm_loadu_si128(reinterpret_cast<const __m128i *>(x));
+    };
+
+    const auto comp = [](const auto x, const auto y)
+    {
+        if constexpr (IS_MAX)
+        {
+            if constexpr (std::is_same<T, uint8_t>::value)
+                return _mm_cmpgt_epi8(
+                    _mm_add_epi8(
+                        x, _mm_set1_epi8(std::numeric_limits<int8_t>::min())),
+                    y);
+            else if constexpr (std::is_same<T, int8_t>::value)
+                return _mm_cmpgt_epi8(x, y);
+            else if constexpr (std::is_same<T, uint16_t>::value)
+                return _mm_cmpgt_epi16(
+                    _mm_add_epi16(
+                        x, _mm_set1_epi16(std::numeric_limits<int16_t>::min())),
+                    y);
+            else if constexpr (std::is_same<T, int16_t>::value)
+                return _mm_cmpgt_epi16(x, y);
+            else if constexpr (std::is_same<T, uint32_t>::value)
+                return _mm_cmpgt_epi32(
+                    _mm_add_epi32(
+                        x, _mm_set1_epi32(std::numeric_limits<int32_t>::min())),
+                    y);
+            else if constexpr (std::is_same<T, int32_t>::value)
+                return _mm_cmpgt_epi32(x, y);
+            // We could use _mm_cmpgt_pX() if there was no NaN values
+            else if constexpr (std::is_same<T, float>::value)
+                return _mm_castps_si128(_mm_cmpnle_ps(x, y));
+            else
+                return _mm_castpd_si128(_mm_cmpnle_pd(x, y));
+        }
+        else
+        {
+            if constexpr (std::is_same<T, uint8_t>::value)
+                return _mm_cmplt_epi8(
+                    _mm_add_epi8(
+                        x, _mm_set1_epi8(std::numeric_limits<int8_t>::min())),
+                    y);
+            else if constexpr (std::is_same<T, int8_t>::value)
+                return _mm_cmplt_epi8(x, y);
+            else if constexpr (std::is_same<T, uint16_t>::value)
+                return _mm_cmplt_epi16(
+                    _mm_add_epi16(
+                        x, _mm_set1_epi16(std::numeric_limits<int16_t>::min())),
+                    y);
+            else if constexpr (std::is_same<T, int16_t>::value)
+                return _mm_cmplt_epi16(x, y);
+            else if constexpr (std::is_same<T, uint32_t>::value)
+                return _mm_cmplt_epi32(
+                    _mm_add_epi32(
+                        x, _mm_set1_epi32(std::numeric_limits<int32_t>::min())),
+                    y);
+            else if constexpr (std::is_same<T, int32_t>::value)
+                return _mm_cmplt_epi32(x, y);
+            // We could use _mm_cmplt_pX() if there was no NaN values
+            else if constexpr (std::is_same<T, float>::value)
+                return _mm_castps_si128(_mm_cmpnge_ps(x, y));
+            else
+                return _mm_castpd_si128(_mm_cmpnge_pd(x, y));
+        }
+    };
+
+    [[maybe_unused]] size_t hits = 0;
+    const auto sse_iter_count = (size / VALS_PER_ITER) * VALS_PER_ITER;
+    for (; i < sse_iter_count; i += VALS_PER_ITER)
+    {
+        // A bit of loop unrolling to save 3/4 of slow movemask operations.
+        const auto sse_val0 = loadv(v + i + 0 * VALS_PER_REG);
+        const auto sse_val1 = loadv(v + i + 1 * VALS_PER_REG);
+        const auto sse_val2 = loadv(v + i + 2 * VALS_PER_REG);
+        const auto sse_val3 = loadv(v + i + 3 * VALS_PER_REG);
+        if (_mm_movemask_epi8(
+                _mm_or_si128(_mm_or_si128(comp(sse_val0, sse_extremum),
+                                          comp(sse_val1, sse_extremum)),
+                             _mm_or_si128(comp(sse_val2, sse_extremum),
+                                          comp(sse_val3, sse_extremum)))) != 0)
+        {
+            if constexpr (!std::is_same<T, uint8_t>::value)
+            {
+                if (++hits == size / 16)
+                {
+                    // If we have an almost sorted array, then using this code path
+                    // will hurt performance. Arbitrary give up if we get here
+                    // more than 1. / 16 of the size of the array.
+                    // fprintf(stderr, "going to non-vector path\n");
+                    break;
+                }
+            }
+            for (size_t j = 0; j < VALS_PER_ITER; j++)
+            {
+                if ((IS_MAX && v[i + j] > extremum) ||
+                    (!IS_MAX && v[i + j] < extremum))
+                {
+                    extremum = v[i + j];
+                    idx_of_extremum = i + j;
+                    extremum_is_nan = false;
+                }
+                else if constexpr (std::is_same<T, float>::value ||
+                                   std::is_same<T, double>::value)
+                {
+                    if (extremum_is_nan && !std::isnan(v[i + j]))
+                    {
+                        extremum = v[i + j];
+                        idx_of_extremum = i + j;
+                        extremum_is_nan = false;
+                    }
+                }
+            }
+            sse_extremum = set1(extremum);
+        }
+    }
+    for (; i < size; ++i)
+    {
+        if ((IS_MAX && v[i] > extremum) || (!IS_MAX && v[i] < extremum))
+        {
+            extremum = v[i];
+            idx_of_extremum = i;
+            extremum_is_nan = false;
+        }
+        else if constexpr (std::is_same<T, float>::value ||
+                           std::is_same<T, double>::value)
+        {
+            if (extremum_is_nan && !std::isnan(v[i]))
+            {
+                extremum = v[i];
+                idx_of_extremum = i;
+                extremum_is_nan = false;
+            }
+        }
+    }
+    return idx_of_extremum;
+}
+
+#else
+
+/************************************************************************/
+/*                    extremum_element_with_nan()                       */
+/************************************************************************/
+
+template <class T, bool IS_MAX>
+inline size_t extremum_element_with_nan(const T *v, size_t size)
+{
+    if (size == 0)
+        return 0;
+    size_t idx_of_extremum = 0;
+    auto extremum = v[0];
+    bool extremum_is_nan = std::isnan(extremum);
+    size_t i = 1;
+    for (; i < size; ++i)
+    {
+        if ((IS_MAX && v[i] > extremum) || (!IS_MAX && v[i] < extremum) ||
+            (extremum_is_nan && !std::isnan(v[i])))
+        {
+            extremum = v[i];
+            idx_of_extremum = i;
+            extremum_is_nan = false;
+        }
+    }
+    return idx_of_extremum;
+}
+#endif
+
+/************************************************************************/
+/*                         extremum_element()                           */
+/************************************************************************/
+
+#if defined(__x86_64) || defined(_M_X64)
+
+template <>
+size_t extremum_element<uint8_t, true>(const uint8_t *v, size_t size)
+{
+    return extremum_element_with_nan<uint8_t, true>(v, size);
+}
+
+template <>
+size_t extremum_element<uint8_t, false>(const uint8_t *v, size_t size)
+{
+    return extremum_element_with_nan<uint8_t, false>(v, size);
+}
+
+template <> size_t extremum_element<int8_t, true>(const int8_t *v, size_t size)
+{
+    return extremum_element_with_nan<int8_t, true>(v, size);
+}
+
+template <> size_t extremum_element<int8_t, false>(const int8_t *v, size_t size)
+{
+    return extremum_element_with_nan<int8_t, false>(v, size);
+}
+
+template <>
+size_t extremum_element<uint16_t, true>(const uint16_t *v, size_t size)
+{
+    return extremum_element_with_nan<uint16_t, true>(v, size);
+}
+
+template <>
+size_t extremum_element<uint16_t, false>(const uint16_t *v, size_t size)
+{
+    return extremum_element_with_nan<uint16_t, false>(v, size);
+}
+
+template <>
+size_t extremum_element<int16_t, true>(const int16_t *v, size_t size)
+{
+    return extremum_element_with_nan<int16_t, true>(v, size);
+}
+
+template <>
+size_t extremum_element<int16_t, false>(const int16_t *v, size_t size)
+{
+    return extremum_element_with_nan<int16_t, false>(v, size);
+}
+
+template <>
+size_t extremum_element<uint32_t, true>(const uint32_t *v, size_t size)
+{
+    return extremum_element_with_nan<uint32_t, true>(v, size);
+}
+
+template <>
+size_t extremum_element<uint32_t, false>(const uint32_t *v, size_t size)
+{
+    return extremum_element_with_nan<uint32_t, false>(v, size);
+}
+
+template <>
+size_t extremum_element<int32_t, true>(const int32_t *v, size_t size)
+{
+    return extremum_element_with_nan<int32_t, true>(v, size);
+}
+
+template <>
+size_t extremum_element<int32_t, false>(const int32_t *v, size_t size)
+{
+    return extremum_element_with_nan<int32_t, false>(v, size);
+}
+
+#endif
+
+template <> size_t extremum_element<float, true>(const float *v, size_t size)
+{
+    return extremum_element_with_nan<float, true>(v, size);
+}
+
+template <> size_t extremum_element<double, true>(const double *v, size_t size)
+{
+    return extremum_element_with_nan<double, true>(v, size);
+}
+
+template <> size_t extremum_element<float, false>(const float *v, size_t size)
+{
+    return extremum_element_with_nan<float, false>(v, size);
+}
+
+template <> size_t extremum_element<double, false>(const double *v, size_t size)
+{
+    return extremum_element_with_nan<double, false>(v, size);
+}
+
+/************************************************************************/
+/*                       extremum_element_with_nan()                    */
+/************************************************************************/
+
+template <class T, bool IS_MAX>
+inline size_t extremum_element_with_nan(const T *v, size_t size, T noDataValue)
+{
+    if (std::isnan(noDataValue))
+        return extremum_element_with_nan<T, IS_MAX>(v, size);
+    if (size == 0)
+        return 0;
+    size_t idx_of_extremum = 0;
+    auto extremum = v[0];
+    bool extremum_is_nan_or_nodata =
+        std::isnan(extremum) || (extremum == noDataValue);
+    size_t i = 1;
+    for (; i < size; ++i)
+    {
+        if (v[i] != noDataValue &&
+            ((IS_MAX && v[i] > extremum) || (!IS_MAX && v[i] < extremum) ||
+             (extremum_is_nan_or_nodata && !std::isnan(v[i]))))
+        {
+            extremum = v[i];
+            idx_of_extremum = i;
+            extremum_is_nan_or_nodata = false;
+        }
+    }
+    return idx_of_extremum;
+}
+
+/************************************************************************/
+/*                            extremum_element()                        */
+/************************************************************************/
+
+template <>
+size_t extremum_element<float, true>(const float *v, size_t size,
+                                     float noDataValue)
+{
+    return extremum_element_with_nan<float, true>(v, size, noDataValue);
+}
+
+template <>
+size_t extremum_element<double, true>(const double *v, size_t size,
+                                      double noDataValue)
+{
+    return extremum_element_with_nan<double, true>(v, size, noDataValue);
+}
+
+template <>
+size_t extremum_element<float, false>(const float *v, size_t size,
+                                      float noDataValue)
+{
+    return extremum_element_with_nan<float, false>(v, size, noDataValue);
+}
+
+template <>
+size_t extremum_element<double, false>(const double *v, size_t size,
+                                       double noDataValue)
+{
+    return extremum_element_with_nan<double, false>(v, size, noDataValue);
+}
+
+template <class T, bool IS_MAX>
+inline size_t extremum_element(const T *buffer, size_t size, bool bHasNoData,
+                               T noDataValue)
+{
+    if (bHasNoData)
+        return extremum_element<T, IS_MAX>(buffer, size, noDataValue);
+    else
+        return extremum_element<T, IS_MAX>(buffer, size);
+}
+
+template <bool IS_MAX>
+size_t extremum_element(const void *buffer, size_t nElts, GDALDataType eDT,
+                        bool bHasNoData, double dfNoDataValue)
+{
+    switch (eDT)
+    {
+        case GDT_Int8:
+        {
+            using T = int8_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Byte:
+        {
+            using T = uint8_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Int16:
+        {
+            using T = int16_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_UInt16:
+        {
+            using T = uint16_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Int32:
+        {
+            using T = int32_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_UInt32:
+        {
+            using T = uint32_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Int64:
+        {
+            using T = int64_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_UInt64:
+        {
+            using T = uint64_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Float32:
+        {
+            using T = float;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Float64:
+        {
+            using T = double;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return extremum_element<T, IS_MAX>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        default:
+            break;
+    }
+    CPLError(CE_Failure, CPLE_NotSupported,
+             "%s not supported for this data type.", __FUNCTION__);
+    return 0;
+}
+
+}  // namespace detail
+
+/************************************************************************/
+/*                            max_element()                             */
+/************************************************************************/
+
+/** Return the index of the element where the maximum value is hit.
+ *
+ * If it is hit in several locations, it is not specified which one will be
+ * returned.
+ *
+ * @param buffer Vector of nElts elements of type eDT.
+ * @param nElts Number of elements in buffer.
+ * @param eDT Data type of the elements of buffer.
+ * @param bHasNoData Whether dfNoDataValue is valid.
+ * @param dfNoDataValue Nodata value, only taken into account if bHasNoData == true
+ *
+ * @since GDAL 3.11
+ */
+inline size_t max_element(const void *buffer, size_t nElts, GDALDataType eDT,
+                          bool bHasNoData, double dfNoDataValue)
+{
+    return detail::extremum_element<true>(buffer, nElts, eDT, bHasNoData,
+                                          dfNoDataValue);
+}
+
+/************************************************************************/
+/*                            min_element()                             */
+/************************************************************************/
+
+/** Return the index of the element where the minimum value is hit.
+ *
+ * If it is hit in several locations, it is not specified which one will be
+ * returned.
+ *
+ * @param buffer Vector of nElts elements of type eDT.
+ * @param nElts Number of elements in buffer.
+ * @param eDT Data type of the elements of buffer.
+ * @param bHasNoData Whether dfNoDataValue is valid.
+ * @param dfNoDataValue Nodata value, only taken into account if bHasNoData == true
+ *
+ * @since GDAL 3.11
+ */
+inline size_t min_element(const void *buffer, size_t nElts, GDALDataType eDT,
+                          bool bHasNoData, double dfNoDataValue)
+{
+    return detail::extremum_element<false>(buffer, nElts, eDT, bHasNoData,
+                                           dfNoDataValue);
+}
+
+namespace detail
+{
+
+#ifdef NOT_EFFICIENT
+
+/************************************************************************/
+/*                         minmax_element()                             */
+/************************************************************************/
+
+template <class T>
+std::pair<size_t, size_t> minmax_element(const T *v, size_t size, T noDataValue)
+{
+    static_assert(
+        !(std::is_same<T, float>::value || std::is_same<T, double>::value));
+    if (size == 0)
+        return std::pair(0, 0);
+    size_t idx_of_min = 0;
+    size_t idx_of_max = 0;
+    T vmin = v[0];
+    T vmax = v[0];
+    bool extremum_is_nodata = vmin == noDataValue;
+    size_t i = 1;
+    for (; i < size; ++i)
+    {
+        if (v[i] != noDataValue && (v[i] < vmin || extremum_is_nodata))
+        {
+            vmin = v[i];
+            idx_of_min = i;
+            extremum_is_nodata = false;
+        }
+        if (v[i] != noDataValue && (v[i] > vmax || extremum_is_nodata))
+        {
+            vmax = v[i];
+            idx_of_max = i;
+            extremum_is_nodata = false;
+        }
+    }
+    return std::pair(idx_of_min, idx_of_max);
+}
+
+template <class T>
+std::pair<size_t, size_t> minmax_element(const T *v, size_t size)
+{
+    static_assert(
+        !(std::is_same<T, float>::value || std::is_same<T, double>::value));
+    if (size == 0)
+        return std::pair(0, 0);
+    size_t idx_of_min = 0;
+    size_t idx_of_max = 0;
+    T vmin = v[0];
+    T vmax = v[0];
+    size_t i = 1;
+    for (; i < size; ++i)
+    {
+        if (v[i] < vmin)
+        {
+            vmin = v[i];
+            idx_of_min = i;
+        }
+        if (v[i] > vmax)
+        {
+            vmax = v[i];
+            idx_of_max = i;
+        }
+    }
+    return std::pair(idx_of_min, idx_of_max);
+}
+
+template <class T>
+inline std::pair<size_t, size_t> minmax_element_with_nan(const T *v,
+                                                         size_t size)
+{
+    if (size == 0)
+        return std::pair(0, 0);
+    size_t idx_of_min = 0;
+    size_t idx_of_max = 0;
+    T vmin = v[0];
+    T vmax = v[0];
+    bool extremum_is_nan = std::isnan(v[0]);
+    size_t i = 1;
+    if (extremum_is_nan)
+    {
+        for (; i < size; ++i)
+        {
+            if (!std::isnan(v[i]))
+            {
+                vmin = v[i];
+                idx_of_min = i;
+                vmax = v[i];
+                idx_of_max = i;
+                break;
+            }
+        }
+    }
+    for (; i < size; ++i)
+    {
+        {
+            if (v[i] < vmin)
+            {
+                vmin = v[i];
+                idx_of_min = i;
+            }
+            if (v[i] > vmax)
+            {
+                vmax = v[i];
+                idx_of_max = i;
+            }
+        }
+    }
+    return std::pair(idx_of_min, idx_of_max);
+}
+
+template <>
+std::pair<size_t, size_t> minmax_element<float>(const float *v, size_t size)
+{
+    return minmax_element_with_nan<float>(v, size);
+}
+
+template <>
+std::pair<size_t, size_t> minmax_element<double>(const double *v, size_t size)
+{
+    return minmax_element_with_nan<double>(v, size);
+}
+
+template <class T>
+inline std::pair<size_t, size_t> minmax_element(const T *buffer, size_t size,
+                                                bool bHasNoData, T noDataValue)
+{
+    if (bHasNoData)
+    {
+        return minmax_element<T>(buffer, size, noDataValue);
+    }
+    else
+    {
+        return minmax_element<T>(buffer, size);
+    }
+}
+#else
+
+/************************************************************************/
+/*                         minmax_element()                             */
+/************************************************************************/
+
+template <class T>
+inline std::pair<size_t, size_t> minmax_element(const T *buffer, size_t size,
+                                                bool bHasNoData, T noDataValue)
+{
+#ifdef NOT_EFFICIENT
+    if (bHasNoData)
+    {
+        return minmax_element<T>(buffer, size, noDataValue);
+    }
+    else
+    {
+        return minmax_element<T>(buffer, size);
+        //auto [imin, imax] = std::minmax_element(buffer, buffer + size);
+        //return std::pair(imin - buffer, imax - buffer);
+    }
+#else
+    // Using separately min and max is more efficient than computing them
+    // within the same loop
+    return std::pair(
+        extremum_element<T, false>(buffer, size, bHasNoData, noDataValue),
+        extremum_element<T, true>(buffer, size, bHasNoData, noDataValue));
+#endif
+}
+#endif
+
+}  // namespace detail
+
+/************************************************************************/
+/*                          minmax_element()                            */
+/************************************************************************/
+
+/** Return the index of the elements where the minimum and maximum values are hit.
+ *
+ * If they are hit in several locations, it is not specified which one will be
+ * returned (contrary to std::minmax_element).
+ *
+ * @param buffer Vector of nElts elements of type eDT.
+ * @param nElts Number of elements in buffer.
+ * @param eDT Data type of the elements of buffer.
+ * @param bHasNoData Whether dfNoDataValue is valid.
+ * @param dfNoDataValue Nodata value, only taken into account if bHasNoData == true
+ *
+ * @since GDAL 3.11
+ */
+inline std::pair<size_t, size_t> minmax_element(const void *buffer,
+                                                size_t nElts, GDALDataType eDT,
+                                                bool bHasNoData,
+                                                double dfNoDataValue)
+{
+    switch (eDT)
+    {
+        case GDT_Int8:
+        {
+            using T = int8_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Byte:
+        {
+            using T = uint8_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Int16:
+        {
+            using T = int16_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_UInt16:
+        {
+            using T = uint16_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Int32:
+        {
+            using T = int32_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_UInt32:
+        {
+            using T = uint32_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Int64:
+        {
+            using T = int64_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_UInt64:
+        {
+            using T = uint64_t;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Float32:
+        {
+            using T = float;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        case GDT_Float64:
+        {
+            using T = double;
+            bHasNoData = bHasNoData && GDALIsValueExactAs<T>(dfNoDataValue);
+            return detail::minmax_element<T>(
+                static_cast<const T *>(buffer), nElts, bHasNoData,
+                bHasNoData ? static_cast<T>(dfNoDataValue) : 0);
+        }
+        default:
+            break;
+    }
+    CPLError(CE_Failure, CPLE_NotSupported,
+             "%s not supported for this data type.", __FUNCTION__);
+    return std::pair(0, 0);
+}
+
+}  // namespace GDAL_MINMAXELT_NS
+
+#endif  // GDAL_MINMAX_ELEMENT_INCLUDED
diff --git a/perftests/CMakeLists.txt b/perftests/CMakeLists.txt
index 4f365f34a326..e020f471d2b9 100644
--- a/perftests/CMakeLists.txt
+++ b/perftests/CMakeLists.txt
@@ -10,3 +10,7 @@ target_link_libraries(bench_ogr_batch PRIVATE $<TARGET_NAME:${GDAL_LIB_TARGET_NA
 add_executable(bench_ogr_c_api bench_ogr_c_api.cpp)
 gdal_standard_includes(bench_ogr_c_api)
 target_link_libraries(bench_ogr_c_api PRIVATE $<TARGET_NAME:${GDAL_LIB_TARGET_NAME}>)
+
+add_executable(testperf_gdal_minmax_element testperf_gdal_minmax_element.cpp)
+gdal_standard_includes(testperf_gdal_minmax_element)
+target_link_libraries(testperf_gdal_minmax_element PRIVATE $<TARGET_NAME:${GDAL_LIB_TARGET_NAME}>)
diff --git a/perftests/testperf_gdal_minmax_element.cpp b/perftests/testperf_gdal_minmax_element.cpp
new file mode 100644
index 000000000000..3ac6f01a7151
--- /dev/null
+++ b/perftests/testperf_gdal_minmax_element.cpp
@@ -0,0 +1,363 @@
+/******************************************************************************
+ * Project:  GDAL Core
+ * Purpose:  Test performance of gdal_minmax_element.hpp
+ * Author:   Even Rouault, <even dot rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2024, Even Rouault <even dot rouault at spatialys.com>
+ *
+ * SPDX-License-Identifier: MIT
+ ****************************************************************************/
+
+#include "gdal_minmax_element.hpp"
+
+#include <chrono>
+#include <random>
+
+template <class T> void randomFill(T *v, size_t size, bool withNaN = true)
+{
+    std::random_device rd;
+    std::mt19937 gen{rd()};
+    std::normal_distribution<> dist{127, 30};
+    for (size_t i = 0; i < size; i++)
+    {
+        v[i] = static_cast<T>(dist(gen));
+        if constexpr (std::is_same<T, float>::value ||
+                      std::is_same<T, double>::value)
+        {
+            if (withNaN && (i == 0 || (i > 10 && ((i + 1) % 1024) <= 4)))
+                v[i] = std::numeric_limits<float>::quiet_NaN();
+        }
+    }
+}
+
+int main(int /* argc */, char * /* argv */[])
+{
+    constexpr size_t SIZE = 10 * 1000 * 1000 + 1;
+    {
+        printf("uint8:\n");
+        std::vector<uint8_t> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size());
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Byte, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end()))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (nodata case)\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Byte, true, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("int8\n");
+        std::vector<int8_t> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size());
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Int8, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end()))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (nodata case)\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Int8, true, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("uint16\n");
+        std::vector<uint16_t> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size());
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_UInt16, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end()))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (nodata case)\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_UInt16, true, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("int16\n");
+        std::vector<int16_t> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size());
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Int16, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end()))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (nodata case)\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Int16, true, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("uint32\n");
+        std::vector<uint32_t> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size());
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_UInt32, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end()))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (nodata case)\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_UInt32, true, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("int32\n");
+        std::vector<int32_t> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size());
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Int32, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end()))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (nodata case)\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Int32, true, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("float with NaN\n");
+        std::vector<float> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size());
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Float32, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element with NaN aware "
+                   "comparison)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end(),
+                                                   [](float x, float y) {
+                                                       return std::isnan(y)
+                                                                  ? true
+                                                              : std::isnan(x)
+                                                                  ? false
+                                                                  : x < y;
+                                                   }))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            int idx = static_cast<int>(
+                gdal::min_element(x.data(), x.size(), GDT_Float32, true, 0));
+            printf("min at idx %d = %f (nodata case)\n", idx, x[idx]);
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("float without NaN\n");
+        std::vector<float> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size(), false);
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Float32, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end()))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            int idx = static_cast<int>(
+                gdal::min_element(x.data(), x.size(), GDT_Float32, true, 0));
+            printf("min at idx %d = %f (nodata case)\n", idx, x[idx]);
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("double\n");
+        std::vector<double> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size());
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Float64, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element with NaN aware "
+                   "comparison)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end(),
+                                                   [](double x, double y) {
+                                                       return std::isnan(y)
+                                                                  ? true
+                                                              : std::isnan(x)
+                                                                  ? false
+                                                                  : x < y;
+                                                   }))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            int idx = static_cast<int>(
+                gdal::min_element(x.data(), x.size(), GDT_Float64, true, 0));
+            printf("min at idx %d = %f (nodata case)\n", idx, x[idx]);
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    printf("--------------------\n");
+    {
+        printf("double without NaN\n");
+        std::vector<double> x;
+        x.resize(SIZE);
+        randomFill(x.data(), x.size(), false);
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d\n",
+                   static_cast<int>(gdal::min_element(x.data(), x.size(),
+                                                      GDT_Float64, false, 0)));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            printf("min at idx %d (using std::min_element)\n",
+                   static_cast<int>(std::distance(
+                       x.begin(), std::min_element(x.begin(), x.end()))));
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+        {
+            auto start = std::chrono::steady_clock::now();
+            int idx = static_cast<int>(
+                gdal::min_element(x.data(), x.size(), GDT_Float64, true, 0));
+            printf("min at idx %d = %f (nodata case)\n", idx, x[idx]);
+            auto end = std::chrono::steady_clock::now();
+            printf("-> elapsed=%d\n", static_cast<int>((end - start).count()));
+        }
+    }
+    return 0;
+}