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