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 484711 : template <class T, bool IS_MAX> inline static bool compScalar(T x, T y)
86 : {
87 : if constexpr (IS_MAX)
88 193463 : return x > y;
89 : else
90 291248 : return x < y;
91 : }
92 :
93 14834 : 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 14834 : return isnan(x);
102 : }
103 :
104 6882449 : template <class T> inline static bool compEqual(T x, T y)
105 : {
106 6882449 : 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 10662 : template <> bool IsNan<GFloat16>(GFloat16 x)
114 : {
115 : uint16_t iX;
116 10662 : memcpy(&iX, &x, sizeof(x));
117 : // Check that the 5 bits of the exponent are set and the mantissa is not zero.
118 10662 : return (iX & 0x7fff) > 0x7c00;
119 : }
120 :
121 1826 : template <> bool compEqual<GFloat16>(GFloat16 x, GFloat16 y)
122 : {
123 : uint16_t iX, iY;
124 1826 : memcpy(&iX, &x, sizeof(x));
125 1826 : memcpy(&iY, &y, sizeof(y));
126 : // Given our usage where y cannot be NaN we can skip the IsNan tests
127 1826 : assert(!IsNan(y));
128 3644 : return (iX == iY ||
129 : // Also check +0 == -0
130 3644 : ((iX | iY) == (1 << 15)));
131 : }
132 :
133 4410 : template <> bool compScalar<GFloat16, true>(GFloat16 x, GFloat16 y)
134 : {
135 : uint16_t iX, iY;
136 4410 : memcpy(&iX, &x, sizeof(x));
137 4410 : memcpy(&iY, &y, sizeof(y));
138 : bool ret;
139 4410 : if (IsNan(x) || IsNan(y))
140 : {
141 831 : ret = false;
142 : }
143 3579 : else if (!(iX >> 15))
144 : {
145 : // +0 considered > -0. We don't really care
146 167 : ret = (iY >> 15) || iX > iY;
147 : }
148 : else
149 : {
150 3412 : ret = (iY >> 15) && iX < iY;
151 : }
152 4410 : return ret;
153 : }
154 :
155 1869 : template <> bool compScalar<GFloat16, false>(GFloat16 x, GFloat16 y)
156 : {
157 1869 : 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 90 : static inline int8_t Shift8(uint8_t x)
357 : {
358 90 : return static_cast<int8_t>(x + std::numeric_limits<int8_t>::min());
359 : }
360 :
361 46 : static inline int16_t Shift16(uint16_t x)
362 : {
363 46 : return static_cast<int16_t>(x + std::numeric_limits<int16_t>::min());
364 : }
365 :
366 : CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW
367 663 : static inline int32_t Shift32(uint32_t x)
368 : {
369 663 : x += static_cast<uint32_t>(std::numeric_limits<int32_t>::min());
370 : int32_t ret;
371 663 : memcpy(&ret, &x, sizeof(x));
372 663 : return ret;
373 : }
374 :
375 : CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW
376 136 : static inline int64_t Shift64(uint64_t x)
377 : {
378 136 : x += static_cast<uint64_t>(std::numeric_limits<int64_t>::min());
379 : int64_t ret;
380 136 : memcpy(&ret, &x, sizeof(x));
381 136 : return ret;
382 : }
383 :
384 : // Return a _mm128[i|d] register with all its elements set to x
385 33669 : template <class T> static inline auto set1(T x)
386 : {
387 : if constexpr (std::is_same_v<T, uint8_t>)
388 180 : return _mm_set1_epi8(Shift8(x));
389 : else if constexpr (std::is_same_v<T, int8_t>)
390 470 : return _mm_set1_epi8(x);
391 : else if constexpr (std::is_same_v<T, uint16_t>)
392 92 : return _mm_set1_epi16(Shift16(x));
393 : else if constexpr (std::is_same_v<T, int16_t>)
394 9178 : return _mm_set1_epi16(x);
395 : else if constexpr (std::is_same_v<T, uint32_t>)
396 1326 : return _mm_set1_epi32(Shift32(x));
397 : else if constexpr (std::is_same_v<T, int32_t>)
398 3170 : return _mm_set1_epi32(x);
399 : else if constexpr (std::is_same_v<T, uint64_t>)
400 272 : return _mm_set1_epi64x(Shift64(x));
401 : else if constexpr (std::is_same_v<T, int64_t>)
402 187 : 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 139 : memcpy(&iX, &x, sizeof(x));
408 278 : return _mm_set1_epi16(iX);
409 : }
410 : #endif
411 : else if constexpr (std::is_same_v<T, float>)
412 2713 : return _mm_set1_ps(x);
413 : else
414 : {
415 : static_assert(std::is_same_v<T, double>);
416 21701 : return _mm_set1_pd(x);
417 : }
418 : }
419 :
420 : // Return a _mm128[i|d] register with all its elements set to x
421 38329 : template <class T> static inline auto set1_unshifted(T x)
422 : {
423 : if constexpr (std::is_same_v<T, uint8_t>)
424 : {
425 : int8_t xSigned;
426 111 : memcpy(&xSigned, &x, sizeof(xSigned));
427 222 : return _mm_set1_epi8(xSigned);
428 : }
429 : else if constexpr (std::is_same_v<T, int8_t>)
430 928 : return _mm_set1_epi8(x);
431 : else if constexpr (std::is_same_v<T, uint16_t>)
432 : {
433 : int16_t xSigned;
434 54 : memcpy(&xSigned, &x, sizeof(xSigned));
435 108 : return _mm_set1_epi16(xSigned);
436 : }
437 : else if constexpr (std::is_same_v<T, int16_t>)
438 19278 : return _mm_set1_epi16(x);
439 : else if constexpr (std::is_same_v<T, uint32_t>)
440 : {
441 : int32_t xSigned;
442 928 : memcpy(&xSigned, &x, sizeof(xSigned));
443 1856 : return _mm_set1_epi32(xSigned);
444 : }
445 : else if constexpr (std::is_same_v<T, int32_t>)
446 5284 : return _mm_set1_epi32(x);
447 : else if constexpr (std::is_same_v<T, uint64_t>)
448 : {
449 : int64_t xSigned;
450 124 : memcpy(&xSigned, &x, sizeof(xSigned));
451 248 : return _mm_set1_epi64x(xSigned);
452 : }
453 : else if constexpr (std::is_same_v<T, int64_t>)
454 146 : return _mm_set1_epi64x(x);
455 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
456 : else if constexpr (std::is_same_v<T, GFloat16>)
457 : {
458 : int16_t iX;
459 162 : memcpy(&iX, &x, sizeof(x));
460 324 : return _mm_set1_epi16(iX);
461 : }
462 : #endif
463 : else if constexpr (std::is_same_v<T, float>)
464 3925 : return _mm_set1_ps(x);
465 : else
466 : {
467 : static_assert(std::is_same_v<T, double>);
468 17492 : return _mm_set1_pd(x);
469 : }
470 : }
471 :
472 : // Load as many values of type T at a _mm128[i|d] register can contain from x
473 145541938 : template <class T> static inline auto loadv(const T *x)
474 : {
475 : if constexpr (std::is_same_v<T, float>)
476 20045636 : return _mm_loadu_ps(x);
477 : else if constexpr (std::is_same_v<T, double>)
478 40287388 : return _mm_loadu_pd(x);
479 : else
480 85208914 : return _mm_loadu_si128(reinterpret_cast<const __m128i *>(x));
481 : }
482 :
483 15000144 : inline __m128i IsNanGFloat16(__m128i x)
484 : {
485 : // (iX & 0x7fff) > 0x7c00
486 45000332 : return _mm_cmpgt_epi16(_mm_and_si128(x, _mm_set1_epi16(0x7fff)),
487 15000144 : _mm_set1_epi16(0x7c00));
488 : }
489 :
490 : template <class T, class SSE_T>
491 82662930 : static inline SSE_T blendv(SSE_T a, SSE_T b, SSE_T mask)
492 : {
493 : if constexpr (std::is_same_v<T, float>)
494 : {
495 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
496 : return _mm_blendv_ps(a, b, mask);
497 : #else
498 30000180 : return _mm_or_ps(_mm_andnot_ps(mask, a), _mm_and_ps(mask, b));
499 : #endif
500 : }
501 : else if constexpr (std::is_same_v<T, double>)
502 : {
503 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
504 : return _mm_blendv_pd(a, b, mask);
505 : #else
506 60064020 : return _mm_or_pd(_mm_andnot_pd(mask, a), _mm_and_pd(mask, b));
507 : #endif
508 : }
509 : else
510 : {
511 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
512 : return _mm_blendv_epi8(a, b, mask);
513 : #else
514 157924760 : return _mm_or_si128(_mm_andnot_si128(mask, a), _mm_and_si128(mask, b));
515 : #endif
516 : }
517 : }
518 :
519 10000058 : inline __m128i cmpgt_ph(__m128i x, __m128i y)
520 : {
521 : #ifdef slow
522 : GFloat16 vx[8], vy[8];
523 : int16_t res[8];
524 : _mm_storeu_si128(reinterpret_cast<__m128i *>(vx), x);
525 : _mm_storeu_si128(reinterpret_cast<__m128i *>(vy), y);
526 : for (int i = 0; i < 8; ++i)
527 : {
528 : res[i] = vx[i] > vy[i] ? -1 : 0;
529 : }
530 : return _mm_loadu_si128(reinterpret_cast<const __m128i *>(res));
531 : #else
532 10000058 : const auto x_is_negative = _mm_srai_epi16(x, 15);
533 10000058 : const auto y_is_negative = _mm_srai_epi16(y, 15);
534 40000252 : return _mm_andnot_si128(
535 : // only x can be NaN given how we use this method
536 : IsNanGFloat16(x),
537 : blendv<int16_t>(_mm_or_si128(y_is_negative, _mm_cmpgt_epi16(x, y)),
538 : _mm_and_si128(y_is_negative, _mm_cmpgt_epi16(y, x)),
539 10000058 : x_is_negative));
540 : #endif
541 : }
542 :
543 40007320 : inline __m128i cmpgt_epi64(__m128i x, __m128i y)
544 : {
545 : #if defined(__SSE4_2__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
546 : return _mm_cmpgt_epi64(x, y);
547 : #else
548 120021960 : auto tmp = _mm_and_si128(_mm_sub_epi64(y, x), _mm_cmpeq_epi32(x, y));
549 40007320 : tmp = _mm_or_si128(tmp, _mm_cmpgt_epi32(x, y));
550 : // Replicate the 2 odd-indexed (hi) 32-bit word into the 2 even-indexed (lo)
551 : // ones
552 40007320 : return _mm_shuffle_epi32(tmp, _MM_SHUFFLE(3, 3, 1, 1));
553 : #endif
554 : }
555 :
556 : // Return a __m128i register with bits set when x[i] < y[i] when !IS_MAX
557 : // or x[i] > y[i] when IS_MAX
558 : template <class T, bool IS_MAX, class SSE_T>
559 145541848 : static inline __m128i comp(SSE_T x, SSE_T y)
560 : {
561 : if constexpr (IS_MAX)
562 : {
563 : if constexpr (std::is_same_v<T, uint8_t>)
564 3750472 : return _mm_cmpgt_epi8(
565 : _mm_add_epi8(x,
566 1250154 : _mm_set1_epi8(std::numeric_limits<int8_t>::min())),
567 1250154 : y);
568 : else if constexpr (std::is_same_v<T, int8_t>)
569 1250422 : return _mm_cmpgt_epi8(x, y);
570 : else if constexpr (std::is_same_v<T, uint16_t>)
571 7500064 : return _mm_cmpgt_epi16(
572 : _mm_add_epi16(
573 2500018 : x, _mm_set1_epi16(std::numeric_limits<int16_t>::min())),
574 2500018 : y);
575 : else if constexpr (std::is_same_v<T, int16_t>)
576 2523946 : return _mm_cmpgt_epi16(x, y);
577 : else if constexpr (std::is_same_v<T, uint32_t>)
578 15010920 : return _mm_cmpgt_epi32(
579 : _mm_add_epi32(
580 : x, _mm_set1_epi32(std::numeric_limits<int32_t>::min())),
581 5003630 : y);
582 : else if constexpr (std::is_same_v<T, int32_t>)
583 5072498 : return _mm_cmpgt_epi32(x, y);
584 : else if constexpr (std::is_same_v<T, uint64_t>)
585 : {
586 20002736 : return cmpgt_epi64(
587 : _mm_add_epi64(
588 10001358 : x, _mm_set1_epi64x(std::numeric_limits<int64_t>::min())),
589 10001358 : y);
590 : }
591 : else if constexpr (std::is_same_v<T, int64_t>)
592 : {
593 10002158 : return cmpgt_epi64(x, y);
594 : }
595 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
596 : else if constexpr (std::is_same_v<T, GFloat16>)
597 : {
598 5000024 : return cmpgt_ph(x, y);
599 : }
600 : #endif
601 : else if constexpr (std::is_same_v<T, float>)
602 20045576 : return _mm_castps_si128(_mm_cmpgt_ps(x, y));
603 : else
604 : {
605 : static_assert(std::is_same_v<T, double>);
606 40292264 : return _mm_castpd_si128(_mm_cmpgt_pd(x, y));
607 : }
608 : }
609 : else
610 : {
611 : if constexpr (std::is_same_v<T, uint8_t>)
612 3750508 : return _mm_cmplt_epi8(
613 : _mm_add_epi8(x,
614 1250166 : _mm_set1_epi8(std::numeric_limits<int8_t>::min())),
615 1250166 : y);
616 : else if constexpr (std::is_same_v<T, int8_t>)
617 1250434 : return _mm_cmplt_epi8(x, y);
618 : else if constexpr (std::is_same_v<T, uint16_t>)
619 7500148 : return _mm_cmplt_epi16(
620 : _mm_add_epi16(
621 2500046 : x, _mm_set1_epi16(std::numeric_limits<int16_t>::min())),
622 2500046 : y);
623 : else if constexpr (std::is_same_v<T, int16_t>)
624 2523974 : return _mm_cmplt_epi16(x, y);
625 : else if constexpr (std::is_same_v<T, uint32_t>)
626 15011100 : return _mm_cmplt_epi32(
627 : _mm_add_epi32(
628 : x, _mm_set1_epi32(std::numeric_limits<int32_t>::min())),
629 5003690 : y);
630 : else if constexpr (std::is_same_v<T, int32_t>)
631 5072558 : return _mm_cmplt_epi32(x, y);
632 : else if constexpr (std::is_same_v<T, uint64_t>)
633 : {
634 20002984 : return cmpgt_epi64(
635 : y, _mm_add_epi64(x, _mm_set1_epi64x(
636 20002984 : std::numeric_limits<int64_t>::min())));
637 : }
638 : else if constexpr (std::is_same_v<T, int64_t>)
639 : {
640 10002282 : return cmpgt_epi64(y, x);
641 : }
642 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
643 : else if constexpr (std::is_same_v<T, GFloat16>)
644 : {
645 5000024 : return cmpgt_ph(y, x);
646 : }
647 : #endif
648 : else if constexpr (std::is_same_v<T, float>)
649 20045696 : return _mm_castps_si128(_mm_cmplt_ps(x, y));
650 : else
651 : {
652 : static_assert(std::is_same_v<T, double>);
653 40282512 : return _mm_castpd_si128(_mm_cmplt_pd(x, y));
654 : }
655 : }
656 : }
657 :
658 : template <class T, class SSE_T> static inline SSE_T compeq(SSE_T a, SSE_T b);
659 :
660 1250018 : template <> __m128i compeq<uint8_t>(__m128i a, __m128i b)
661 : {
662 1250018 : return _mm_cmpeq_epi8(a, b);
663 : }
664 :
665 1250002 : template <> __m128i compeq<int8_t>(__m128i a, __m128i b)
666 : {
667 1250002 : return _mm_cmpeq_epi8(a, b);
668 : }
669 :
670 2500018 : template <> __m128i compeq<uint16_t>(__m128i a, __m128i b)
671 : {
672 2500018 : return _mm_cmpeq_epi16(a, b);
673 : }
674 :
675 2544114 : template <> __m128i compeq<int16_t>(__m128i a, __m128i b)
676 : {
677 2544114 : return _mm_cmpeq_epi16(a, b);
678 : }
679 :
680 5000058 : template <> __m128i compeq<uint32_t>(__m128i a, __m128i b)
681 : {
682 5000058 : return _mm_cmpeq_epi32(a, b);
683 : }
684 :
685 5096898 : template <> __m128i compeq<int32_t>(__m128i a, __m128i b)
686 : {
687 5096898 : return _mm_cmpeq_epi32(a, b);
688 : }
689 :
690 20000328 : template <> __m128i compeq<uint64_t>(__m128i a, __m128i b)
691 : {
692 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
693 : return _mm_cmpeq_epi64(a, b);
694 : #else
695 20000328 : auto tmp = _mm_cmpeq_epi32(a, b);
696 : // The shuffle swaps hi-lo 32-bit words
697 40000656 : return _mm_and_si128(tmp, _mm_shuffle_epi32(tmp, _MM_SHUFFLE(2, 3, 0, 1)));
698 : #endif
699 : }
700 :
701 : template <>
702 : #if defined(__INTEL_CLANG_COMPILER) && \
703 : !(defined(__SSE4_1__) || defined(__AVX__))
704 : // ICC 2024 has a bug with the following code when -fiopenmp is enabled.
705 : // Disabling inlining works around it...
706 : __attribute__((noinline))
707 : #endif
708 : __m128i
709 10000170 : compeq<int64_t>(__m128i a, __m128i b)
710 : {
711 10000170 : return compeq<uint64_t>(a, b);
712 : }
713 :
714 10000040 : template <> __m128 compeq<float>(__m128 a, __m128 b)
715 : {
716 10000040 : return _mm_cmpeq_ps(a, b);
717 : }
718 :
719 20021340 : template <> __m128d compeq<double>(__m128d a, __m128d b)
720 : {
721 20021340 : return _mm_cmpeq_pd(a, b);
722 : }
723 :
724 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
725 5000036 : template <> __m128i compeq<GFloat16>(__m128i a, __m128i b)
726 : {
727 : // !isnan(a) && !isnan(b) && (a == b || (a|b) == 0x8000)
728 30000236 : return _mm_andnot_si128(
729 : IsNanGFloat16(a), // b cannot be NaN given how we use this method
730 : _mm_or_si128(_mm_cmpeq_epi16(a, b),
731 : _mm_cmpeq_epi16(
732 : _mm_or_si128(a, b),
733 10000082 : _mm_set1_epi16(std::numeric_limits<int16_t>::min()))));
734 : }
735 : #endif
736 :
737 : // Using SSE2
738 : template <class T, bool IS_MAX, bool HAS_NODATA>
739 : #if defined(__GNUC__)
740 : __attribute__((noinline))
741 : #endif
742 : size_t
743 13966 : extremum_element_sse2(const T *v, size_t size, T noDataValue)
744 : {
745 : static_assert(std::is_same_v<T, uint8_t> || std::is_same_v<T, int8_t> ||
746 : std::is_same_v<T, uint16_t> || std::is_same_v<T, int16_t> ||
747 : std::is_same_v<T, uint32_t> || std::is_same_v<T, int32_t> ||
748 : std::is_same_v<T, uint64_t> || std::is_same_v<T, int64_t> ||
749 : is_floating_point_v<T>);
750 13966 : if (size == 0)
751 22 : return 0;
752 13944 : size_t idx_of_extremum = 0;
753 13944 : T extremum = v[0];
754 13944 : [[maybe_unused]] bool extremum_is_invalid = false;
755 : if constexpr (is_floating_point_v<T>)
756 : {
757 : if constexpr (HAS_NODATA)
758 : {
759 5346 : if (IsNan(noDataValue))
760 : {
761 6 : return extremum_element_sse2<T, IS_MAX, false>(
762 4 : v, size, static_cast<T>(0));
763 : }
764 : }
765 8228 : extremum_is_invalid = IsNan(extremum);
766 : }
767 : if constexpr (HAS_NODATA)
768 : {
769 8793 : if (compEqual(extremum, noDataValue))
770 4235 : extremum_is_invalid = true;
771 : }
772 13940 : size_t i = 1;
773 :
774 13940 : constexpr size_t VALS_PER_REG = sizeof(set1(extremum)) / sizeof(extremum);
775 13940 : constexpr int LOOP_UNROLLING = 4;
776 : // If changing the value, then we need to adjust the number of sse_valX
777 : // loading in the loop.
778 : static_assert(LOOP_UNROLLING == 4);
779 13940 : constexpr size_t VALS_PER_ITER = VALS_PER_REG * LOOP_UNROLLING;
780 :
781 8389405 : const auto update = [v, noDataValue, &extremum, &idx_of_extremum,
782 : &extremum_is_invalid](size_t idx)
783 : {
784 : if constexpr (HAS_NODATA)
785 : {
786 6875482 : if (compEqual(v[idx], noDataValue))
787 6649363 : return;
788 226121 : if (extremum_is_invalid)
789 : {
790 : if constexpr (is_floating_point_v<T>)
791 : {
792 32 : if (IsNan(v[idx]))
793 6 : return;
794 : }
795 50 : extremum = v[idx];
796 50 : idx_of_extremum = idx;
797 50 : extremum_is_invalid = false;
798 50 : return;
799 : }
800 : }
801 : else
802 : {
803 263056 : CPL_IGNORE_RET_VAL(noDataValue);
804 : }
805 489121 : if (compScalar<T, IS_MAX>(v[idx], extremum))
806 : {
807 155488 : extremum = v[idx];
808 155488 : idx_of_extremum = idx;
809 155488 : extremum_is_invalid = false;
810 : }
811 : else if constexpr (is_floating_point_v<T>)
812 : {
813 100973 : if (extremum_is_invalid && !IsNan(v[idx]))
814 : {
815 28 : extremum = v[idx];
816 28 : idx_of_extremum = idx;
817 28 : extremum_is_invalid = false;
818 : }
819 : }
820 : };
821 :
822 200419 : for (; i < VALS_PER_ITER && i < size; ++i)
823 : {
824 186479 : update(i);
825 : }
826 :
827 13940 : [[maybe_unused]] auto sse_neutral = set1_unshifted(static_cast<T>(0));
828 13940 : [[maybe_unused]] auto sse_nodata = set1_unshifted(noDataValue);
829 : if constexpr (HAS_NODATA || is_floating_point_v<T>)
830 : {
831 6633076 : for (; i < size && extremum_is_invalid; ++i)
832 : {
833 6621399 : update(i);
834 : }
835 11679 : if (!extremum_is_invalid)
836 : {
837 7624 : for (; i < size && (i % VALS_PER_ITER) != 0; ++i)
838 : {
839 148 : update(i);
840 : }
841 7476 : sse_neutral = set1_unshifted(extremum);
842 : }
843 : }
844 :
845 13940 : auto sse_extremum = set1(extremum);
846 :
847 13940 : [[maybe_unused]] size_t hits = 0;
848 13940 : const auto sse_iter_count = (size / VALS_PER_ITER) * VALS_PER_ITER;
849 36398833 : for (; i < sse_iter_count; i += VALS_PER_ITER)
850 : {
851 : // A bit of loop unrolling to save 3/4 of slow movemask operations.
852 36385490 : auto sse_val0 = loadv(v + i + 0 * VALS_PER_REG);
853 36385490 : auto sse_val1 = loadv(v + i + 1 * VALS_PER_REG);
854 36385490 : auto sse_val2 = loadv(v + i + 2 * VALS_PER_REG);
855 36385490 : auto sse_val3 = loadv(v + i + 3 * VALS_PER_REG);
856 :
857 : if constexpr (HAS_NODATA)
858 : {
859 : // Replace all components that are at the nodata value by a
860 : // neutral value (current minimum)
861 18165722 : const auto replaceNoDataByNeutral =
862 145325664 : [sse_neutral, sse_nodata](auto sse_val)
863 : {
864 72662852 : const auto eq_nodata = compeq<T>(sse_val, sse_nodata);
865 72662852 : return blendv<T>(sse_val, sse_neutral, eq_nodata);
866 : };
867 :
868 18165722 : sse_val0 = replaceNoDataByNeutral(sse_val0);
869 18165722 : sse_val1 = replaceNoDataByNeutral(sse_val1);
870 18165722 : sse_val2 = replaceNoDataByNeutral(sse_val2);
871 18165722 : sse_val3 = replaceNoDataByNeutral(sse_val3);
872 : }
873 :
874 145541936 : if (_mm_movemask_epi8(_mm_or_si128(
875 : _mm_or_si128(comp<T, IS_MAX>(sse_val0, sse_extremum),
876 : comp<T, IS_MAX>(sse_val1, sse_extremum)),
877 : _mm_or_si128(comp<T, IS_MAX>(sse_val2, sse_extremum),
878 36385490 : comp<T, IS_MAX>(sse_val3, sse_extremum)))) != 0)
879 : {
880 : if constexpr (!std::is_same_v<T, int8_t> &&
881 : !std::is_same_v<T, uint8_t>)
882 : {
883 : // The above tests excluding int8_t/uint8_t is due to the fact
884 : // with those small ranges of values we will quickly converge
885 : // to the minimum, so no need to do the below "smart" test.
886 :
887 20245 : if (++hits == size / 16)
888 : {
889 : // If we have an almost sorted array, then using this code path
890 : // will hurt performance. Arbitrary give up if we get here
891 : // more than 1. / 16 of the size of the array.
892 : // fprintf(stderr, "going to non-vector path\n");
893 573 : break;
894 : }
895 : }
896 248473 : for (size_t j = 0; j < VALS_PER_ITER; j++)
897 : {
898 228744 : update(i + j);
899 : }
900 :
901 19729 : sse_extremum = set1(extremum);
902 : if constexpr (HAS_NODATA)
903 : {
904 2973 : sse_neutral = set1_unshifted(extremum);
905 : }
906 : }
907 : }
908 115712 : for (; i < size; ++i)
909 : {
910 101772 : update(i);
911 : }
912 13940 : return idx_of_extremum;
913 : }
914 :
915 : /************************************************************************/
916 : /* extremum_element() */
917 : /************************************************************************/
918 :
919 : template <class T, bool IS_MAX>
920 8808 : inline size_t extremum_element(const T *buffer, size_t size, T noDataValue)
921 : {
922 8808 : return extremum_element_sse2<T, IS_MAX, true>(buffer, size, noDataValue);
923 : }
924 :
925 : template <class T, bool IS_MAX>
926 5154 : inline size_t extremum_element(const T *buffer, size_t size)
927 : {
928 5177 : return extremum_element_sse2<T, IS_MAX, false>(buffer, size,
929 5177 : static_cast<T>(0));
930 : }
931 :
932 : #else
933 :
934 : template <class T, bool IS_MAX>
935 : inline size_t extremum_element(const T *buffer, size_t size, T noDataValue)
936 : {
937 : if constexpr (is_floating_point_v<T>)
938 : return extremum_element_with_nan_generic<T, IS_MAX>(buffer, size,
939 : noDataValue);
940 : else
941 : return extremum_element_generic<T, IS_MAX>(buffer, size, true,
942 : noDataValue);
943 : }
944 :
945 : template <class T, bool IS_MAX>
946 : inline size_t extremum_element(const T *buffer, size_t size)
947 : {
948 : if constexpr (is_floating_point_v<T>)
949 : return extremum_element_with_nan_generic<T, IS_MAX>(buffer, size);
950 : else
951 : return extremum_element_generic<T, IS_MAX>(buffer, size, false,
952 : static_cast<T>(0));
953 : }
954 :
955 : #endif
956 :
957 : /************************************************************************/
958 : /* extremum_element() */
959 : /************************************************************************/
960 :
961 : template <class T, bool IS_MAX>
962 13962 : inline size_t extremum_element(const T *buffer, size_t size, bool bHasNoData,
963 : T noDataValue)
964 : {
965 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0) && \
966 : !defined(GDAL_MINMAX_ELEMENT_USE_SSE2)
967 : if constexpr (std::is_same_v<T, GFloat16>)
968 : {
969 : if (bHasNoData)
970 : return extremum_element_with_nan_generic<T, IS_MAX>(buffer, size,
971 : noDataValue);
972 : else
973 : return extremum_element_with_nan_generic<T, IS_MAX>(buffer, size);
974 : }
975 : else
976 : #endif
977 13962 : if (bHasNoData)
978 8808 : return extremum_element<T, IS_MAX>(buffer, size, noDataValue);
979 : else
980 5154 : return extremum_element<T, IS_MAX>(buffer, size);
981 : }
982 :
983 : template <class T, class NODATAType>
984 4458 : inline bool IsValueExactAs([[maybe_unused]] NODATAType noDataValue)
985 : {
986 : if constexpr (std::is_same_v<T, NODATAType>)
987 2753 : return true;
988 : else
989 1705 : return GDALIsValueExactAs<T>(static_cast<double>(noDataValue));
990 : }
991 :
992 : template <bool IS_MAX, class NODATAType>
993 226 : size_t extremum_element(const void *buffer, size_t nElts, GDALDataType eDT,
994 : bool bHasNoData, NODATAType noDataValue)
995 : {
996 226 : switch (eDT)
997 : {
998 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 7, 0)
999 15 : case GDT_Int8:
1000 : {
1001 : using T = int8_t;
1002 15 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1003 15 : return extremum_element<T, IS_MAX>(
1004 : static_cast<const T *>(buffer), nElts, bHasNoData,
1005 15 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1006 : }
1007 : #endif
1008 35 : case GDT_Byte:
1009 : {
1010 : using T = uint8_t;
1011 35 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1012 35 : return extremum_element<T, IS_MAX>(
1013 : static_cast<const T *>(buffer), nElts, bHasNoData,
1014 35 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1015 : }
1016 15 : case GDT_Int16:
1017 : {
1018 : using T = int16_t;
1019 15 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1020 15 : return extremum_element<T, IS_MAX>(
1021 : static_cast<const T *>(buffer), nElts, bHasNoData,
1022 15 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1023 : }
1024 15 : case GDT_UInt16:
1025 : {
1026 : using T = uint16_t;
1027 15 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1028 15 : return extremum_element<T, IS_MAX>(
1029 : static_cast<const T *>(buffer), nElts, bHasNoData,
1030 15 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1031 : }
1032 15 : case GDT_Int32:
1033 : {
1034 : using T = int32_t;
1035 15 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1036 15 : return extremum_element<T, IS_MAX>(
1037 : static_cast<const T *>(buffer), nElts, bHasNoData,
1038 15 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1039 : }
1040 15 : case GDT_UInt32:
1041 : {
1042 : using T = uint32_t;
1043 15 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1044 15 : return extremum_element<T, IS_MAX>(
1045 : static_cast<const T *>(buffer), nElts, bHasNoData,
1046 15 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1047 : }
1048 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 5, 0)
1049 15 : case GDT_Int64:
1050 : {
1051 : using T = int64_t;
1052 15 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1053 15 : return extremum_element<T, IS_MAX>(
1054 : static_cast<const T *>(buffer), nElts, bHasNoData,
1055 15 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1056 : }
1057 15 : case GDT_UInt64:
1058 : {
1059 : using T = uint64_t;
1060 15 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1061 15 : return extremum_element<T, IS_MAX>(
1062 : static_cast<const T *>(buffer), nElts, bHasNoData,
1063 15 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1064 : }
1065 : #endif
1066 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
1067 30 : case GDT_Float16:
1068 : {
1069 : using T = GFloat16;
1070 30 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1071 30 : return extremum_element<T, IS_MAX>(
1072 : static_cast<const T *>(buffer), nElts, bHasNoData,
1073 30 : bHasNoData ? static_cast<T>(noDataValue) : static_cast<T>(0));
1074 : }
1075 : #endif
1076 29 : case GDT_Float32:
1077 : {
1078 : using T = float;
1079 29 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1080 29 : return extremum_element<T, IS_MAX>(
1081 : static_cast<const T *>(buffer), nElts, bHasNoData,
1082 29 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1083 : }
1084 25 : case GDT_Float64:
1085 : {
1086 : using T = double;
1087 25 : bHasNoData = bHasNoData && IsValueExactAs<T>(noDataValue);
1088 25 : return extremum_element<T, IS_MAX>(
1089 : static_cast<const T *>(buffer), nElts, bHasNoData,
1090 25 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1091 : }
1092 2 : case GDT_CInt16:
1093 : case GDT_CInt32:
1094 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
1095 : case GDT_CFloat16:
1096 : #endif
1097 : case GDT_CFloat32:
1098 : case GDT_CFloat64:
1099 : case GDT_Unknown:
1100 : case GDT_TypeCount:
1101 2 : break;
1102 : }
1103 2 : CPLError(CE_Failure, CPLE_NotSupported,
1104 : "%s not supported for this data type.", __FUNCTION__);
1105 2 : return 0;
1106 : }
1107 :
1108 : } // namespace detail
1109 :
1110 : /************************************************************************/
1111 : /* max_element() */
1112 : /************************************************************************/
1113 :
1114 : /** Return the index of the element where the maximum value is hit.
1115 : *
1116 : * If it is hit in several locations, it is not specified which one will be
1117 : * returned.
1118 : *
1119 : * @param buffer Vector of nElts elements of type eDT.
1120 : * @param nElts Number of elements in buffer.
1121 : * @param eDT Data type of the elements of buffer.
1122 : * @param bHasNoData Whether noDataValue is valid.
1123 : * @param noDataValue Nodata value, only taken into account if bHasNoData == true
1124 : *
1125 : * @since GDAL 3.11
1126 : */
1127 : template <class NODATAType>
1128 97 : inline size_t max_element(const void *buffer, size_t nElts, GDALDataType eDT,
1129 : bool bHasNoData, NODATAType noDataValue)
1130 : {
1131 97 : return detail::extremum_element<true>(buffer, nElts, eDT, bHasNoData,
1132 97 : noDataValue);
1133 : }
1134 :
1135 : /************************************************************************/
1136 : /* min_element() */
1137 : /************************************************************************/
1138 :
1139 : /** Return the index of the element where the minimum value is hit.
1140 : *
1141 : * If it is hit in several locations, it is not specified which one will be
1142 : * returned.
1143 : *
1144 : * @param buffer Vector of nElts elements of type eDT.
1145 : * @param nElts Number of elements in buffer.
1146 : * @param eDT Data type of the elements of buffer.
1147 : * @param bHasNoData Whether noDataValue is valid.
1148 : * @param noDataValue Nodata value, only taken into account if bHasNoData == true
1149 : *
1150 : * @since GDAL 3.11
1151 : */
1152 : template <class NODATAType>
1153 129 : inline size_t min_element(const void *buffer, size_t nElts, GDALDataType eDT,
1154 : bool bHasNoData, NODATAType noDataValue)
1155 : {
1156 129 : return detail::extremum_element<false>(buffer, nElts, eDT, bHasNoData,
1157 129 : noDataValue);
1158 : }
1159 :
1160 : namespace detail
1161 : {
1162 :
1163 : #ifdef NOT_EFFICIENT
1164 :
1165 : /************************************************************************/
1166 : /* minmax_element() */
1167 : /************************************************************************/
1168 :
1169 : template <class T>
1170 : std::pair<size_t, size_t> minmax_element(const T *v, size_t size, T noDataValue)
1171 : {
1172 : static_assert(!(std::is_floating_point_v<T>));
1173 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
1174 : static_assert(!std::is_same_v<T, GFloat16>);
1175 : #endif
1176 : if (size == 0)
1177 : return std::pair(0, 0);
1178 : size_t idx_of_min = 0;
1179 : size_t idx_of_max = 0;
1180 : T vmin = v[0];
1181 : T vmax = v[0];
1182 : bool extremum_is_nodata = vmin == noDataValue;
1183 : size_t i = 1;
1184 : for (; i < size; ++i)
1185 : {
1186 : if (v[i] != noDataValue && (v[i] < vmin || extremum_is_nodata))
1187 : {
1188 : vmin = v[i];
1189 : idx_of_min = i;
1190 : extremum_is_nodata = false;
1191 : }
1192 : if (v[i] != noDataValue && (v[i] > vmax || extremum_is_nodata))
1193 : {
1194 : vmax = v[i];
1195 : idx_of_max = i;
1196 : extremum_is_nodata = false;
1197 : }
1198 : }
1199 : return std::pair(idx_of_min, idx_of_max);
1200 : }
1201 :
1202 : template <class T>
1203 : std::pair<size_t, size_t> minmax_element(const T *v, size_t size)
1204 : {
1205 : static_assert(!(std::is_floating_point_v<T>));
1206 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
1207 : static_assert(!std::is_same_v<T, GFloat16>);
1208 : #endif
1209 : if (size == 0)
1210 : return std::pair(0, 0);
1211 : size_t idx_of_min = 0;
1212 : size_t idx_of_max = 0;
1213 : T vmin = v[0];
1214 : T vmax = v[0];
1215 : size_t i = 1;
1216 : for (; i < size; ++i)
1217 : {
1218 : if (v[i] < vmin)
1219 : {
1220 : vmin = v[i];
1221 : idx_of_min = i;
1222 : }
1223 : if (v[i] > vmax)
1224 : {
1225 : vmax = v[i];
1226 : idx_of_max = i;
1227 : }
1228 : }
1229 : return std::pair(idx_of_min, idx_of_max);
1230 : }
1231 :
1232 : template <class T>
1233 : inline std::pair<size_t, size_t> minmax_element_with_nan(const T *v,
1234 : size_t size)
1235 : {
1236 : if (size == 0)
1237 : return std::pair(0, 0);
1238 : size_t idx_of_min = 0;
1239 : size_t idx_of_max = 0;
1240 : T vmin = v[0];
1241 : T vmax = v[0];
1242 : size_t i = 1;
1243 : if (IsNan(v[0]))
1244 : {
1245 : for (; i < size; ++i)
1246 : {
1247 : if (!IsNan(v[i]))
1248 : {
1249 : vmin = v[i];
1250 : idx_of_min = i;
1251 : vmax = v[i];
1252 : idx_of_max = i;
1253 : break;
1254 : }
1255 : }
1256 : }
1257 : for (; i < size; ++i)
1258 : {
1259 : if (v[i] < vmin)
1260 : {
1261 : vmin = v[i];
1262 : idx_of_min = i;
1263 : }
1264 : if (v[i] > vmax)
1265 : {
1266 : vmax = v[i];
1267 : idx_of_max = i;
1268 : }
1269 : }
1270 : return std::pair(idx_of_min, idx_of_max);
1271 : }
1272 :
1273 : template <>
1274 : std::pair<size_t, size_t> minmax_element<float>(const float *v, size_t size)
1275 : {
1276 : return minmax_element_with_nan<float>(v, size);
1277 : }
1278 :
1279 : template <>
1280 : std::pair<size_t, size_t> minmax_element<double>(const double *v, size_t size)
1281 : {
1282 : return minmax_element_with_nan<double>(v, size);
1283 : }
1284 :
1285 : template <class T>
1286 : inline std::pair<size_t, size_t> minmax_element(const T *buffer, size_t size,
1287 : bool bHasNoData, T noDataValue)
1288 : {
1289 : if (bHasNoData)
1290 : {
1291 : return minmax_element<T>(buffer, size, noDataValue);
1292 : }
1293 : else
1294 : {
1295 : return minmax_element<T>(buffer, size);
1296 : }
1297 : }
1298 : #else
1299 :
1300 : /************************************************************************/
1301 : /* minmax_element() */
1302 : /************************************************************************/
1303 :
1304 : template <class T>
1305 6869 : inline std::pair<size_t, size_t> minmax_element(const T *buffer, size_t size,
1306 : bool bHasNoData, T noDataValue)
1307 : {
1308 : #ifdef NOT_EFFICIENT
1309 : if (bHasNoData)
1310 : {
1311 : return minmax_element<T>(buffer, size, noDataValue);
1312 : }
1313 : else
1314 : {
1315 : return minmax_element<T>(buffer, size);
1316 : //auto [imin, imax] = std::minmax_element(buffer, buffer + size);
1317 : //return std::pair(imin - buffer, imax - buffer);
1318 : }
1319 : #else
1320 :
1321 : #if !defined(GDAL_MINMAX_ELEMENT_USE_SSE2)
1322 : if constexpr (!std::is_floating_point_v<T>
1323 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
1324 : && !std::is_same_v<T, GFloat16>
1325 : #endif
1326 : )
1327 : {
1328 : if (!bHasNoData)
1329 : {
1330 : auto [min_iter, max_iter] =
1331 : std::minmax_element(buffer, buffer + size);
1332 : return std::pair(min_iter - buffer, max_iter - buffer);
1333 : }
1334 : }
1335 : #endif
1336 :
1337 : // Using separately min and max is more efficient than computing them
1338 : // within the same loop
1339 : return std::pair(
1340 20607 : extremum_element<T, false>(buffer, size, bHasNoData, noDataValue),
1341 6869 : extremum_element<T, true>(buffer, size, bHasNoData, noDataValue));
1342 :
1343 : #endif
1344 : }
1345 : #endif
1346 :
1347 : } // namespace detail
1348 :
1349 : /************************************************************************/
1350 : /* minmax_element() */
1351 : /************************************************************************/
1352 :
1353 : /** Return the index of the elements where the minimum and maximum values are hit.
1354 : *
1355 : * If they are hit in several locations, it is not specified which one will be
1356 : * returned (contrary to std::minmax_element).
1357 : *
1358 : * @param buffer Vector of nElts elements of type eDT.
1359 : * @param nElts Number of elements in buffer.
1360 : * @param eDT Data type of the elements of buffer.
1361 : * @param bHasNoData Whether noDataValue is valid.
1362 : * @param noDataValue Nodata value, only taken into account if bHasNoData == true
1363 : *
1364 : * @since GDAL 3.11
1365 : */
1366 : template <class NODATAType>
1367 : inline std::pair<size_t, size_t>
1368 6870 : minmax_element(const void *buffer, size_t nElts, GDALDataType eDT,
1369 : bool bHasNoData, NODATAType noDataValue)
1370 : {
1371 6870 : switch (eDT)
1372 : {
1373 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 7, 0)
1374 106 : case GDT_Int8:
1375 : {
1376 : using T = int8_t;
1377 106 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1378 106 : return detail::minmax_element<T>(
1379 : static_cast<const T *>(buffer), nElts, bHasNoData,
1380 106 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1381 : }
1382 : #endif
1383 5 : case GDT_Byte:
1384 : {
1385 : using T = uint8_t;
1386 5 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1387 5 : return detail::minmax_element<T>(
1388 : static_cast<const T *>(buffer), nElts, bHasNoData,
1389 5 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1390 : }
1391 1399 : case GDT_Int16:
1392 : {
1393 : using T = int16_t;
1394 1399 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1395 1399 : return detail::minmax_element<T>(
1396 : static_cast<const T *>(buffer), nElts, bHasNoData,
1397 1399 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1398 : }
1399 1 : case GDT_UInt16:
1400 : {
1401 : using T = uint16_t;
1402 1 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1403 1 : return detail::minmax_element<T>(
1404 : static_cast<const T *>(buffer), nElts, bHasNoData,
1405 1 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1406 : }
1407 1068 : case GDT_Int32:
1408 : {
1409 : using T = int32_t;
1410 1068 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1411 1068 : return detail::minmax_element<T>(
1412 : static_cast<const T *>(buffer), nElts, bHasNoData,
1413 1068 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1414 : }
1415 177 : case GDT_UInt32:
1416 : {
1417 : using T = uint32_t;
1418 177 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1419 177 : return detail::minmax_element<T>(
1420 : static_cast<const T *>(buffer), nElts, bHasNoData,
1421 177 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1422 : }
1423 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 5, 0)
1424 21 : case GDT_Int64:
1425 : {
1426 : using T = int64_t;
1427 21 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1428 21 : return detail::minmax_element<T>(
1429 : static_cast<const T *>(buffer), nElts, bHasNoData,
1430 21 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1431 : }
1432 17 : case GDT_UInt64:
1433 : {
1434 : using T = uint64_t;
1435 17 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1436 17 : return detail::minmax_element<T>(
1437 : static_cast<const T *>(buffer), nElts, bHasNoData,
1438 17 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1439 : }
1440 : #endif
1441 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
1442 5 : case GDT_Float16:
1443 : {
1444 : using T = GFloat16;
1445 5 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1446 5 : return detail::minmax_element<T>(
1447 : static_cast<const T *>(buffer), nElts, bHasNoData,
1448 5 : bHasNoData ? static_cast<T>(noDataValue) : static_cast<T>(0));
1449 : }
1450 : #endif
1451 634 : case GDT_Float32:
1452 : {
1453 : using T = float;
1454 634 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1455 634 : return detail::minmax_element<T>(
1456 : static_cast<const T *>(buffer), nElts, bHasNoData,
1457 634 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1458 : }
1459 3436 : case GDT_Float64:
1460 : {
1461 : using T = double;
1462 3436 : bHasNoData = bHasNoData && detail::IsValueExactAs<T>(noDataValue);
1463 3436 : return detail::minmax_element<T>(
1464 : static_cast<const T *>(buffer), nElts, bHasNoData,
1465 3436 : bHasNoData ? static_cast<T>(noDataValue) : 0);
1466 : }
1467 1 : case GDT_CInt16:
1468 : case GDT_CInt32:
1469 : #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 11, 0)
1470 : case GDT_CFloat16:
1471 : #endif
1472 : case GDT_CFloat32:
1473 : case GDT_CFloat64:
1474 : case GDT_Unknown:
1475 : case GDT_TypeCount:
1476 1 : break;
1477 : }
1478 1 : CPLError(CE_Failure, CPLE_NotSupported,
1479 : "%s not supported for this data type.", __FUNCTION__);
1480 1 : return std::pair(0, 0);
1481 : }
1482 :
1483 : } // namespace GDAL_MINMAXELT_NS
1484 :
1485 : #endif // GDAL_MINMAX_ELEMENT_INCLUDED
|
