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