LCOV - code coverage report
Current view: top level - gcore - gdal_minmax_element.hpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 341 342 99.7 %
Date: 2025-09-10 17:48:50 Functions: 374 481 77.8 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  * Project:  GDAL Core
       3             :  * Purpose:  Utility functions to find minimum and maximum values in a buffer
       4             :  * Author:   Even Rouault, <even dot rouault at spatialys.com>
       5             :  *
       6             :  ******************************************************************************
       7             :  * Copyright (c) 2024, Even Rouault <even dot rouault at spatialys.com>
       8             :  *
       9             :  * SPDX-License-Identifier: MIT
      10             :  ****************************************************************************/
      11             : 
      12             : #ifndef GDAL_MINMAX_ELEMENT_INCLUDED
      13             : #define GDAL_MINMAX_ELEMENT_INCLUDED
      14             : 
      15             : // NOTE: This header requires C++17
      16             : 
      17             : // This file may be vendored by other applications than GDAL
      18             : // WARNING: if modifying this file, please also update the upstream GDAL version
      19             : // at https://github.com/OSGeo/gdal/blob/master/gcore/gdal_minmax_element.hpp
      20             : 
      21             : #include <algorithm>
      22             : #include <cassert>
      23             : #include <cmath>
      24             : #include <cstdint>
      25             : #include <limits>
      26             : #include <type_traits>
      27             : #include <utility>
      28             : 
      29             : #include "gdal.h"
      30             : 
      31             : // Just to please cppcheck
      32             : #ifndef GDAL_COMPUTE_VERSION
      33             : #define GDAL_COMPUTE_VERSION(maj, min, rev)                                    \
      34             :     ((maj)*1000000 + (min)*10000 + (rev)*100)
      35             : #endif
      36             : 
      37             : #ifdef GDAL_COMPILATION
      38             : #include "cpl_float.h"
      39             : #define GDAL_MINMAXELT_NS gdal
      40             : #elif !defined(GDAL_MINMAXELT_NS)
      41             : #error "Please define the GDAL_MINMAXELT_NS macro to define the namespace"
      42             : #elif GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
      43             : #include "cpl_float.h"
      44             : #endif
      45             : 
      46             : #ifdef USE_NEON_OPTIMIZATIONS
      47             : #include "include_sse2neon.h"
      48             : #define GDAL_MINMAX_ELEMENT_USE_SSE2
      49             : #else
      50             : #if defined(__x86_64) || defined(_M_X64)
      51             : #define GDAL_MINMAX_ELEMENT_USE_SSE2
      52             : #endif
      53             : #ifdef GDAL_MINMAX_ELEMENT_USE_SSE2
      54             : // SSE2 header
      55             : #include <emmintrin.h>
      56             : #if defined(__SSE4_1__) || defined(__AVX__)
      57             : #include <smmintrin.h>
      58             : #endif
      59             : #endif
      60             : #endif
      61             : 
      62             : #include "gdal_priv_templates.hpp"
      63             : 
      64             : namespace GDAL_MINMAXELT_NS
      65             : {
      66             : namespace detail
      67             : {
      68             : 
      69             : template <class T> struct is_floating_point
      70             : {
      71             :     static constexpr bool value = std::is_floating_point_v<T>
      72             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
      73             :                                   || std::is_same_v<T, GFloat16>
      74             : #endif
      75             :         ;
      76             : };
      77             : 
      78             : template <class T>
      79             : constexpr bool is_floating_point_v = is_floating_point<T>::value;
      80             : 
      81             : /************************************************************************/
      82             : /*                            compScalar()                              */
      83             : /************************************************************************/
      84             : 
      85       15849 : template <class T, bool IS_MAX> inline static bool compScalar(T x, T y)
      86             : {
      87             :     if constexpr (IS_MAX)
      88        8355 :         return x > y;
      89             :     else
      90        7494 :         return x < y;
      91             : }
      92             : 
      93        1390 : template <typename T> inline static bool IsNan(T x)
      94             : {
      95             :     // We need to write `using std::isnan` instead of directly using
      96             :     // `std::isnan` because `std::isnan` only supports the types
      97             :     // `float` and `double`. The `isnan` for `cpl::Float16` is found in the
      98             :     // `cpl` namespace via argument-dependent lookup
      99             :     // <https://en.cppreference.com/w/cpp/language/adl>.
     100             :     using std::isnan;
     101        1390 :     return isnan(x);
     102             : }
     103             : 
     104        6557 : template <class T> inline static bool compEqual(T x, T y)
     105             : {
     106        6557 :     return x == y;
     107             : }
     108             : 
     109             : // On Intel/Neon, we do comparisons on uint16_t instead of casting to float for
     110             : // faster execution.
     111             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0) &&                      \
     112             :     defined(GDAL_MINMAX_ELEMENT_USE_SSE2)
     113       11946 : template <> bool IsNan<GFloat16>(GFloat16 x)
     114             : {
     115             :     uint16_t iX;
     116       11946 :     memcpy(&iX, &x, sizeof(x));
     117             :     // Check that the 5 bits of the exponent are set and the mantissa is not zero.
     118       11946 :     return (iX & 0x7fff) > 0x7c00;
     119             : }
     120             : 
     121        2082 : template <> bool compEqual<GFloat16>(GFloat16 x, GFloat16 y)
     122             : {
     123             :     uint16_t iX, iY;
     124        2082 :     memcpy(&iX, &x, sizeof(x));
     125        2082 :     memcpy(&iY, &y, sizeof(y));
     126             :     // Given our usage where y cannot be NaN we can skip the IsNan tests
     127        2082 :     assert(!IsNan(y));
     128        4156 :     return (iX == iY ||
     129             :             // Also check +0 == -0
     130        4156 :             ((iX | iY) == (1 << 15)));
     131             : }
     132             : 
     133        4920 : template <> bool compScalar<GFloat16, true>(GFloat16 x, GFloat16 y)
     134             : {
     135             :     uint16_t iX, iY;
     136        4920 :     memcpy(&iX, &x, sizeof(x));
     137        4920 :     memcpy(&iY, &y, sizeof(y));
     138             :     bool ret;
     139        4920 :     if (IsNan(x) || IsNan(y))
     140             :     {
     141         807 :         ret = false;
     142             :     }
     143        4113 :     else if (!(iX >> 15))
     144             :     {
     145             :         // +0 considered > -0. We don't really care
     146         173 :         ret = (iY >> 15) || iX > iY;
     147             :     }
     148             :     else
     149             :     {
     150        3940 :         ret = (iY >> 15) && iX < iY;
     151             :     }
     152        4920 :     return ret;
     153             : }
     154             : 
     155        2252 : template <> bool compScalar<GFloat16, false>(GFloat16 x, GFloat16 y)
     156             : {
     157        2252 :     return compScalar<GFloat16, true>(y, x);
     158             : }
     159             : #endif
     160             : 
     161             : /************************************************************************/
     162             : /*                 extremum_element_with_nan_generic()                  */
     163             : /************************************************************************/
     164             : 
     165             : template <class T, bool IS_MAX>
     166             : inline size_t extremum_element_with_nan_generic(const T *v, size_t size)
     167             : {
     168             :     if (size == 0)
     169             :         return 0;
     170             :     size_t idx_of_extremum = 0;
     171             :     auto extremum = v[0];
     172             :     bool extremum_is_nan = IsNan(extremum);
     173             :     size_t i = 1;
     174             :     for (; i < size; ++i)
     175             :     {
     176             :         if (compScalar<T, IS_MAX>(v[i], extremum) ||
     177             :             (extremum_is_nan && !IsNan(v[i])))
     178             :         {
     179             :             extremum = v[i];
     180             :             idx_of_extremum = i;
     181             :             extremum_is_nan = false;
     182             :         }
     183             :     }
     184             :     return idx_of_extremum;
     185             : }
     186             : 
     187             : /************************************************************************/
     188             : /*                 extremum_element_with_nan_generic()                  */
     189             : /************************************************************************/
     190             : 
     191             : template <class T, bool IS_MAX>
     192             : inline size_t extremum_element_with_nan_generic(const T *v, size_t size,
     193             :                                                 T noDataValue)
     194             : {
     195             :     if (IsNan(noDataValue))
     196             :         return extremum_element_with_nan_generic<T, IS_MAX>(v, size);
     197             :     if (size == 0)
     198             :         return 0;
     199             :     size_t idx_of_extremum = 0;
     200             :     auto extremum = v[0];
     201             :     bool extremum_is_nan_or_nodata =
     202             :         IsNan(extremum) || compEqual(extremum, noDataValue);
     203             :     size_t i = 1;
     204             :     for (; i < size; ++i)
     205             :     {
     206             :         if (!compEqual(v[i], noDataValue) &&
     207             :             (compScalar<T, IS_MAX>(v[i], extremum) ||
     208             :              (extremum_is_nan_or_nodata && !IsNan(v[i]))))
     209             :         {
     210             :             extremum = v[i];
     211             :             idx_of_extremum = i;
     212             :             extremum_is_nan_or_nodata = false;
     213             :         }
     214             :     }
     215             :     return idx_of_extremum;
     216             : }
     217             : 
     218             : /************************************************************************/
     219             : /*                       extremum_element_generic()                     */
     220             : /************************************************************************/
     221             : 
     222             : template <class T, bool IS_MAX>
     223             : inline size_t extremum_element_generic(const T *buffer, size_t size,
     224             :                                        bool bHasNoData, T noDataValue)
     225             : {
     226             :     if (bHasNoData)
     227             :     {
     228             :         if constexpr (is_floating_point_v<T>)
     229             :         {
     230             :             if (IsNan(noDataValue))
     231             :             {
     232             :                 if constexpr (IS_MAX)
     233             :                 {
     234             :                     return std::max_element(buffer, buffer + size,
     235             :                                             [](T a, T b) {
     236             :                                                 return IsNan(b)   ? false
     237             :                                                        : IsNan(a) ? true
     238             :                                                                   : a < b;
     239             :                                             }) -
     240             :                            buffer;
     241             :                 }
     242             :                 else
     243             :                 {
     244             :                     return std::min_element(buffer, buffer + size,
     245             :                                             [](T a, T b) {
     246             :                                                 return IsNan(b)   ? true
     247             :                                                        : IsNan(a) ? false
     248             :                                                                   : a < b;
     249             :                                             }) -
     250             :                            buffer;
     251             :                 }
     252             :             }
     253             :             else
     254             :             {
     255             :                 if constexpr (IS_MAX)
     256             :                 {
     257             :                     return std::max_element(buffer, buffer + size,
     258             :                                             [noDataValue](T a, T b)
     259             :                                             {
     260             :                                                 return IsNan(b)   ? false
     261             :                                                        : IsNan(a) ? true
     262             :                                                        : (b == noDataValue)
     263             :                                                            ? false
     264             :                                                        : (a == noDataValue)
     265             :                                                            ? true
     266             :                                                            : a < b;
     267             :                                             }) -
     268             :                            buffer;
     269             :                 }
     270             :                 else
     271             :                 {
     272             :                     return std::min_element(buffer, buffer + size,
     273             :                                             [noDataValue](T a, T b)
     274             :                                             {
     275             :                                                 return IsNan(b)   ? true
     276             :                                                        : IsNan(a) ? false
     277             :                                                        : (b == noDataValue)
     278             :                                                            ? true
     279             :                                                        : (a == noDataValue)
     280             :                                                            ? false
     281             :                                                            : a < b;
     282             :                                             }) -
     283             :                            buffer;
     284             :                 }
     285             :             }
     286             :         }
     287             :         else
     288             :         {
     289             :             if constexpr (IS_MAX)
     290             :             {
     291             :                 return std::max_element(buffer, buffer + size,
     292             :                                         [noDataValue](T a, T b) {
     293             :                                             return (b == noDataValue)   ? false
     294             :                                                    : (a == noDataValue) ? true
     295             :                                                                         : a < b;
     296             :                                         }) -
     297             :                        buffer;
     298             :             }
     299             :             else
     300             :             {
     301             :                 return std::min_element(buffer, buffer + size,
     302             :                                         [noDataValue](T a, T b) {
     303             :                                             return (b == noDataValue)   ? true
     304             :                                                    : (a == noDataValue) ? false
     305             :                                                                         : a < b;
     306             :                                         }) -
     307             :                        buffer;
     308             :             }
     309             :         }
     310             :     }
     311             :     else
     312             :     {
     313             :         if constexpr (is_floating_point_v<T>)
     314             :         {
     315             :             if constexpr (IS_MAX)
     316             :             {
     317             :                 return std::max_element(buffer, buffer + size,
     318             :                                         [](T a, T b) {
     319             :                                             return IsNan(b)   ? false
     320             :                                                    : IsNan(a) ? true
     321             :                                                               : a < b;
     322             :                                         }) -
     323             :                        buffer;
     324             :             }
     325             :             else
     326             :             {
     327             :                 return std::min_element(buffer, buffer + size,
     328             :                                         [](T a, T b) {
     329             :                                             return IsNan(b)   ? true
     330             :                                                    : IsNan(a) ? false
     331             :                                                               : a < b;
     332             :                                         }) -
     333             :                        buffer;
     334             :             }
     335             :         }
     336             :         else
     337             :         {
     338             :             if constexpr (IS_MAX)
     339             :             {
     340             :                 return std::max_element(buffer, buffer + size) - buffer;
     341             :             }
     342             :             else
     343             :             {
     344             :                 return std::min_element(buffer, buffer + size) - buffer;
     345             :             }
     346             :         }
     347             :     }
     348             : }
     349             : 
     350             : #ifdef GDAL_MINMAX_ELEMENT_USE_SSE2
     351             : 
     352             : /************************************************************************/
     353             : /*                     extremum_element_sse2()                          */
     354             : /************************************************************************/
     355             : 
     356         104 : static inline int8_t Shift8(uint8_t x)
     357             : {
     358         104 :     return static_cast<int8_t>(x + std::numeric_limits<int8_t>::min());
     359             : }
     360             : 
     361          42 : static inline int16_t Shift16(uint16_t x)
     362             : {
     363          42 :     return static_cast<int16_t>(x + std::numeric_limits<int16_t>::min());
     364             : }
     365             : 
     366             : CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW
     367          53 : static inline int32_t Shift32(uint32_t x)
     368             : {
     369          53 :     x += static_cast<uint32_t>(std::numeric_limits<int32_t>::min());
     370             :     int32_t ret;
     371          53 :     memcpy(&ret, &x, sizeof(x));
     372          53 :     return ret;
     373             : }
     374             : 
     375             : CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW
     376          59 : static inline int64_t Shift64(uint64_t x)
     377             : {
     378          59 :     x += static_cast<uint64_t>(std::numeric_limits<int64_t>::min());
     379             :     int64_t ret;
     380          59 :     memcpy(&ret, &x, sizeof(x));
     381          59 :     return ret;
     382             : }
     383             : 
     384             : // Return a _mm128[i|d] register with all its elements set to x
     385         897 : template <class T> static inline auto set1(T x)
     386             : {
     387             :     if constexpr (std::is_same_v<T, uint8_t>)
     388         208 :         return _mm_set1_epi8(Shift8(x));
     389             :     else if constexpr (std::is_same_v<T, int8_t>)
     390          42 :         return _mm_set1_epi8(x);
     391             :     else if constexpr (std::is_same_v<T, uint16_t>)
     392          84 :         return _mm_set1_epi16(Shift16(x));
     393             :     else if constexpr (std::is_same_v<T, int16_t>)
     394         102 :         return _mm_set1_epi16(x);
     395             :     else if constexpr (std::is_same_v<T, uint32_t>)
     396         106 :         return _mm_set1_epi32(Shift32(x));
     397             :     else if constexpr (std::is_same_v<T, int32_t>)
     398          74 :         return _mm_set1_epi32(x);
     399             :     else if constexpr (std::is_same_v<T, uint64_t>)
     400         118 :         return _mm_set1_epi64x(Shift64(x));
     401             :     else if constexpr (std::is_same_v<T, int64_t>)
     402          68 :         return _mm_set1_epi64x(x);
     403             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
     404             :     else if constexpr (std::is_same_v<T, GFloat16>)
     405             :     {
     406             :         int16_t iX;
     407         147 :         memcpy(&iX, &x, sizeof(x));
     408         294 :         return _mm_set1_epi16(iX);
     409             :     }
     410             : #endif
     411             :     else if constexpr (std::is_same_v<T, float>)
     412         134 :         return _mm_set1_ps(x);
     413             :     else
     414         144 :         return _mm_set1_pd(x);
     415             : }
     416             : 
     417             : // Return a _mm128[i|d] register with all its elements set to x
     418         944 : template <class T> static inline auto set1_unshifted(T x)
     419             : {
     420             :     if constexpr (std::is_same_v<T, uint8_t>)
     421             :     {
     422             :         int8_t xSigned;
     423         118 :         memcpy(&xSigned, &x, sizeof(xSigned));
     424         236 :         return _mm_set1_epi8(xSigned);
     425             :     }
     426             :     else if constexpr (std::is_same_v<T, int8_t>)
     427          84 :         return _mm_set1_epi8(x);
     428             :     else if constexpr (std::is_same_v<T, uint16_t>)
     429             :     {
     430             :         int16_t xSigned;
     431          52 :         memcpy(&xSigned, &x, sizeof(xSigned));
     432         104 :         return _mm_set1_epi16(xSigned);
     433             :     }
     434             :     else if constexpr (std::is_same_v<T, int16_t>)
     435         114 :         return _mm_set1_epi16(x);
     436             :     else if constexpr (std::is_same_v<T, uint32_t>)
     437             :     {
     438             :         int32_t xSigned;
     439          58 :         memcpy(&xSigned, &x, sizeof(xSigned));
     440         116 :         return _mm_set1_epi32(xSigned);
     441             :     }
     442             :     else if constexpr (std::is_same_v<T, int32_t>)
     443          68 :         return _mm_set1_epi32(x);
     444             :     else if constexpr (std::is_same_v<T, uint64_t>)
     445             :     {
     446             :         int64_t xSigned;
     447          60 :         memcpy(&xSigned, &x, sizeof(xSigned));
     448         120 :         return _mm_set1_epi64x(xSigned);
     449             :     }
     450             :     else if constexpr (std::is_same_v<T, int64_t>)
     451          65 :         return _mm_set1_epi64x(x);
     452             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
     453             :     else if constexpr (std::is_same_v<T, GFloat16>)
     454             :     {
     455             :         int16_t iX;
     456         148 :         memcpy(&iX, &x, sizeof(x));
     457         296 :         return _mm_set1_epi16(iX);
     458             :     }
     459             : #endif
     460             :     else if constexpr (std::is_same_v<T, float>)
     461         139 :         return _mm_set1_ps(x);
     462             :     else
     463         137 :         return _mm_set1_pd(x);
     464             : }
     465             : 
     466             : // Load as many values of type T at a _mm128[i|d] register can contain from x
     467   145001866 : template <class T> static inline auto loadv(const T *x)
     468             : {
     469             :     if constexpr (std::is_same_v<T, float>)
     470    20000132 :         return _mm_loadu_ps(x);
     471             :     else if constexpr (std::is_same_v<T, double>)
     472    40000148 :         return _mm_loadu_pd(x);
     473             :     else
     474    85001586 :         return _mm_loadu_si128(reinterpret_cast<const __m128i *>(x));
     475             : }
     476             : 
     477    15000144 : inline __m128i IsNanGFloat16(__m128i x)
     478             : {
     479             :     // (iX & 0x7fff) > 0x7c00
     480    45000332 :     return _mm_cmpgt_epi16(_mm_and_si128(x, _mm_set1_epi16(0x7fff)),
     481    15000144 :                            _mm_set1_epi16(0x7c00));
     482             : }
     483             : 
     484             : template <class T, class SSE_T>
     485    82500682 : static inline SSE_T blendv(SSE_T a, SSE_T b, SSE_T mask)
     486             : {
     487             :     if constexpr (std::is_same_v<T, float>)
     488             :     {
     489             : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
     490             :         return _mm_blendv_ps(a, b, mask);
     491             : #else
     492    30000180 :         return _mm_or_ps(_mm_andnot_ps(mask, a), _mm_and_ps(mask, b));
     493             : #endif
     494             :     }
     495             :     else if constexpr (std::is_same_v<T, double>)
     496             :     {
     497             : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
     498             :         return _mm_blendv_pd(a, b, mask);
     499             : #else
     500    60000372 :         return _mm_or_pd(_mm_andnot_pd(mask, a), _mm_and_pd(mask, b));
     501             : #endif
     502             :     }
     503             :     else
     504             :     {
     505             : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
     506             :         return _mm_blendv_epi8(a, b, mask);
     507             : #else
     508   157501664 :         return _mm_or_si128(_mm_andnot_si128(mask, a), _mm_and_si128(mask, b));
     509             : #endif
     510             :     }
     511             : }
     512             : 
     513    10000058 : inline __m128i cmpgt_ph(__m128i x, __m128i y)
     514             : {
     515             : #ifdef slow
     516             :     GFloat16 vx[8], vy[8];
     517             :     int16_t res[8];
     518             :     _mm_storeu_si128(reinterpret_cast<__m128i *>(vx), x);
     519             :     _mm_storeu_si128(reinterpret_cast<__m128i *>(vy), y);
     520             :     for (int i = 0; i < 8; ++i)
     521             :     {
     522             :         res[i] = vx[i] > vy[i] ? -1 : 0;
     523             :     }
     524             :     return _mm_loadu_si128(reinterpret_cast<const __m128i *>(res));
     525             : #else
     526    10000058 :     const auto x_is_negative = _mm_srai_epi16(x, 15);
     527    10000058 :     const auto y_is_negative = _mm_srai_epi16(y, 15);
     528    40000252 :     return _mm_andnot_si128(
     529             :         // only x can be NaN given how we use this method
     530             :         IsNanGFloat16(x),
     531             :         blendv<int16_t>(_mm_or_si128(y_is_negative, _mm_cmpgt_epi16(x, y)),
     532             :                         _mm_and_si128(y_is_negative, _mm_cmpgt_epi16(y, x)),
     533    10000058 :                         x_is_negative));
     534             : #endif
     535             : }
     536             : 
     537    40000744 : inline __m128i cmpgt_epi64(__m128i x, __m128i y)
     538             : {
     539             : #if defined(__SSE4_2__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
     540             :     return _mm_cmpgt_epi64(x, y);
     541             : #else
     542   120002232 :     auto tmp = _mm_and_si128(_mm_sub_epi64(y, x), _mm_cmpeq_epi32(x, y));
     543    40000744 :     tmp = _mm_or_si128(tmp, _mm_cmpgt_epi32(x, y));
     544             :     // Replicate the 2 odd-indexed (hi) 32-bit word into the 2 even-indexed (lo)
     545             :     // ones
     546    40000744 :     return _mm_shuffle_epi32(tmp, _MM_SHUFFLE(3, 3, 1, 1));
     547             : #endif
     548             : }
     549             : 
     550             : // Return a __m128i register with bits set when x[i] < y[i] when !IS_MAX
     551             : // or x[i] > y[i] when IS_MAX
     552             : template <class T, bool IS_MAX, class SSE_T>
     553   145001776 : static inline __m128i comp(SSE_T x, SSE_T y)
     554             : {
     555             :     if constexpr (IS_MAX)
     556             :     {
     557             :         if constexpr (std::is_same_v<T, uint8_t>)
     558     3750472 :             return _mm_cmpgt_epi8(
     559             :                 _mm_add_epi8(x,
     560     1250154 :                              _mm_set1_epi8(std::numeric_limits<int8_t>::min())),
     561     1250154 :                 y);
     562             :         else if constexpr (std::is_same_v<T, int8_t>)
     563     1250002 :             return _mm_cmpgt_epi8(x, y);
     564             :         else if constexpr (std::is_same_v<T, uint16_t>)
     565     7500064 :             return _mm_cmpgt_epi16(
     566             :                 _mm_add_epi16(
     567     2500018 :                     x, _mm_set1_epi16(std::numeric_limits<int16_t>::min())),
     568     2500018 :                 y);
     569             :         else if constexpr (std::is_same_v<T, int16_t>)
     570     2500018 :             return _mm_cmpgt_epi16(x, y);
     571             :         else if constexpr (std::is_same_v<T, uint32_t>)
     572    15000180 :             return _mm_cmpgt_epi32(
     573             :                 _mm_add_epi32(
     574             :                     x, _mm_set1_epi32(std::numeric_limits<int32_t>::min())),
     575     5000050 :                 y);
     576             :         else if constexpr (std::is_same_v<T, int32_t>)
     577     5000050 :             return _mm_cmpgt_epi32(x, y);
     578             :         else if constexpr (std::is_same_v<T, uint64_t>)
     579             :         {
     580    20000248 :             return cmpgt_epi64(
     581             :                 _mm_add_epi64(
     582    10000114 :                     x, _mm_set1_epi64x(std::numeric_limits<int64_t>::min())),
     583    10000114 :                 y);
     584             :         }
     585             :         else if constexpr (std::is_same_v<T, int64_t>)
     586             :         {
     587    10000114 :             return cmpgt_epi64(x, y);
     588             :         }
     589             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
     590             :         else if constexpr (std::is_same_v<T, GFloat16>)
     591             :         {
     592     5000024 :             return cmpgt_ph(x, y);
     593             :         }
     594             : #endif
     595             :         else if constexpr (std::is_same_v<T, float>)
     596    20000072 :             return _mm_castps_si128(_mm_cmpgt_ps(x, y));
     597             :         else
     598    40000024 :             return _mm_castpd_si128(_mm_cmpgt_pd(x, y));
     599             :     }
     600             :     else
     601             :     {
     602             :         if constexpr (std::is_same_v<T, uint8_t>)
     603     3750508 :             return _mm_cmplt_epi8(
     604             :                 _mm_add_epi8(x,
     605     1250166 :                              _mm_set1_epi8(std::numeric_limits<int8_t>::min())),
     606     1250166 :                 y);
     607             :         else if constexpr (std::is_same_v<T, int8_t>)
     608     1250014 :             return _mm_cmplt_epi8(x, y);
     609             :         else if constexpr (std::is_same_v<T, uint16_t>)
     610     7500148 :             return _mm_cmplt_epi16(
     611             :                 _mm_add_epi16(
     612     2500046 :                     x, _mm_set1_epi16(std::numeric_limits<int16_t>::min())),
     613     2500046 :                 y);
     614             :         else if constexpr (std::is_same_v<T, int16_t>)
     615     2500046 :             return _mm_cmplt_epi16(x, y);
     616             :         else if constexpr (std::is_same_v<T, uint32_t>)
     617    15000360 :             return _mm_cmplt_epi32(
     618             :                 _mm_add_epi32(
     619             :                     x, _mm_set1_epi32(std::numeric_limits<int32_t>::min())),
     620     5000110 :                 y);
     621             :         else if constexpr (std::is_same_v<T, int32_t>)
     622     5000110 :             return _mm_cmplt_epi32(x, y);
     623             :         else if constexpr (std::is_same_v<T, uint64_t>)
     624             :         {
     625    20000496 :             return cmpgt_epi64(
     626             :                 y, _mm_add_epi64(x, _mm_set1_epi64x(
     627    20000496 :                                         std::numeric_limits<int64_t>::min())));
     628             :         }
     629             :         else if constexpr (std::is_same_v<T, int64_t>)
     630             :         {
     631    10000238 :             return cmpgt_epi64(y, x);
     632             :         }
     633             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
     634             :         else if constexpr (std::is_same_v<T, GFloat16>)
     635             :         {
     636     5000024 :             return cmpgt_ph(y, x);
     637             :         }
     638             : #endif
     639             :         else if constexpr (std::is_same_v<T, float>)
     640    20000192 :             return _mm_castps_si128(_mm_cmplt_ps(x, y));
     641             :         else
     642    40000272 :             return _mm_castpd_si128(_mm_cmplt_pd(x, y));
     643             :     }
     644             : }
     645             : 
     646             : template <class T, class SSE_T> static inline SSE_T compeq(SSE_T a, SSE_T b);
     647             : 
     648     1250018 : template <> __m128i compeq<uint8_t>(__m128i a, __m128i b)
     649             : {
     650     1250018 :     return _mm_cmpeq_epi8(a, b);
     651             : }
     652             : 
     653     1250002 : template <> __m128i compeq<int8_t>(__m128i a, __m128i b)
     654             : {
     655     1250002 :     return _mm_cmpeq_epi8(a, b);
     656             : }
     657             : 
     658     2500018 : template <> __m128i compeq<uint16_t>(__m128i a, __m128i b)
     659             : {
     660     2500018 :     return _mm_cmpeq_epi16(a, b);
     661             : }
     662             : 
     663     2500018 : template <> __m128i compeq<int16_t>(__m128i a, __m128i b)
     664             : {
     665     2500018 :     return _mm_cmpeq_epi16(a, b);
     666             : }
     667             : 
     668     5000050 : template <> __m128i compeq<uint32_t>(__m128i a, __m128i b)
     669             : {
     670     5000050 :     return _mm_cmpeq_epi32(a, b);
     671             : }
     672             : 
     673     5000050 : template <> __m128i compeq<int32_t>(__m128i a, __m128i b)
     674             : {
     675     5000050 :     return _mm_cmpeq_epi32(a, b);
     676             : }
     677             : 
     678    20000248 : template <> __m128i compeq<uint64_t>(__m128i a, __m128i b)
     679             : {
     680             : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
     681             :     return _mm_cmpeq_epi64(a, b);
     682             : #else
     683    20000248 :     auto tmp = _mm_cmpeq_epi32(a, b);
     684             :     // The shuffle swaps hi-lo 32-bit words
     685    40000496 :     return _mm_and_si128(tmp, _mm_shuffle_epi32(tmp, _MM_SHUFFLE(2, 3, 0, 1)));
     686             : #endif
     687             : }
     688             : 
     689    10000114 : template <> __m128i compeq<int64_t>(__m128i a, __m128i b)
     690             : {
     691    10000114 :     return compeq<uint64_t>(a, b);
     692             : }
     693             : 
     694    10000040 : template <> __m128 compeq<float>(__m128 a, __m128 b)
     695             : {
     696    10000040 :     return _mm_cmpeq_ps(a, b);
     697             : }
     698             : 
     699    20000124 : template <> __m128d compeq<double>(__m128d a, __m128d b)
     700             : {
     701    20000124 :     return _mm_cmpeq_pd(a, b);
     702             : }
     703             : 
     704             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
     705     5000036 : template <> __m128i compeq<GFloat16>(__m128i a, __m128i b)
     706             : {
     707             :     // !isnan(a) && !isnan(b) && (a == b || (a|b) == 0x8000)
     708    30000236 :     return _mm_andnot_si128(
     709             :         IsNanGFloat16(a),  // b cannot be NaN given how we use this method
     710             :         _mm_or_si128(_mm_cmpeq_epi16(a, b),
     711             :                      _mm_cmpeq_epi16(
     712             :                          _mm_or_si128(a, b),
     713    10000082 :                          _mm_set1_epi16(std::numeric_limits<int16_t>::min()))));
     714             : }
     715             : #endif
     716             : 
     717             : // Using SSE2
     718             : template <class T, bool IS_MAX, bool HAS_NODATA>
     719             : #if defined(__GNUC__)
     720             : __attribute__((noinline))
     721             : #endif
     722             : size_t
     723         260 : extremum_element_sse2(const T *v, size_t size, T noDataValue)
     724             : {
     725             :     static_assert(std::is_same_v<T, uint8_t> || std::is_same_v<T, int8_t> ||
     726             :                   std::is_same_v<T, uint16_t> || std::is_same_v<T, int16_t> ||
     727             :                   std::is_same_v<T, uint32_t> || std::is_same_v<T, int32_t> ||
     728             :                   std::is_same_v<T, uint64_t> || std::is_same_v<T, int64_t> ||
     729             :                   is_floating_point_v<T>);
     730         260 :     if (size == 0)
     731          22 :         return 0;
     732         238 :     size_t idx_of_extremum = 0;
     733         238 :     T extremum = v[0];
     734         238 :     [[maybe_unused]] bool extremum_is_invalid = false;
     735             :     if constexpr (is_floating_point_v<T>)
     736             :     {
     737             :         if constexpr (HAS_NODATA)
     738             :         {
     739          44 :             if (IsNan(noDataValue))
     740             :             {
     741           6 :                 return extremum_element_sse2<T, IS_MAX, false>(
     742           4 :                     v, size, static_cast<T>(0));
     743             :             }
     744             :         }
     745          86 :         extremum_is_invalid = IsNan(extremum);
     746             :     }
     747             :     if constexpr (HAS_NODATA)
     748             :     {
     749         115 :         if (compEqual(extremum, noDataValue))
     750          25 :             extremum_is_invalid = true;
     751             :     }
     752         234 :     size_t i = 1;
     753             : 
     754         234 :     constexpr size_t VALS_PER_REG = sizeof(set1(extremum)) / sizeof(extremum);
     755         234 :     constexpr int LOOP_UNROLLING = 4;
     756             :     // If changing the value, then we need to adjust the number of sse_valX
     757             :     // loading in the loop.
     758             :     static_assert(LOOP_UNROLLING == 4);
     759         234 :     constexpr size_t VALS_PER_ITER = VALS_PER_REG * LOOP_UNROLLING;
     760             : 
     761       75106 :     const auto update = [v, noDataValue, &extremum, &idx_of_extremum,
     762             :                          &extremum_is_invalid](size_t idx)
     763             :     {
     764             :         if constexpr (HAS_NODATA)
     765             :         {
     766        8524 :             if (compEqual(v[idx], noDataValue))
     767         693 :                 return;
     768        7831 :             if (extremum_is_invalid)
     769             :             {
     770             :                 if constexpr (is_floating_point_v<T>)
     771             :                 {
     772          24 :                     if (IsNan(v[idx]))
     773           6 :                         return;
     774             :                 }
     775          36 :                 extremum = v[idx];
     776          36 :                 idx_of_extremum = idx;
     777          36 :                 extremum_is_invalid = false;
     778          36 :                 return;
     779             :             }
     780             :         }
     781             :         else
     782             :         {
     783       12980 :             CPL_IGNORE_RET_VAL(noDataValue);
     784             :         }
     785       20769 :         if (compScalar<T, IS_MAX>(v[idx], extremum))
     786             :         {
     787         954 :             extremum = v[idx];
     788         954 :             idx_of_extremum = idx;
     789         954 :             extremum_is_invalid = false;
     790             :         }
     791             :         else if constexpr (is_floating_point_v<T>)
     792             :         {
     793        8688 :             if (extremum_is_invalid && !IsNan(v[idx]))
     794             :             {
     795          28 :                 extremum = v[idx];
     796          28 :                 idx_of_extremum = idx;
     797          28 :                 extremum_is_invalid = false;
     798             :             }
     799             :         }
     800             :     };
     801             : 
     802        4381 :     for (; i < VALS_PER_ITER && i < size; ++i)
     803             :     {
     804        4147 :         update(i);
     805             :     }
     806             : 
     807         234 :     [[maybe_unused]] auto sse_neutral = set1_unshifted(static_cast<T>(0));
     808         234 :     [[maybe_unused]] auto sse_nodata = set1_unshifted(noDataValue);
     809             :     if constexpr (HAS_NODATA || is_floating_point_v<T>)
     810             :     {
     811        2240 :         for (; i < size && extremum_is_invalid; ++i)
     812             :         {
     813        2079 :             update(i);
     814             :         }
     815         161 :         if (!extremum_is_invalid)
     816             :         {
     817         234 :             for (; i < size && (i % VALS_PER_ITER) != 0; ++i)
     818             :             {
     819          76 :                 update(i);
     820             :             }
     821         158 :             sse_neutral = set1_unshifted(extremum);
     822             :         }
     823             :     }
     824             : 
     825         234 :     auto sse_extremum = set1(extremum);
     826             : 
     827         234 :     [[maybe_unused]] size_t hits = 0;
     828         234 :     const auto sse_iter_count = (size / VALS_PER_ITER) * VALS_PER_ITER;
     829    36250682 :     for (; i < sse_iter_count; i += VALS_PER_ITER)
     830             :     {
     831             :         // A bit of loop unrolling to save 3/4 of slow movemask operations.
     832    36250472 :         auto sse_val0 = loadv(v + i + 0 * VALS_PER_REG);
     833    36250472 :         auto sse_val1 = loadv(v + i + 1 * VALS_PER_REG);
     834    36250472 :         auto sse_val2 = loadv(v + i + 2 * VALS_PER_REG);
     835    36250472 :         auto sse_val3 = loadv(v + i + 3 * VALS_PER_REG);
     836             : 
     837             :         if constexpr (HAS_NODATA)
     838             :         {
     839             :             // Replace all components that are at the nodata value by a
     840             :             // neutral value (current minimum)
     841    18125160 :             const auto replaceNoDataByNeutral =
     842   145001168 :                 [sse_neutral, sse_nodata](auto sse_val)
     843             :             {
     844    72500604 :                 const auto eq_nodata = compeq<T>(sse_val, sse_nodata);
     845    72500604 :                 return blendv<T>(sse_val, sse_neutral, eq_nodata);
     846             :             };
     847             : 
     848    18125160 :             sse_val0 = replaceNoDataByNeutral(sse_val0);
     849    18125160 :             sse_val1 = replaceNoDataByNeutral(sse_val1);
     850    18125160 :             sse_val2 = replaceNoDataByNeutral(sse_val2);
     851    18125160 :             sse_val3 = replaceNoDataByNeutral(sse_val3);
     852             :         }
     853             : 
     854   145001864 :         if (_mm_movemask_epi8(_mm_or_si128(
     855             :                 _mm_or_si128(comp<T, IS_MAX>(sse_val0, sse_extremum),
     856             :                              comp<T, IS_MAX>(sse_val1, sse_extremum)),
     857             :                 _mm_or_si128(comp<T, IS_MAX>(sse_val2, sse_extremum),
     858    36250472 :                              comp<T, IS_MAX>(sse_val3, sse_extremum)))) != 0)
     859             :         {
     860             :             if constexpr (!std::is_same_v<T, int8_t> &&
     861             :                           !std::is_same_v<T, uint8_t>)
     862             :             {
     863             :                 // The above tests excluding int8_t/uint8_t is due to the fact
     864             :                 // with those small ranges of values we will quickly converge
     865             :                 // to the minimum, so no need to do the below "smart" test.
     866             : 
     867         596 :                 if (++hits == size / 16)
     868             :                 {
     869             :                     // If we have an almost sorted array, then using this code path
     870             :                     // will hurt performance. Arbitrary give up if we get here
     871             :                     // more than 1. / 16 of the size of the array.
     872             :                     // fprintf(stderr, "going to non-vector path\n");
     873           0 :                     break;
     874             :                 }
     875             :             }
     876       15655 :             for (size_t j = 0; j < VALS_PER_ITER; j++)
     877             :             {
     878       14992 :                 update(i + j);
     879             :             }
     880             : 
     881         663 :             sse_extremum = set1(extremum);
     882             :             if constexpr (HAS_NODATA)
     883             :             {
     884         318 :                 sse_neutral = set1_unshifted(extremum);
     885             :             }
     886             :         }
     887             :     }
     888         444 :     for (; i < size; ++i)
     889             :     {
     890         210 :         update(i);
     891             :     }
     892         234 :     return idx_of_extremum;
     893             : }
     894             : 
     895             : /************************************************************************/
     896             : /*                         extremum_element()                           */
     897             : /************************************************************************/
     898             : 
     899             : template <class T, bool IS_MAX>
     900         130 : inline size_t extremum_element(const T *buffer, size_t size, T noDataValue)
     901             : {
     902         130 :     return extremum_element_sse2<T, IS_MAX, true>(buffer, size, noDataValue);
     903             : }
     904             : 
     905             : template <class T, bool IS_MAX>
     906         126 : inline size_t extremum_element(const T *buffer, size_t size)
     907             : {
     908         141 :     return extremum_element_sse2<T, IS_MAX, false>(buffer, size,
     909         141 :                                                    static_cast<T>(0));
     910             : }
     911             : 
     912             : #else
     913             : 
     914             : template <class T, bool IS_MAX>
     915             : inline size_t extremum_element(const T *buffer, size_t size, T noDataValue)
     916             : {
     917             :     if constexpr (is_floating_point_v<T>)
     918             :         return extremum_element_with_nan_generic<T, IS_MAX>(buffer, size,
     919             :                                                             noDataValue);
     920             :     else
     921             :         return extremum_element_generic<T, IS_MAX>(buffer, size, true,
     922             :                                                    noDataValue);
     923             : }
     924             : 
     925             : template <class T, bool IS_MAX>
     926             : inline size_t extremum_element(const T *buffer, size_t size)
     927             : {
     928             :     if constexpr (is_floating_point_v<T>)
     929             :         return extremum_element_with_nan_generic<T, IS_MAX>(buffer, size);
     930             :     else
     931             :         return extremum_element_generic<T, IS_MAX>(buffer, size, false,
     932             :                                                    static_cast<T>(0));
     933             : }
     934             : 
     935             : #endif
     936             : 
     937             : /************************************************************************/
     938             : /*                            extremum_element()                        */
     939             : /************************************************************************/
     940             : 
     941             : template <class T, bool IS_MAX>
     942         256 : inline size_t extremum_element(const T *buffer, size_t size, bool bHasNoData,
     943             :                                T noDataValue)
     944             : {
     945             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0) &&                      \
     946             :     !defined(GDAL_MINMAX_ELEMENT_USE_SSE2)
     947             :     if constexpr (std::is_same_v<T, GFloat16>)
     948             :     {
     949             :         if (bHasNoData)
     950             :             return extremum_element_with_nan_generic<T, IS_MAX>(buffer, size,
     951             :                                                                 noDataValue);
     952             :         else
     953             :             return extremum_element_with_nan_generic<T, IS_MAX>(buffer, size);
     954             :     }
     955             :     else
     956             : #endif
     957         256 :         if (bHasNoData)
     958         130 :         return extremum_element<T, IS_MAX>(buffer, size, noDataValue);
     959             :     else
     960         126 :         return extremum_element<T, IS_MAX>(buffer, size);
     961             : }
     962             : 
     963             : template <class T, class NODATAType>
     964         119 : inline bool IsValueExactAs([[maybe_unused]] NODATAType noDataValue)
     965             : {
     966             :     if constexpr (std::is_same_v<T, NODATAType>)
     967         102 :         return true;
     968             :     else
     969          17 :         return GDALIsValueExactAs<T>(static_cast<double>(noDataValue));
     970             : }
     971             : 
     972             : template <bool IS_MAX, class NODATAType>
     973         226 : size_t extremum_element(const void *buffer, size_t nElts, GDALDataType eDT,
     974             :                         bool bHasNoData, NODATAType noDataValue)
     975             : {
     976         226 :     switch (eDT)
     977             :     {
     978             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 7, 0)
     979          15 :         case GDT_Int8:
     980             :         {
     981             :             using T = int8_t;
     982          15 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
     983          15 :             return extremum_element<T, IS_MAX>(
     984             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
     985          15 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
     986             :         }
     987             : #endif
     988          35 :         case GDT_Byte:
     989             :         {
     990             :             using T = uint8_t;
     991          35 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
     992          35 :             return extremum_element<T, IS_MAX>(
     993             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
     994          35 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
     995             :         }
     996          15 :         case GDT_Int16:
     997             :         {
     998             :             using T = int16_t;
     999          15 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1000          15 :             return extremum_element<T, IS_MAX>(
    1001             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1002          15 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1003             :         }
    1004          15 :         case GDT_UInt16:
    1005             :         {
    1006             :             using T = uint16_t;
    1007          15 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1008          15 :             return extremum_element<T, IS_MAX>(
    1009             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1010          15 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1011             :         }
    1012          15 :         case GDT_Int32:
    1013             :         {
    1014             :             using T = int32_t;
    1015          15 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1016          15 :             return extremum_element<T, IS_MAX>(
    1017             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1018          15 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1019             :         }
    1020          15 :         case GDT_UInt32:
    1021             :         {
    1022             :             using T = uint32_t;
    1023          15 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1024          15 :             return extremum_element<T, IS_MAX>(
    1025             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1026          15 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1027             :         }
    1028             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 5, 0)
    1029          15 :         case GDT_Int64:
    1030             :         {
    1031             :             using T = int64_t;
    1032          15 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1033          15 :             return extremum_element<T, IS_MAX>(
    1034             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1035          15 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1036             :         }
    1037          15 :         case GDT_UInt64:
    1038             :         {
    1039             :             using T = uint64_t;
    1040          15 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1041          15 :             return extremum_element<T, IS_MAX>(
    1042             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1043          15 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1044             :         }
    1045             : #endif
    1046             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
    1047          30 :         case GDT_Float16:
    1048             :         {
    1049             :             using T = GFloat16;
    1050          30 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1051          30 :             return extremum_element<T, IS_MAX>(
    1052             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1053          30 :                 bHasNoData ? static_cast<T>(noDataValue) : static_cast<T>(0));
    1054             :         }
    1055             : #endif
    1056          29 :         case GDT_Float32:
    1057             :         {
    1058             :             using T = float;
    1059          29 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1060          29 :             return extremum_element<T, IS_MAX>(
    1061             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1062          29 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1063             :         }
    1064          25 :         case GDT_Float64:
    1065             :         {
    1066             :             using T = double;
    1067          25 :             bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
    1068          25 :             return extremum_element<T, IS_MAX>(
    1069             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1070          25 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1071             :         }
    1072           2 :         case GDT_CInt16:
    1073             :         case GDT_CInt32:
    1074             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
    1075             :         case GDT_CFloat16:
    1076             : #endif
    1077             :         case GDT_CFloat32:
    1078             :         case GDT_CFloat64:
    1079             :         case GDT_Unknown:
    1080             :         case GDT_TypeCount:
    1081           2 :             break;
    1082             :     }
    1083           2 :     CPLError(CE_Failure, CPLE_NotSupported,
    1084             :              "%s not supported for this data type.", __FUNCTION__);
    1085           2 :     return 0;
    1086             : }
    1087             : 
    1088             : }  // namespace detail
    1089             : 
    1090             : /************************************************************************/
    1091             : /*                            max_element()                             */
    1092             : /************************************************************************/
    1093             : 
    1094             : /** Return the index of the element where the maximum value is hit.
    1095             :  *
    1096             :  * If it is hit in several locations, it is not specified which one will be
    1097             :  * returned.
    1098             :  *
    1099             :  * @param buffer Vector of nElts elements of type eDT.
    1100             :  * @param nElts Number of elements in buffer.
    1101             :  * @param eDT Data type of the elements of buffer.
    1102             :  * @param bHasNoData Whether noDataValue is valid.
    1103             :  * @param noDataValue Nodata value, only taken into account if bHasNoData == true
    1104             :  *
    1105             :  * @since GDAL 3.11
    1106             :  */
    1107             : template <class NODATAType>
    1108          97 : inline size_t max_element(const void *buffer, size_t nElts, GDALDataType eDT,
    1109             :                           bool bHasNoData, NODATAType noDataValue)
    1110             : {
    1111          97 :     return detail::extremum_element<true>(buffer, nElts, eDT, bHasNoData,
    1112          97 :                                           noDataValue);
    1113             : }
    1114             : 
    1115             : /************************************************************************/
    1116             : /*                            min_element()                             */
    1117             : /************************************************************************/
    1118             : 
    1119             : /** Return the index of the element where the minimum value is hit.
    1120             :  *
    1121             :  * If it is hit in several locations, it is not specified which one will be
    1122             :  * returned.
    1123             :  *
    1124             :  * @param buffer Vector of nElts elements of type eDT.
    1125             :  * @param nElts Number of elements in buffer.
    1126             :  * @param eDT Data type of the elements of buffer.
    1127             :  * @param bHasNoData Whether noDataValue is valid.
    1128             :  * @param noDataValue Nodata value, only taken into account if bHasNoData == true
    1129             :  *
    1130             :  * @since GDAL 3.11
    1131             :  */
    1132             : template <class NODATAType>
    1133         129 : inline size_t min_element(const void *buffer, size_t nElts, GDALDataType eDT,
    1134             :                           bool bHasNoData, NODATAType noDataValue)
    1135             : {
    1136         129 :     return detail::extremum_element<false>(buffer, nElts, eDT, bHasNoData,
    1137         129 :                                            noDataValue);
    1138             : }
    1139             : 
    1140             : namespace detail
    1141             : {
    1142             : 
    1143             : #ifdef NOT_EFFICIENT
    1144             : 
    1145             : /************************************************************************/
    1146             : /*                         minmax_element()                             */
    1147             : /************************************************************************/
    1148             : 
    1149             : template <class T>
    1150             : std::pair<size_t, size_t> minmax_element(const T *v, size_t size, T noDataValue)
    1151             : {
    1152             :     static_assert(!(std::is_floating_point_v<T>));
    1153             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
    1154             :     static_assert(!std::is_same_v<T, GFloat16>);
    1155             : #endif
    1156             :     if (size == 0)
    1157             :         return std::pair(0, 0);
    1158             :     size_t idx_of_min = 0;
    1159             :     size_t idx_of_max = 0;
    1160             :     T vmin = v[0];
    1161             :     T vmax = v[0];
    1162             :     bool extremum_is_nodata = vmin == noDataValue;
    1163             :     size_t i = 1;
    1164             :     for (; i < size; ++i)
    1165             :     {
    1166             :         if (v[i] != noDataValue && (v[i] < vmin || extremum_is_nodata))
    1167             :         {
    1168             :             vmin = v[i];
    1169             :             idx_of_min = i;
    1170             :             extremum_is_nodata = false;
    1171             :         }
    1172             :         if (v[i] != noDataValue && (v[i] > vmax || extremum_is_nodata))
    1173             :         {
    1174             :             vmax = v[i];
    1175             :             idx_of_max = i;
    1176             :             extremum_is_nodata = false;
    1177             :         }
    1178             :     }
    1179             :     return std::pair(idx_of_min, idx_of_max);
    1180             : }
    1181             : 
    1182             : template <class T>
    1183             : std::pair<size_t, size_t> minmax_element(const T *v, size_t size)
    1184             : {
    1185             :     static_assert(!(std::is_floating_point_v<T>));
    1186             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
    1187             :     static_assert(!std::is_same_v<T, GFloat16>);
    1188             : #endif
    1189             :     if (size == 0)
    1190             :         return std::pair(0, 0);
    1191             :     size_t idx_of_min = 0;
    1192             :     size_t idx_of_max = 0;
    1193             :     T vmin = v[0];
    1194             :     T vmax = v[0];
    1195             :     size_t i = 1;
    1196             :     for (; i < size; ++i)
    1197             :     {
    1198             :         if (v[i] < vmin)
    1199             :         {
    1200             :             vmin = v[i];
    1201             :             idx_of_min = i;
    1202             :         }
    1203             :         if (v[i] > vmax)
    1204             :         {
    1205             :             vmax = v[i];
    1206             :             idx_of_max = i;
    1207             :         }
    1208             :     }
    1209             :     return std::pair(idx_of_min, idx_of_max);
    1210             : }
    1211             : 
    1212             : template <class T>
    1213             : inline std::pair<size_t, size_t> minmax_element_with_nan(const T *v,
    1214             :                                                          size_t size)
    1215             : {
    1216             :     if (size == 0)
    1217             :         return std::pair(0, 0);
    1218             :     size_t idx_of_min = 0;
    1219             :     size_t idx_of_max = 0;
    1220             :     T vmin = v[0];
    1221             :     T vmax = v[0];
    1222             :     size_t i = 1;
    1223             :     if (IsNan(v[0]))
    1224             :     {
    1225             :         for (; i < size; ++i)
    1226             :         {
    1227             :             if (!IsNan(v[i]))
    1228             :             {
    1229             :                 vmin = v[i];
    1230             :                 idx_of_min = i;
    1231             :                 vmax = v[i];
    1232             :                 idx_of_max = i;
    1233             :                 break;
    1234             :             }
    1235             :         }
    1236             :     }
    1237             :     for (; i < size; ++i)
    1238             :     {
    1239             :         if (v[i] < vmin)
    1240             :         {
    1241             :             vmin = v[i];
    1242             :             idx_of_min = i;
    1243             :         }
    1244             :         if (v[i] > vmax)
    1245             :         {
    1246             :             vmax = v[i];
    1247             :             idx_of_max = i;
    1248             :         }
    1249             :     }
    1250             :     return std::pair(idx_of_min, idx_of_max);
    1251             : }
    1252             : 
    1253             : template <>
    1254             : std::pair<size_t, size_t> minmax_element<float>(const float *v, size_t size)
    1255             : {
    1256             :     return minmax_element_with_nan<float>(v, size);
    1257             : }
    1258             : 
    1259             : template <>
    1260             : std::pair<size_t, size_t> minmax_element<double>(const double *v, size_t size)
    1261             : {
    1262             :     return minmax_element_with_nan<double>(v, size);
    1263             : }
    1264             : 
    1265             : template <class T>
    1266             : inline std::pair<size_t, size_t> minmax_element(const T *buffer, size_t size,
    1267             :                                                 bool bHasNoData, T noDataValue)
    1268             : {
    1269             :     if (bHasNoData)
    1270             :     {
    1271             :         return minmax_element<T>(buffer, size, noDataValue);
    1272             :     }
    1273             :     else
    1274             :     {
    1275             :         return minmax_element<T>(buffer, size);
    1276             :     }
    1277             : }
    1278             : #else
    1279             : 
    1280             : /************************************************************************/
    1281             : /*                         minmax_element()                             */
    1282             : /************************************************************************/
    1283             : 
    1284             : template <class T>
    1285          16 : inline std::pair<size_t, size_t> minmax_element(const T *buffer, size_t size,
    1286             :                                                 bool bHasNoData, T noDataValue)
    1287             : {
    1288             : #ifdef NOT_EFFICIENT
    1289             :     if (bHasNoData)
    1290             :     {
    1291             :         return minmax_element<T>(buffer, size, noDataValue);
    1292             :     }
    1293             :     else
    1294             :     {
    1295             :         return minmax_element<T>(buffer, size);
    1296             :         //auto [imin, imax] = std::minmax_element(buffer, buffer + size);
    1297             :         //return std::pair(imin - buffer, imax - buffer);
    1298             :     }
    1299             : #else
    1300             : 
    1301             : #if !defined(GDAL_MINMAX_ELEMENT_USE_SSE2)
    1302             :     if constexpr (!std::is_floating_point_v<T>
    1303             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
    1304             :                   && !std::is_same_v<T, GFloat16>
    1305             : #endif
    1306             :     )
    1307             :     {
    1308             :         if (!bHasNoData)
    1309             :         {
    1310             :             auto [min_iter, max_iter] =
    1311             :                 std::minmax_element(buffer, buffer + size);
    1312             :             return std::pair(min_iter - buffer, max_iter - buffer);
    1313             :         }
    1314             :     }
    1315             : #endif
    1316             : 
    1317             :     // Using separately min and max is more efficient than computing them
    1318             :     // within the same loop
    1319             :     return std::pair(
    1320          48 :         extremum_element<T, false>(buffer, size, bHasNoData, noDataValue),
    1321          16 :         extremum_element<T, true>(buffer, size, bHasNoData, noDataValue));
    1322             : 
    1323             : #endif
    1324             : }
    1325             : #endif
    1326             : 
    1327             : }  // namespace detail
    1328             : 
    1329             : /************************************************************************/
    1330             : /*                          minmax_element()                            */
    1331             : /************************************************************************/
    1332             : 
    1333             : /** Return the index of the elements where the minimum and maximum values are hit.
    1334             :  *
    1335             :  * If they are hit in several locations, it is not specified which one will be
    1336             :  * returned (contrary to std::minmax_element).
    1337             :  *
    1338             :  * @param buffer Vector of nElts elements of type eDT.
    1339             :  * @param nElts Number of elements in buffer.
    1340             :  * @param eDT Data type of the elements of buffer.
    1341             :  * @param bHasNoData Whether noDataValue is valid.
    1342             :  * @param noDataValue Nodata value, only taken into account if bHasNoData == true
    1343             :  *
    1344             :  * @since GDAL 3.11
    1345             :  */
    1346             : template <class NODATAType>
    1347             : inline std::pair<size_t, size_t>
    1348          17 : minmax_element(const void *buffer, size_t nElts, GDALDataType eDT,
    1349             :                bool bHasNoData, NODATAType noDataValue)
    1350             : {
    1351          17 :     switch (eDT)
    1352             :     {
    1353             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 7, 0)
    1354           1 :         case GDT_Int8:
    1355             :         {
    1356             :             using T = int8_t;
    1357           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1358           1 :             return detail::minmax_element<T>(
    1359             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1360           1 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1361             :         }
    1362             : #endif
    1363           5 :         case GDT_Byte:
    1364             :         {
    1365             :             using T = uint8_t;
    1366           5 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1367           5 :             return detail::minmax_element<T>(
    1368             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1369           5 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1370             :         }
    1371           1 :         case GDT_Int16:
    1372             :         {
    1373             :             using T = int16_t;
    1374           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1375           1 :             return detail::minmax_element<T>(
    1376             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1377           1 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1378             :         }
    1379           1 :         case GDT_UInt16:
    1380             :         {
    1381             :             using T = uint16_t;
    1382           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1383           1 :             return detail::minmax_element<T>(
    1384             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1385           1 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1386             :         }
    1387           1 :         case GDT_Int32:
    1388             :         {
    1389             :             using T = int32_t;
    1390           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1391           1 :             return detail::minmax_element<T>(
    1392             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1393           1 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1394             :         }
    1395           1 :         case GDT_UInt32:
    1396             :         {
    1397             :             using T = uint32_t;
    1398           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1399           1 :             return detail::minmax_element<T>(
    1400             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1401           1 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1402             :         }
    1403             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 5, 0)
    1404           1 :         case GDT_Int64:
    1405             :         {
    1406             :             using T = int64_t;
    1407           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1408           1 :             return detail::minmax_element<T>(
    1409             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1410           1 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1411             :         }
    1412           1 :         case GDT_UInt64:
    1413             :         {
    1414             :             using T = uint64_t;
    1415           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1416           1 :             return detail::minmax_element<T>(
    1417             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1418           1 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1419             :         }
    1420             : #endif
    1421             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
    1422           1 :         case GDT_Float16:
    1423             :         {
    1424             :             using T = GFloat16;
    1425           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1426           1 :             return detail::minmax_element<T>(
    1427             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1428           1 :                 bHasNoData ? static_cast<T>(noDataValue) : static_cast<T>(0));
    1429             :         }
    1430             : #endif
    1431           1 :         case GDT_Float32:
    1432             :         {
    1433             :             using T = float;
    1434           1 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1435           1 :             return detail::minmax_element<T>(
    1436             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1437           1 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1438             :         }
    1439           2 :         case GDT_Float64:
    1440             :         {
    1441             :             using T = double;
    1442           2 :             bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
    1443           2 :             return detail::minmax_element<T>(
    1444             :                 static_cast<const T *>(buffer), nElts, bHasNoData,
    1445           2 :                 bHasNoData ? static_cast<T>(noDataValue) : 0);
    1446             :         }
    1447           1 :         case GDT_CInt16:
    1448             :         case GDT_CInt32:
    1449             : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
    1450             :         case GDT_CFloat16:
    1451             : #endif
    1452             :         case GDT_CFloat32:
    1453             :         case GDT_CFloat64:
    1454             :         case GDT_Unknown:
    1455             :         case GDT_TypeCount:
    1456           1 :             break;
    1457             :     }
    1458           1 :     CPLError(CE_Failure, CPLE_NotSupported,
    1459             :              "%s not supported for this data type.", __FUNCTION__);
    1460           1 :     return std::pair(0, 0);
    1461             : }
    1462             : 
    1463             : }  // namespace GDAL_MINMAXELT_NS
    1464             : 
    1465             : #endif  // GDAL_MINMAX_ELEMENT_INCLUDED

Generated by: LCOV version 1.14