Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL Core
4 : * Purpose: Inline C++ templates
5 : * Author: Phil Vachon, <philippe at cowpig.ca>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2009, Phil Vachon, <philippe at cowpig.ca>
9 : * Copyright (c) 2025, Even Rouault, <even.rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #ifndef GDAL_PRIV_TEMPLATES_HPP_INCLUDED
15 : #define GDAL_PRIV_TEMPLATES_HPP_INCLUDED
16 :
17 : #include "cpl_port.h"
18 :
19 : #include <algorithm>
20 : #include <cmath>
21 : #include <cstdint>
22 : #include <limits>
23 : #include <type_traits>
24 :
25 : #include "cpl_float.h"
26 :
27 : // Needs SSE2
28 : #if defined(__x86_64) || defined(_M_X64) || defined(USE_SSE2) || \
29 : defined(USE_NEON_OPTIMIZATIONS)
30 :
31 : #ifdef USE_NEON_OPTIMIZATIONS
32 : #include "include_sse2neon.h"
33 : #else
34 : #include <immintrin.h>
35 : #endif
36 :
37 44330595 : static inline void GDALCopyXMMToInt32(const __m128i xmm, void *pDest)
38 : {
39 44330595 : int n32 = _mm_cvtsi128_si32(xmm); // Extract lower 32 bit word
40 44330595 : memcpy(pDest, &n32, sizeof(n32));
41 44330595 : }
42 :
43 94076911 : static inline void GDALCopyXMMToInt64(const __m128i xmm, void *pDest)
44 : {
45 : _mm_storel_epi64(reinterpret_cast<__m128i *>(pDest), xmm);
46 94076911 : }
47 :
48 : #if __SSSE3__
49 : #include <tmmintrin.h>
50 : #endif
51 :
52 : #if defined(__SSE4_1__) || defined(__AVX__)
53 : #include <smmintrin.h>
54 : #endif
55 :
56 : #ifdef __F16C__
57 : #include <immintrin.h>
58 : #endif
59 :
60 : #endif
61 :
62 : /************************************************************************/
63 : /* GDALGetDataLimits() */
64 : /************************************************************************/
65 : /**
66 : * Compute the limits of values that can be placed in Tout in terms of
67 : * Tin. Usually used for output clamping, when the output data type's
68 : * limits are stable relative to the input type (i.e. no roundoff error).
69 : *
70 : * @param tMaxValue the returned maximum value
71 : * @param tMinValue the returned minimum value
72 : */
73 :
74 : template <class Tin, class Tout>
75 367022900 : inline void GDALGetDataLimits(Tin &tMaxValue, Tin &tMinValue)
76 : {
77 367022900 : tMaxValue = cpl::NumericLimits<Tin>::max();
78 366781590 : tMinValue = cpl::NumericLimits<Tin>::lowest();
79 :
80 : // Compute the actual minimum value of Tout in terms of Tin.
81 : if constexpr (cpl::NumericLimits<Tout>::is_signed &&
82 : cpl::NumericLimits<Tout>::is_integer)
83 : {
84 : // the minimum value is less than zero
85 : // cppcheck-suppress knownConditionTrueFalse
86 : if constexpr (cpl::NumericLimits<Tout>::digits <
87 : cpl::NumericLimits<Tin>::digits ||
88 : !cpl::NumericLimits<Tin>::is_integer)
89 : {
90 : // Tout is smaller than Tin, so we need to clamp values in input
91 : // to the range of Tout's min/max values
92 : if constexpr (cpl::NumericLimits<Tin>::is_signed)
93 : {
94 85290070 : tMinValue =
95 85290070 : static_cast<Tin>(cpl::NumericLimits<Tout>::lowest());
96 : }
97 85810726 : tMaxValue = static_cast<Tin>(cpl::NumericLimits<Tout>::max());
98 : }
99 : }
100 : else if constexpr (cpl::NumericLimits<Tout>::is_integer)
101 : {
102 : // the output is unsigned, so we just need to determine the max
103 : if constexpr (!std::is_same_v<Tin, Tout> &&
104 : cpl::NumericLimits<Tout>::digits <=
105 : cpl::NumericLimits<Tin>::digits)
106 : {
107 : // Tout is smaller than Tin, so we need to clamp the input values
108 : // to the range of Tout's max
109 133035163 : tMaxValue = static_cast<Tin>(cpl::NumericLimits<Tout>::max());
110 : }
111 134227090 : tMinValue = 0;
112 : }
113 366414620 : }
114 :
115 : /************************************************************************/
116 : /* GDALClampValue() */
117 : /************************************************************************/
118 : /**
119 : * Clamp values of type T to a specified range
120 : *
121 : * @param tValue the value
122 : * @param tMax the max value
123 : * @param tMin the min value
124 : */
125 : template <class T>
126 365887433 : inline T GDALClampValue(const T tValue, const T tMax, const T tMin)
127 : {
128 365887433 : return tValue > tMax ? tMax : tValue < tMin ? tMin : tValue;
129 : }
130 :
131 : /************************************************************************/
132 : /* GDALClampDoubleValue() */
133 : /************************************************************************/
134 : /**
135 : * Clamp double values to a specified range, this uses the same
136 : * argument ordering as std::clamp, returns TRUE if the value was clamped.
137 : *
138 : * @param tValue the value
139 : * @param tMin the min value
140 : * @param tMax the max value
141 : *
142 : */
143 : template <class T2, class T3>
144 269 : inline bool GDALClampDoubleValue(double &tValue, const T2 tMin, const T3 tMax)
145 : {
146 269 : const double tMin2{static_cast<double>(tMin)};
147 269 : const double tMax2{static_cast<double>(tMax)};
148 269 : if (tValue > tMax2 || tValue < tMin2)
149 : {
150 26 : tValue = tValue > tMax2 ? tMax2 : tValue < tMin2 ? tMin2 : tValue;
151 26 : return true;
152 : }
153 : else
154 : {
155 243 : return false;
156 : }
157 : }
158 :
159 : /************************************************************************/
160 : /* GDALIsValueInRange() */
161 : /************************************************************************/
162 : /**
163 : * Returns whether a value is in the type range.
164 : * NaN is considered not to be in type range.
165 : *
166 : * @param dfValue the value
167 : * @return whether the value is in the type range.
168 : */
169 152044 : template <class T> inline bool GDALIsValueInRange(double dfValue)
170 : {
171 304017 : return dfValue >= static_cast<double>(cpl::NumericLimits<T>::lowest()) &&
172 304005 : dfValue <= static_cast<double>(cpl::NumericLimits<T>::max());
173 : }
174 :
175 26 : template <> inline bool GDALIsValueInRange<double>(double dfValue)
176 : {
177 26 : return !CPLIsNan(dfValue);
178 : }
179 :
180 42482 : template <> inline bool GDALIsValueInRange<float>(double dfValue)
181 : {
182 84887 : return CPLIsInf(dfValue) ||
183 42405 : (dfValue >= -std::numeric_limits<float>::max() &&
184 84874 : dfValue <= std::numeric_limits<float>::max());
185 : }
186 :
187 1 : template <> inline bool GDALIsValueInRange<GFloat16>(double dfValue)
188 : {
189 3 : return CPLIsInf(dfValue) ||
190 2 : (dfValue >= -cpl::NumericLimits<GFloat16>::max() &&
191 2 : dfValue <= cpl::NumericLimits<GFloat16>::max());
192 : }
193 :
194 16899 : template <> inline bool GDALIsValueInRange<int64_t>(double dfValue)
195 : {
196 : // Values in the range [INT64_MAX - 1023, INT64_MAX - 1]
197 : // get converted to a double that once cast to int64_t is
198 : // INT64_MAX + 1, hence the < strict comparison.
199 : return dfValue >=
200 33797 : static_cast<double>(cpl::NumericLimits<int64_t>::lowest()) &&
201 33797 : dfValue < static_cast<double>(cpl::NumericLimits<int64_t>::max());
202 : }
203 :
204 8914 : template <> inline bool GDALIsValueInRange<uint64_t>(double dfValue)
205 : {
206 : // Values in the range [UINT64_MAX - 2047, UINT64_MAX - 1]
207 : // get converted to a double that once cast to uint64_t is
208 : // UINT64_MAX + 1, hence the < strict comparison.
209 17825 : return dfValue >= 0 &&
210 17825 : dfValue < static_cast<double>(cpl::NumericLimits<uint64_t>::max());
211 : }
212 :
213 : /************************************************************************/
214 : /* GDALIsValueExactAs() */
215 : /************************************************************************/
216 : /**
217 : * Returns whether a value can be exactly represented on type T.
218 : *
219 : * That is static_cast\<double\>(static_cast\<T\>(dfValue)) is legal and is
220 : * equal to dfValue.
221 : *
222 : * Note: for T=float or double, a NaN input leads to true
223 : *
224 : * @param dfValue the value
225 : * @return whether the value can be exactly represented on type T.
226 : */
227 5344 : template <class T> inline bool GDALIsValueExactAs(double dfValue)
228 : {
229 10628 : return GDALIsValueInRange<T>(dfValue) &&
230 10618 : static_cast<double>(static_cast<T>(dfValue)) == dfValue;
231 : }
232 :
233 146 : template <> inline bool GDALIsValueExactAs<float>(double dfValue)
234 : {
235 284 : return CPLIsNan(dfValue) ||
236 138 : (GDALIsValueInRange<float>(dfValue) &&
237 280 : static_cast<double>(static_cast<float>(dfValue)) == dfValue);
238 : }
239 :
240 0 : template <> inline bool GDALIsValueExactAs<GFloat16>(double dfValue)
241 : {
242 0 : return CPLIsNan(dfValue) ||
243 0 : (GDALIsValueInRange<GFloat16>(dfValue) &&
244 0 : static_cast<double>(static_cast<GFloat16>(dfValue)) == dfValue);
245 : }
246 :
247 16 : template <> inline bool GDALIsValueExactAs<double>(double)
248 : {
249 16 : return true;
250 : }
251 :
252 : /************************************************************************/
253 : /* GDALCopyWord() */
254 : /************************************************************************/
255 :
256 : // Integer input and output: clamp the input
257 :
258 : template <class Tin, class Tout> struct sGDALCopyWord
259 : {
260 205990750 : static inline void f(const Tin tValueIn, Tout &tValueOut)
261 : {
262 : Tin tMaxVal, tMinVal;
263 205990750 : GDALGetDataLimits<Tin, Tout>(tMaxVal, tMinVal);
264 206659650 : tValueOut =
265 205739550 : static_cast<Tout>(GDALClampValue(tValueIn, tMaxVal, tMinVal));
266 206659650 : }
267 : };
268 :
269 : // Integer input and floating point output: simply convert
270 :
271 : template <class Tin> struct sGDALCopyWord<Tin, GFloat16>
272 : {
273 270264 : static inline void f(const Tin tValueIn, GFloat16 &hfValueOut)
274 : {
275 270264 : hfValueOut = static_cast<GFloat16>(tValueIn);
276 270264 : }
277 : };
278 :
279 : template <class Tin> struct sGDALCopyWord<Tin, float>
280 : {
281 7588275 : static inline void f(const Tin tValueIn, float &fValueOut)
282 : {
283 7588275 : fValueOut = static_cast<float>(tValueIn);
284 7588275 : }
285 : };
286 :
287 : template <class Tin> struct sGDALCopyWord<Tin, double>
288 : {
289 71813575 : static inline void f(const Tin tValueIn, double &dfValueOut)
290 : {
291 71813575 : dfValueOut = static_cast<double>(tValueIn);
292 71813575 : }
293 : };
294 :
295 : // Floating point input and output, converting between identical types: simply copy
296 :
297 : template <> struct sGDALCopyWord<GFloat16, GFloat16>
298 : {
299 : static inline void f(const GFloat16 hfValueIn, GFloat16 &hfValueOut)
300 : {
301 : hfValueOut = hfValueIn;
302 : }
303 : };
304 :
305 : template <> struct sGDALCopyWord<float, float>
306 : {
307 : static inline void f(const float fValueIn, float &fValueOut)
308 : {
309 : fValueOut = fValueIn;
310 : }
311 : };
312 :
313 : template <> struct sGDALCopyWord<double, double>
314 : {
315 : static inline void f(const double dfValueIn, double &dfValueOut)
316 : {
317 : dfValueOut = dfValueIn;
318 : }
319 : };
320 :
321 : // Floating point input and output, converting to a larger type: use implicit conversion
322 :
323 : template <> struct sGDALCopyWord<GFloat16, float>
324 : {
325 6189 : static inline void f(const GFloat16 hfValueIn, float &dfValueOut)
326 : {
327 6189 : dfValueOut = hfValueIn;
328 6189 : }
329 : };
330 :
331 : template <> struct sGDALCopyWord<GFloat16, double>
332 : {
333 2530 : static inline void f(const GFloat16 hfValueIn, double &dfValueOut)
334 : {
335 2530 : dfValueOut = hfValueIn;
336 2530 : }
337 : };
338 :
339 : template <> struct sGDALCopyWord<float, double>
340 : {
341 81515 : static inline void f(const float fValueIn, double &dfValueOut)
342 : {
343 81515 : dfValueOut = fValueIn;
344 81515 : }
345 : };
346 :
347 : // Floating point input and out, converting to a smaller type: ensure overflow results in infinity
348 :
349 : template <> struct sGDALCopyWord<float, GFloat16>
350 : {
351 612 : static inline void f(const float fValueIn, GFloat16 &hfValueOut)
352 : {
353 : // Our custom implementation when std::float16_t is not
354 : // available ensures proper behavior.
355 : #if !defined(HAVE_STD_FLOAT16_T)
356 612 : if (fValueIn > cpl::NumericLimits<GFloat16>::max())
357 : {
358 69 : hfValueOut = cpl::NumericLimits<GFloat16>::infinity();
359 69 : return;
360 : }
361 543 : if (fValueIn < -cpl::NumericLimits<GFloat16>::max())
362 : {
363 68 : hfValueOut = -cpl::NumericLimits<GFloat16>::infinity();
364 68 : return;
365 : }
366 : #endif
367 475 : hfValueOut = static_cast<GFloat16>(fValueIn);
368 : }
369 : };
370 :
371 : template <> struct sGDALCopyWord<double, GFloat16>
372 : {
373 2216 : static inline void f(const double dfValueIn, GFloat16 &hfValueOut)
374 : {
375 : // Our custom implementation when std::float16_t is not
376 : // available ensures proper behavior.
377 : #if !defined(HAVE_STD_FLOAT16_T)
378 2216 : if (dfValueIn > cpl::NumericLimits<GFloat16>::max())
379 : {
380 69 : hfValueOut = cpl::NumericLimits<GFloat16>::infinity();
381 69 : return;
382 : }
383 2147 : if (dfValueIn < -cpl::NumericLimits<GFloat16>::max())
384 : {
385 69 : hfValueOut = -cpl::NumericLimits<GFloat16>::infinity();
386 69 : return;
387 : }
388 : #endif
389 2078 : hfValueOut = static_cast<GFloat16>(dfValueIn);
390 : }
391 : };
392 :
393 : template <> struct sGDALCopyWord<double, float>
394 : {
395 129370 : static inline void f(const double dfValueIn, float &fValueOut)
396 : {
397 : #if defined(__x86_64) || defined(_M_X64) || defined(USE_SSE2)
398 : // We could just write fValueOut = static_cast<float>(dfValueIn);
399 : // but a sanitizer might complain with values above FLT_MAX
400 388110 : _mm_store_ss(&fValueOut,
401 : _mm_cvtsd_ss(_mm_undefined_ps(), _mm_load_sd(&dfValueIn)));
402 : #else
403 : if (dfValueIn > std::numeric_limits<float>::max())
404 : {
405 : fValueOut = std::numeric_limits<float>::infinity();
406 : return;
407 : }
408 : if (dfValueIn < -std::numeric_limits<float>::max())
409 : {
410 : fValueOut = -std::numeric_limits<float>::infinity();
411 : return;
412 : }
413 :
414 : fValueOut = static_cast<float>(dfValueIn);
415 : #endif
416 129370 : }
417 : };
418 :
419 : // Floating point input to a small unsigned integer type: nan becomes zero, otherwise round and clamp
420 :
421 : template <class Tout> struct sGDALCopyWord<GFloat16, Tout>
422 : {
423 5032 : static inline void f(const GFloat16 hfValueIn, Tout &tValueOut)
424 : {
425 5032 : if (CPLIsNan(hfValueIn))
426 : {
427 1 : tValueOut = 0;
428 1 : return;
429 : }
430 : GFloat16 hfMaxVal, hfMinVal;
431 5031 : GDALGetDataLimits<GFloat16, Tout>(hfMaxVal, hfMinVal);
432 5031 : tValueOut = static_cast<Tout>(
433 10062 : GDALClampValue(hfValueIn + GFloat16(0.5f), hfMaxVal, hfMinVal));
434 : }
435 : };
436 :
437 : template <class Tout> struct sGDALCopyWord<float, Tout>
438 : {
439 4485240 : static inline void f(const float fValueIn, Tout &tValueOut)
440 : {
441 4485240 : if (CPLIsNan(fValueIn))
442 : {
443 0 : tValueOut = 0;
444 0 : return;
445 : }
446 : float fMaxVal, fMinVal;
447 4485240 : GDALGetDataLimits<float, Tout>(fMaxVal, fMinVal);
448 4485150 : tValueOut = static_cast<Tout>(
449 4485210 : GDALClampValue(fValueIn + 0.5f, fMaxVal, fMinVal));
450 : }
451 : };
452 :
453 : template <class Tout> struct sGDALCopyWord<double, Tout>
454 : {
455 71653110 : static inline void f(const double dfValueIn, Tout &tValueOut)
456 : {
457 71653110 : if (CPLIsNan(dfValueIn))
458 : {
459 6 : tValueOut = 0;
460 6 : return;
461 : }
462 : double dfMaxVal, dfMinVal;
463 71823309 : GDALGetDataLimits<double, Tout>(dfMaxVal, dfMinVal);
464 71127709 : tValueOut = static_cast<Tout>(
465 70340709 : GDALClampValue(dfValueIn + 0.5, dfMaxVal, dfMinVal));
466 : }
467 : };
468 :
469 : // Floating point input to a large unsigned integer type: nan becomes zero, otherwise round and clamp.
470 : // Avoid roundoff while clamping.
471 :
472 : template <> struct sGDALCopyWord<GFloat16, std::uint64_t>
473 : {
474 624 : static inline void f(const GFloat16 hfValueIn, std::uint64_t &nValueOut)
475 : {
476 624 : if (!(hfValueIn > 0))
477 : {
478 4 : nValueOut = 0;
479 : }
480 620 : else if (CPLIsInf(hfValueIn))
481 : {
482 1 : nValueOut = cpl::NumericLimits<std::uint64_t>::max();
483 : }
484 : else
485 : {
486 619 : nValueOut = static_cast<std::uint64_t>(hfValueIn + GFloat16(0.5f));
487 : }
488 624 : }
489 : };
490 :
491 : template <> struct sGDALCopyWord<float, unsigned int>
492 : {
493 5039 : static inline void f(const float fValueIn, unsigned int &nValueOut)
494 : {
495 5039 : if (!(fValueIn > 0))
496 : {
497 149 : nValueOut = 0;
498 : }
499 4890 : else if (fValueIn >=
500 4890 : static_cast<float>(cpl::NumericLimits<unsigned int>::max()))
501 : {
502 134 : nValueOut = cpl::NumericLimits<unsigned int>::max();
503 : }
504 : else
505 : {
506 4756 : nValueOut = static_cast<unsigned int>(fValueIn + 0.5f);
507 : }
508 5039 : }
509 : };
510 :
511 : template <> struct sGDALCopyWord<float, std::uint64_t>
512 : {
513 624 : static inline void f(const float fValueIn, std::uint64_t &nValueOut)
514 : {
515 624 : if (!(fValueIn > 0))
516 : {
517 4 : nValueOut = 0;
518 : }
519 620 : else if (fValueIn >=
520 620 : static_cast<float>(cpl::NumericLimits<std::uint64_t>::max()))
521 : {
522 2 : nValueOut = cpl::NumericLimits<std::uint64_t>::max();
523 : }
524 : else
525 : {
526 618 : nValueOut = static_cast<std::uint64_t>(fValueIn + 0.5f);
527 : }
528 624 : }
529 : };
530 :
531 : template <> struct sGDALCopyWord<double, std::uint64_t>
532 : {
533 1062 : static inline void f(const double dfValueIn, std::uint64_t &nValueOut)
534 : {
535 1062 : if (!(dfValueIn > 0))
536 : {
537 166 : nValueOut = 0;
538 : }
539 896 : else if (dfValueIn >
540 896 : static_cast<double>(cpl::NumericLimits<uint64_t>::max()))
541 : {
542 4 : nValueOut = cpl::NumericLimits<uint64_t>::max();
543 : }
544 : else
545 : {
546 892 : nValueOut = static_cast<std::uint64_t>(dfValueIn + 0.5);
547 : }
548 1062 : }
549 : };
550 :
551 : // Floating point input to a very large unsigned integer type: nan becomes zero, otherwise round and clamp.
552 : // Avoid infinity while clamping when the maximum integer is too large for the floating-point type.
553 : // Avoid roundoff while clamping.
554 :
555 : template <> struct sGDALCopyWord<GFloat16, unsigned short>
556 : {
557 4706 : static inline void f(const GFloat16 hfValueIn, unsigned short &nValueOut)
558 : {
559 4706 : if (!(hfValueIn > 0))
560 : {
561 216 : nValueOut = 0;
562 : }
563 4490 : else if (CPLIsInf(hfValueIn))
564 : {
565 67 : nValueOut = cpl::NumericLimits<unsigned short>::max();
566 : }
567 : else
568 : {
569 4423 : nValueOut = static_cast<unsigned short>(hfValueIn + GFloat16(0.5f));
570 : }
571 4706 : }
572 : };
573 :
574 : template <> struct sGDALCopyWord<GFloat16, unsigned int>
575 : {
576 4639 : static inline void f(const GFloat16 hfValueIn, unsigned int &nValueOut)
577 : {
578 4639 : if (!(hfValueIn > 0))
579 : {
580 149 : nValueOut = 0;
581 : }
582 4490 : else if (CPLIsInf(hfValueIn))
583 : {
584 0 : nValueOut = cpl::NumericLimits<unsigned int>::max();
585 : }
586 : else
587 : {
588 4490 : nValueOut = static_cast<unsigned int>(hfValueIn + GFloat16(0.5f));
589 : }
590 4639 : }
591 : };
592 :
593 : // Floating point input to a small signed integer type: nan becomes zero, otherwise round and clamp.
594 : // Rounding for signed integers is different than for the unsigned integers above.
595 :
596 : template <> struct sGDALCopyWord<GFloat16, signed char>
597 : {
598 1088 : static inline void f(const GFloat16 hfValueIn, signed char &nValueOut)
599 : {
600 1088 : if (CPLIsNan(hfValueIn))
601 : {
602 0 : nValueOut = 0;
603 0 : return;
604 : }
605 : GFloat16 hfMaxVal, hfMinVal;
606 1088 : GDALGetDataLimits<GFloat16, signed char>(hfMaxVal, hfMinVal);
607 1088 : GFloat16 hfValue = hfValueIn >= GFloat16(0.0f)
608 619 : ? hfValueIn + GFloat16(0.5f)
609 1088 : : hfValueIn - GFloat16(0.5f);
610 1088 : nValueOut = static_cast<signed char>(
611 2176 : GDALClampValue(hfValue, hfMaxVal, hfMinVal));
612 : }
613 : };
614 :
615 : template <> struct sGDALCopyWord<float, signed char>
616 : {
617 1172 : static inline void f(const float fValueIn, signed char &nValueOut)
618 : {
619 1172 : if (CPLIsNan(fValueIn))
620 : {
621 0 : nValueOut = 0;
622 0 : return;
623 : }
624 : float fMaxVal, fMinVal;
625 1172 : GDALGetDataLimits<float, signed char>(fMaxVal, fMinVal);
626 1172 : float fValue = fValueIn >= 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f;
627 1172 : nValueOut =
628 1172 : static_cast<signed char>(GDALClampValue(fValue, fMaxVal, fMinVal));
629 : }
630 : };
631 :
632 : template <> struct sGDALCopyWord<float, short>
633 : {
634 2939800 : static inline void f(const float fValueIn, short &nValueOut)
635 : {
636 2939800 : if (CPLIsNan(fValueIn))
637 : {
638 0 : nValueOut = 0;
639 0 : return;
640 : }
641 : float fMaxVal, fMinVal;
642 2939800 : GDALGetDataLimits<float, short>(fMaxVal, fMinVal);
643 2939800 : float fValue = fValueIn >= 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f;
644 2939800 : nValueOut =
645 2939800 : static_cast<short>(GDALClampValue(fValue, fMaxVal, fMinVal));
646 : }
647 : };
648 :
649 : template <> struct sGDALCopyWord<double, signed char>
650 : {
651 1320 : static inline void f(const double dfValueIn, signed char &nValueOut)
652 : {
653 1320 : if (CPLIsNan(dfValueIn))
654 : {
655 0 : nValueOut = 0;
656 0 : return;
657 : }
658 : double dfMaxVal, dfMinVal;
659 1320 : GDALGetDataLimits<double, signed char>(dfMaxVal, dfMinVal);
660 1320 : double dfValue = dfValueIn > 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5;
661 1320 : nValueOut = static_cast<signed char>(
662 1320 : GDALClampValue(dfValue, dfMaxVal, dfMinVal));
663 : }
664 : };
665 :
666 : template <> struct sGDALCopyWord<double, short>
667 : {
668 5156900 : static inline void f(const double dfValueIn, short &nValueOut)
669 : {
670 5156900 : if (CPLIsNan(dfValueIn))
671 : {
672 2 : nValueOut = 0;
673 2 : return;
674 : }
675 : double dfMaxVal, dfMinVal;
676 5156890 : GDALGetDataLimits<double, short>(dfMaxVal, dfMinVal);
677 5156870 : double dfValue = dfValueIn > 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5;
678 5156880 : nValueOut =
679 5156870 : static_cast<short>(GDALClampValue(dfValue, dfMaxVal, dfMinVal));
680 : }
681 : };
682 :
683 : template <> struct sGDALCopyWord<double, int>
684 : {
685 77098100 : static inline void f(const double dfValueIn, int &nValueOut)
686 : {
687 77098100 : if (CPLIsNan(dfValueIn))
688 : {
689 0 : nValueOut = 0;
690 0 : return;
691 : }
692 : double dfMaxVal, dfMinVal;
693 77098100 : GDALGetDataLimits<double, int>(dfMaxVal, dfMinVal);
694 77098100 : double dfValue = dfValueIn >= 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5;
695 77098100 : nValueOut =
696 77098100 : static_cast<int>(GDALClampValue(dfValue, dfMaxVal, dfMinVal));
697 : }
698 : };
699 :
700 : // Floating point input to a large signed integer type: nan becomes zero, otherwise round and clamp.
701 : // Rounding for signed integers is different than for the unsigned integers above.
702 : // Avoid roundoff while clamping.
703 :
704 : template <> struct sGDALCopyWord<GFloat16, short>
705 : {
706 5289 : static inline void f(const GFloat16 hfValueIn, short &nValueOut)
707 : {
708 5289 : if (CPLIsNan(hfValueIn))
709 : {
710 0 : nValueOut = 0;
711 : }
712 5289 : else if (hfValueIn >=
713 5289 : static_cast<GFloat16>(cpl::NumericLimits<short>::max()))
714 : {
715 146 : nValueOut = cpl::NumericLimits<short>::max();
716 : }
717 5143 : else if (hfValueIn <=
718 5143 : static_cast<GFloat16>(cpl::NumericLimits<short>::lowest()))
719 : {
720 213 : nValueOut = cpl::NumericLimits<short>::lowest();
721 : }
722 : else
723 : {
724 9860 : nValueOut = static_cast<short>(hfValueIn > GFloat16(0.0f)
725 4930 : ? hfValueIn + GFloat16(0.5f)
726 5451 : : hfValueIn - GFloat16(0.5f));
727 : }
728 5289 : }
729 : };
730 :
731 : template <> struct sGDALCopyWord<float, int>
732 : {
733 150877000 : static inline void f(const float fValueIn, int &nValueOut)
734 : {
735 150877000 : if (CPLIsNan(fValueIn))
736 : {
737 0 : nValueOut = 0;
738 : }
739 153578000 : else if (fValueIn >= static_cast<float>(cpl::NumericLimits<int>::max()))
740 : {
741 274 : nValueOut = cpl::NumericLimits<int>::max();
742 : }
743 150267000 : else if (fValueIn <=
744 153172000 : static_cast<float>(cpl::NumericLimits<int>::lowest()))
745 : {
746 0 : nValueOut = cpl::NumericLimits<int>::lowest();
747 : }
748 : else
749 : {
750 153099000 : nValueOut = static_cast<int>(fValueIn > 0.0f ? fValueIn + 0.5f
751 121386 : : fValueIn - 0.5f);
752 : }
753 152978000 : }
754 : };
755 :
756 : template <> struct sGDALCopyWord<float, std::int64_t>
757 : {
758 1093 : static inline void f(const float fValueIn, std::int64_t &nValueOut)
759 : {
760 1093 : if (CPLIsNan(fValueIn))
761 : {
762 1 : nValueOut = 0;
763 : }
764 1092 : else if (fValueIn >=
765 1092 : static_cast<float>(cpl::NumericLimits<std::int64_t>::max()))
766 : {
767 2 : nValueOut = cpl::NumericLimits<std::int64_t>::max();
768 : }
769 1090 : else if (fValueIn <=
770 1090 : static_cast<float>(cpl::NumericLimits<std::int64_t>::lowest()))
771 : {
772 2 : nValueOut = cpl::NumericLimits<std::int64_t>::lowest();
773 : }
774 : else
775 : {
776 2176 : nValueOut = static_cast<std::int64_t>(
777 1088 : fValueIn > 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f);
778 : }
779 1093 : }
780 : };
781 :
782 : template <> struct sGDALCopyWord<double, std::int64_t>
783 : {
784 1545 : static inline void f(const double dfValueIn, std::int64_t &nValueOut)
785 : {
786 1545 : if (CPLIsNan(dfValueIn))
787 : {
788 1 : nValueOut = 0;
789 : }
790 1544 : else if (dfValueIn >=
791 1544 : static_cast<double>(cpl::NumericLimits<std::int64_t>::max()))
792 : {
793 4 : nValueOut = cpl::NumericLimits<std::int64_t>::max();
794 : }
795 1540 : else if (dfValueIn <=
796 1540 : static_cast<double>(cpl::NumericLimits<std::int64_t>::min()))
797 : {
798 4 : nValueOut = cpl::NumericLimits<std::int64_t>::min();
799 : }
800 : else
801 : {
802 3072 : nValueOut = static_cast<std::int64_t>(
803 1536 : dfValueIn > 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5);
804 : }
805 1545 : }
806 : };
807 :
808 : // Floating point input to a very large signed integer type: nan becomes zero, otherwise round and clamp.
809 : // Rounding for signed integers is different than for the unsigned integers above.
810 : // Avoid infinity while clamping when the maximum integer is too large for the floating-point type.
811 : // Avoid roundoff while clamping.
812 :
813 : template <> struct sGDALCopyWord<GFloat16, int>
814 : {
815 5222 : static inline void f(const GFloat16 hfValueIn, int &nValueOut)
816 : {
817 5222 : if (CPLIsNan(hfValueIn))
818 : {
819 0 : nValueOut = 0;
820 : }
821 5222 : else if (CPLIsInf(hfValueIn))
822 : {
823 0 : nValueOut = hfValueIn > GFloat16(0.0f)
824 0 : ? cpl::NumericLimits<int>::max()
825 0 : : cpl::NumericLimits<int>::lowest();
826 : }
827 : else
828 : {
829 10444 : nValueOut = static_cast<int>(hfValueIn > GFloat16(0.0f)
830 5222 : ? hfValueIn + GFloat16(0.5f)
831 5889 : : hfValueIn - GFloat16(0.5f));
832 : }
833 5222 : }
834 : };
835 :
836 : template <> struct sGDALCopyWord<GFloat16, std::int64_t>
837 : {
838 1093 : static inline void f(const GFloat16 hfValueIn, std::int64_t &nValueOut)
839 : {
840 1093 : if (CPLIsNan(hfValueIn))
841 : {
842 1 : nValueOut = 0;
843 : }
844 1092 : else if (CPLIsInf(hfValueIn))
845 : {
846 4 : nValueOut = hfValueIn > GFloat16(0.0f)
847 2 : ? cpl::NumericLimits<std::int64_t>::max()
848 1 : : cpl::NumericLimits<std::int64_t>::lowest();
849 : }
850 : else
851 : {
852 1090 : nValueOut = static_cast<std::int64_t>(
853 1090 : hfValueIn > GFloat16(0.0f) ? hfValueIn + GFloat16(0.5f)
854 1560 : : hfValueIn - GFloat16(0.5f));
855 : }
856 1093 : }
857 : };
858 :
859 : /**
860 : * Copy a single word, optionally rounding if appropriate (i.e. going
861 : * from the float to the integer case). Note that this is the function
862 : * you should specialize if you're adding a new data type.
863 : *
864 : * @param tValueIn value of type Tin; the input value to be converted
865 : * @param tValueOut value of type Tout; the output value
866 : */
867 :
868 : template <class Tin, class Tout>
869 629549892 : inline void GDALCopyWord(const Tin tValueIn, Tout &tValueOut)
870 : {
871 : if constexpr (std::is_same<Tin, Tout>::value)
872 28541708 : tValueOut = tValueIn;
873 : else
874 601008184 : sGDALCopyWord<Tin, Tout>::f(tValueIn, tValueOut);
875 627887672 : }
876 :
877 : /************************************************************************/
878 : /* GDALCopy4Words() */
879 : /************************************************************************/
880 : /**
881 : * Copy 4 packed words to 4 packed words, optionally rounding if appropriate
882 : * (i.e. going from the float to the integer case).
883 : *
884 : * @param pValueIn pointer to 4 input values of type Tin.
885 : * @param pValueOut pointer to 4 output values of type Tout.
886 : */
887 :
888 : template <class Tin, class Tout>
889 604 : inline void GDALCopy4Words(const Tin *pValueIn, Tout *const pValueOut)
890 : {
891 604 : GDALCopyWord(pValueIn[0], pValueOut[0]);
892 604 : GDALCopyWord(pValueIn[1], pValueOut[1]);
893 604 : GDALCopyWord(pValueIn[2], pValueOut[2]);
894 604 : GDALCopyWord(pValueIn[3], pValueOut[3]);
895 604 : }
896 :
897 : /************************************************************************/
898 : /* GDALCopy8Words() */
899 : /************************************************************************/
900 : /**
901 : * Copy 8 packed words to 8 packed words, optionally rounding if appropriate
902 : * (i.e. going from the float to the integer case).
903 : *
904 : * @param pValueIn pointer to 8 input values of type Tin.
905 : * @param pValueOut pointer to 8 output values of type Tout.
906 : */
907 :
908 : template <class Tin, class Tout>
909 26719967 : inline void GDALCopy8Words(const Tin *pValueIn, Tout *const pValueOut)
910 : {
911 26719967 : GDALCopy4Words(pValueIn, pValueOut);
912 26715740 : GDALCopy4Words(pValueIn + 4, pValueOut + 4);
913 26718994 : }
914 :
915 : // Needs SSE2
916 : #if defined(__x86_64) || defined(_M_X64) || defined(USE_SSE2) || \
917 : defined(USE_NEON_OPTIMIZATIONS)
918 :
919 : template <>
920 34978702 : inline void GDALCopy4Words(const float *pValueIn, GByte *const pValueOut)
921 : {
922 34978702 : __m128 xmm = _mm_loadu_ps(pValueIn);
923 :
924 : // The following clamping would be useless due to the final saturating
925 : // packing if we could guarantee the input range in [INT_MIN,INT_MAX]
926 34978702 : const __m128 p0d5 = _mm_set1_ps(0.5f);
927 34978702 : const __m128 xmm_max = _mm_set1_ps(255);
928 34978702 : xmm = _mm_add_ps(xmm, p0d5);
929 69981904 : xmm = _mm_min_ps(_mm_max_ps(xmm, p0d5), xmm_max);
930 :
931 34993602 : __m128i xmm_i = _mm_cvttps_epi32(xmm);
932 :
933 : #if defined(__SSSE3__) || defined(USE_NEON_OPTIMIZATIONS)
934 : xmm_i = _mm_shuffle_epi8(
935 : xmm_i, _mm_cvtsi32_si128(0 | (4 << 8) | (8 << 16) | (12 << 24)));
936 : #else
937 34992002 : xmm_i = _mm_packs_epi32(xmm_i, xmm_i); // Pack int32 to int16
938 34994002 : xmm_i = _mm_packus_epi16(xmm_i, xmm_i); // Pack int16 to uint8
939 : #endif
940 34994002 : GDALCopyXMMToInt32(xmm_i, pValueOut);
941 34978902 : }
942 :
943 : template <>
944 3363960 : inline void GDALCopy4Words(const float *pValueIn, GInt16 *const pValueOut)
945 : {
946 3363960 : __m128 xmm = _mm_loadu_ps(pValueIn);
947 :
948 3363960 : const __m128 xmm_min = _mm_set1_ps(-32768);
949 3363960 : const __m128 xmm_max = _mm_set1_ps(32767);
950 6727920 : xmm = _mm_min_ps(_mm_max_ps(xmm, xmm_min), xmm_max);
951 :
952 3363960 : const __m128 p0d5 = _mm_set1_ps(0.5f);
953 3363960 : const __m128 m0d5 = _mm_set1_ps(-0.5f);
954 3363960 : const __m128 mask = _mm_cmpge_ps(xmm, p0d5);
955 : // f >= 0.5f ? f + 0.5f : f - 0.5f
956 13455800 : xmm = _mm_add_ps(
957 : xmm, _mm_or_ps(_mm_and_ps(mask, p0d5), _mm_andnot_ps(mask, m0d5)));
958 :
959 3363960 : __m128i xmm_i = _mm_cvttps_epi32(xmm);
960 :
961 3363960 : xmm_i = _mm_packs_epi32(xmm_i, xmm_i); // Pack int32 to int16
962 3363960 : GDALCopyXMMToInt64(xmm_i, pValueOut);
963 3363960 : }
964 :
965 : template <>
966 1 : inline void GDALCopy4Words(const float *pValueIn, GUInt16 *const pValueOut)
967 : {
968 1 : __m128 xmm = _mm_loadu_ps(pValueIn);
969 :
970 1 : const __m128 p0d5 = _mm_set1_ps(0.5f);
971 1 : const __m128 xmm_max = _mm_set1_ps(65535);
972 1 : xmm = _mm_add_ps(xmm, p0d5);
973 2 : xmm = _mm_min_ps(_mm_max_ps(xmm, p0d5), xmm_max);
974 :
975 1 : __m128i xmm_i = _mm_cvttps_epi32(xmm);
976 :
977 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
978 : xmm_i = _mm_packus_epi32(xmm_i, xmm_i); // Pack int32 to uint16
979 : #else
980 : // Translate to int16 range because _mm_packus_epi32 is SSE4.1 only
981 2 : xmm_i = _mm_add_epi32(xmm_i, _mm_set1_epi32(-32768));
982 1 : xmm_i = _mm_packs_epi32(xmm_i, xmm_i); // Pack int32 to int16
983 : // Translate back to uint16 range (actually -32768==32768 in int16)
984 1 : xmm_i = _mm_add_epi16(xmm_i, _mm_set1_epi16(-32768));
985 : #endif
986 1 : GDALCopyXMMToInt64(xmm_i, pValueOut);
987 1 : }
988 :
989 : template <>
990 621041 : inline void GDALCopy4Words(const double *pValueIn, float *const pValueOut)
991 : {
992 621041 : const __m128d val01 = _mm_loadu_pd(pValueIn);
993 1242080 : const __m128d val23 = _mm_loadu_pd(pValueIn + 2);
994 621081 : const __m128 val01_s = _mm_cvtpd_ps(val01);
995 621229 : const __m128 val23_s = _mm_cvtpd_ps(val23);
996 621213 : const __m128 val = _mm_movelh_ps(val01_s, val23_s);
997 : _mm_storeu_ps(pValueOut, val);
998 621213 : }
999 :
1000 : template <>
1001 2767140 : inline void GDALCopy4Words(const double *pValueIn, GByte *const pValueOut)
1002 : {
1003 2767140 : const __m128d p0d5 = _mm_set1_pd(0.5);
1004 2767140 : const __m128d xmm_max = _mm_set1_pd(255);
1005 :
1006 2767140 : __m128d val01 = _mm_loadu_pd(pValueIn);
1007 5534280 : __m128d val23 = _mm_loadu_pd(pValueIn + 2);
1008 2767140 : val01 = _mm_add_pd(val01, p0d5);
1009 5535800 : val01 = _mm_min_pd(_mm_max_pd(val01, p0d5), xmm_max);
1010 2768210 : val23 = _mm_add_pd(val23, p0d5);
1011 5537170 : val23 = _mm_min_pd(_mm_max_pd(val23, p0d5), xmm_max);
1012 :
1013 2768020 : const __m128i val01_u32 = _mm_cvttpd_epi32(val01);
1014 2768460 : const __m128i val23_u32 = _mm_cvttpd_epi32(val23);
1015 :
1016 : // Merge 4 int32 values into a single register
1017 8305030 : auto xmm_i = _mm_castpd_si128(_mm_shuffle_pd(
1018 : _mm_castsi128_pd(val01_u32), _mm_castsi128_pd(val23_u32), 0));
1019 :
1020 : #if defined(__SSSE3__) || defined(USE_NEON_OPTIMIZATIONS)
1021 : xmm_i = _mm_shuffle_epi8(
1022 : xmm_i, _mm_cvtsi32_si128(0 | (4 << 8) | (8 << 16) | (12 << 24)));
1023 : #else
1024 2768040 : xmm_i = _mm_packs_epi32(xmm_i, xmm_i); // Pack int32 to int16
1025 2769380 : xmm_i = _mm_packus_epi16(xmm_i, xmm_i); // Pack int16 to uint8
1026 : #endif
1027 2769380 : GDALCopyXMMToInt32(xmm_i, pValueOut);
1028 2767240 : }
1029 :
1030 : template <>
1031 11688000 : inline void GDALCopy4Words(const float *pValueIn, double *const pValueOut)
1032 : {
1033 11688000 : const __m128 valIn = _mm_loadu_ps(pValueIn);
1034 11688000 : _mm_storeu_pd(pValueOut, _mm_cvtps_pd(valIn));
1035 23376000 : _mm_storeu_pd(pValueOut + 2, _mm_cvtps_pd(_mm_movehl_ps(valIn, valIn)));
1036 11688000 : }
1037 :
1038 : #ifdef __F16C__
1039 : template <>
1040 : inline void GDALCopy4Words(const GFloat16 *pValueIn, float *const pValueOut)
1041 : {
1042 : __m128i xmm = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(pValueIn));
1043 : _mm_storeu_ps(pValueOut, _mm_cvtph_ps(xmm));
1044 : }
1045 :
1046 : template <>
1047 : inline void GDALCopy4Words(const float *pValueIn, GFloat16 *const pValueOut)
1048 : {
1049 : __m128 xmm = _mm_loadu_ps(pValueIn);
1050 : GDALCopyXMMToInt64(_mm_cvtps_ph(xmm, _MM_FROUND_TO_NEAREST_INT), pValueOut);
1051 : }
1052 :
1053 : template <>
1054 : inline void GDALCopy4Words(const GFloat16 *pValueIn, double *const pValueOut)
1055 : {
1056 : float tmp[4];
1057 : GDALCopy4Words(pValueIn, tmp);
1058 : GDALCopy4Words(tmp, pValueOut);
1059 : }
1060 :
1061 : template <>
1062 : inline void GDALCopy4Words(const double *pValueIn, GFloat16 *const pValueOut)
1063 : {
1064 : float tmp[4];
1065 : GDALCopy4Words(pValueIn, tmp);
1066 : GDALCopy4Words(tmp, pValueOut);
1067 : }
1068 : #else // !__F16C__
1069 :
1070 83336 : static inline __m128i GDALIfThenElse(__m128i mask, __m128i thenVal,
1071 : __m128i elseVal)
1072 : {
1073 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
1074 : return _mm_blendv_epi8(elseVal, thenVal, mask);
1075 : #else
1076 166672 : return _mm_or_si128(_mm_and_si128(mask, thenVal),
1077 83336 : _mm_andnot_si128(mask, elseVal));
1078 : #endif
1079 : }
1080 :
1081 : // Convert 4 float16 values to 4 float 32 values
1082 : // xmm must contain 4 float16 values stored in 32 bit each (with upper 16 bits at zero)
1083 41668 : static inline __m128i GDALFourFloat16ToFloat32(__m128i xmm)
1084 : {
1085 : // Ported from https://github.com/simd-everywhere/simde/blob/51743e7920b6e867678cb50e9c62effe28f70b33/simde/simde-f16.h#L242C4-L242C68
1086 : // to SSE2 in a branch-less way
1087 :
1088 : /* This code is CC0, based heavily on code by Fabian Giesen. */
1089 : const auto denorm_magic =
1090 83336 : _mm_castsi128_ps(_mm_set1_epi32((128 - 15) << 23));
1091 : const auto shifted_exp =
1092 41668 : _mm_set1_epi32(0x7c00 << 13); /* exponent mask after shift */
1093 :
1094 : // Shift exponent and mantissa bits to their position in a float32
1095 125004 : auto f32u = _mm_slli_epi32(_mm_and_si128(xmm, _mm_set1_epi32(0x7fff)), 13);
1096 : // Extract the (shifted) exponent
1097 41668 : const auto exp = _mm_and_si128(shifted_exp, f32u);
1098 : // Adjust the exponent
1099 41668 : const auto exp_adjustment = _mm_set1_epi32((127 - 15) << 23);
1100 41668 : f32u = _mm_add_epi32(f32u, exp_adjustment);
1101 :
1102 41668 : const auto is_inf_nan = _mm_cmpeq_epi32(exp, shifted_exp); /* Inf/NaN? */
1103 : // When is_inf_nan is true: extra exponent adjustment
1104 41668 : const auto f32u_inf_nan = _mm_add_epi32(f32u, exp_adjustment);
1105 :
1106 : const auto is_denormal =
1107 83336 : _mm_cmpeq_epi32(exp, _mm_setzero_si128()); /* Zero/Denormal? */
1108 : // When is_denormal is true:
1109 83336 : auto f32u_denormal = _mm_add_epi32(f32u, _mm_set1_epi32(1 << 23));
1110 83336 : f32u_denormal = _mm_castps_si128(
1111 : _mm_sub_ps(_mm_castsi128_ps(f32u_denormal), denorm_magic));
1112 :
1113 41668 : f32u = GDALIfThenElse(is_inf_nan, f32u_inf_nan, f32u);
1114 41668 : f32u = GDALIfThenElse(is_denormal, f32u_denormal, f32u);
1115 :
1116 : // Re-apply sign bit
1117 125004 : f32u = _mm_or_si128(
1118 : f32u, _mm_slli_epi32(_mm_and_si128(xmm, _mm_set1_epi32(0x8000)), 16));
1119 41668 : return f32u;
1120 : }
1121 :
1122 : template <>
1123 484 : inline void GDALCopy8Words(const GFloat16 *pValueIn, float *const pValueOut)
1124 : {
1125 484 : __m128i xmm = _mm_loadu_si128(reinterpret_cast<const __m128i *>(pValueIn));
1126 : const auto xmm_0 =
1127 968 : GDALFourFloat16ToFloat32(_mm_unpacklo_epi16(xmm, _mm_setzero_si128()));
1128 : const auto xmm_1 =
1129 968 : GDALFourFloat16ToFloat32(_mm_unpackhi_epi16(xmm, _mm_setzero_si128()));
1130 484 : _mm_storeu_ps(pValueOut + 0, _mm_castsi128_ps(xmm_0));
1131 484 : _mm_storeu_ps(pValueOut + 4, _mm_castsi128_ps(xmm_1));
1132 484 : }
1133 :
1134 : template <>
1135 20350 : inline void GDALCopy8Words(const GFloat16 *pValueIn, double *const pValueOut)
1136 : {
1137 20350 : __m128i xmm = _mm_loadu_si128(reinterpret_cast<const __m128i *>(pValueIn));
1138 61050 : const auto xmm_0 = _mm_castsi128_ps(
1139 : GDALFourFloat16ToFloat32(_mm_unpacklo_epi16(xmm, _mm_setzero_si128())));
1140 61050 : const auto xmm_1 = _mm_castsi128_ps(
1141 : GDALFourFloat16ToFloat32(_mm_unpackhi_epi16(xmm, _mm_setzero_si128())));
1142 20350 : _mm_storeu_pd(pValueOut + 0, _mm_cvtps_pd(xmm_0));
1143 40700 : _mm_storeu_pd(pValueOut + 2, _mm_cvtps_pd(_mm_movehl_ps(xmm_0, xmm_0)));
1144 20350 : _mm_storeu_pd(pValueOut + 4, _mm_cvtps_pd(xmm_1));
1145 40700 : _mm_storeu_pd(pValueOut + 6, _mm_cvtps_pd(_mm_movehl_ps(xmm_1, xmm_1)));
1146 20350 : }
1147 :
1148 : #endif // __F16C__
1149 :
1150 : #ifdef __AVX2__
1151 :
1152 : #include <immintrin.h>
1153 :
1154 : template <>
1155 : inline void GDALCopy8Words(const double *pValueIn, float *const pValueOut)
1156 : {
1157 : const __m256d val0123 = _mm256_loadu_pd(pValueIn);
1158 : const __m256d val4567 = _mm256_loadu_pd(pValueIn + 4);
1159 : const __m256 val0123_s = _mm256_castps128_ps256(_mm256_cvtpd_ps(val0123));
1160 : const __m256 val4567_s = _mm256_castps128_ps256(_mm256_cvtpd_ps(val4567));
1161 : const __m256 val =
1162 : _mm256_permute2f128_ps(val0123_s, val4567_s, 0 | (2 << 4));
1163 : _mm256_storeu_ps(pValueOut, val);
1164 : }
1165 :
1166 : template <>
1167 : inline void GDALCopy8Words(const float *pValueIn, double *const pValueOut)
1168 : {
1169 : const __m256 valIn = _mm256_loadu_ps(pValueIn);
1170 : _mm256_storeu_pd(pValueOut, _mm256_cvtps_pd(_mm256_castps256_ps128(valIn)));
1171 : _mm256_storeu_pd(pValueOut + 4,
1172 : _mm256_cvtps_pd(_mm256_castps256_ps128(
1173 : _mm256_permute2f128_ps(valIn, valIn, 1))));
1174 : }
1175 :
1176 : #ifdef __F16C__
1177 :
1178 : template <>
1179 : inline void GDALCopy8Words(const GFloat16 *pValueIn, float *const pValueOut)
1180 : {
1181 : __m128i xmm = _mm_loadu_si128(reinterpret_cast<const __m128i *>(pValueIn));
1182 : _mm256_storeu_ps(pValueOut, _mm256_cvtph_ps(xmm));
1183 : }
1184 :
1185 : template <>
1186 : inline void GDALCopy8Words(const float *pValueIn, GFloat16 *const pValueOut)
1187 : {
1188 : __m256 ymm = _mm256_loadu_ps(pValueIn);
1189 : _mm_storeu_si128(reinterpret_cast<__m128i *>(pValueOut),
1190 : _mm256_cvtps_ph(ymm, _MM_FROUND_TO_NEAREST_INT));
1191 : }
1192 :
1193 : template <>
1194 : inline void GDALCopy8Words(const GFloat16 *pValueIn, double *const pValueOut)
1195 : {
1196 : __m128i xmm = _mm_loadu_si128(reinterpret_cast<const __m128i *>(pValueIn));
1197 : const auto ymm = _mm256_cvtph_ps(xmm);
1198 : _mm256_storeu_pd(pValueOut, _mm256_cvtps_pd(_mm256_extractf128_ps(ymm, 0)));
1199 : _mm256_storeu_pd(pValueOut + 4,
1200 : _mm256_cvtps_pd(_mm256_extractf128_ps(ymm, 1)));
1201 : }
1202 :
1203 : template <>
1204 : inline void GDALCopy8Words(const double *pValueIn, GFloat16 *const pValueOut)
1205 : {
1206 : __m256d ymm0 = _mm256_loadu_pd(pValueIn);
1207 : __m256d ymm1 = _mm256_loadu_pd(pValueIn + 4);
1208 : __m256 ymm = _mm256_set_m128(_mm256_cvtpd_ps(ymm1), _mm256_cvtpd_ps(ymm0));
1209 : _mm_storeu_si128(reinterpret_cast<__m128i *>(pValueOut),
1210 : _mm256_cvtps_ph(ymm, _MM_FROUND_TO_NEAREST_INT));
1211 : }
1212 :
1213 : #endif
1214 :
1215 : template <>
1216 : inline void GDALCopy8Words(const float *pValueIn, GByte *const pValueOut)
1217 : {
1218 : __m256 ymm = _mm256_loadu_ps(pValueIn);
1219 :
1220 : const __m256 p0d5 = _mm256_set1_ps(0.5f);
1221 : const __m256 ymm_max = _mm256_set1_ps(255);
1222 : ymm = _mm256_add_ps(ymm, p0d5);
1223 : ymm = _mm256_min_ps(_mm256_max_ps(ymm, p0d5), ymm_max);
1224 :
1225 : __m256i ymm_i = _mm256_cvttps_epi32(ymm);
1226 :
1227 : ymm_i = _mm256_packus_epi32(ymm_i, ymm_i); // Pack int32 to uint16
1228 : ymm_i = _mm256_permute4x64_epi64(ymm_i, 0 | (2 << 2)); // AVX2
1229 :
1230 : __m128i xmm_i = _mm256_castsi256_si128(ymm_i);
1231 : xmm_i = _mm_packus_epi16(xmm_i, xmm_i);
1232 : GDALCopyXMMToInt64(xmm_i, pValueOut);
1233 : }
1234 :
1235 : template <>
1236 : inline void GDALCopy8Words(const float *pValueIn, GUInt16 *const pValueOut)
1237 : {
1238 : __m256 ymm = _mm256_loadu_ps(pValueIn);
1239 :
1240 : const __m256 p0d5 = _mm256_set1_ps(0.5f);
1241 : const __m256 ymm_max = _mm256_set1_ps(65535);
1242 : ymm = _mm256_add_ps(ymm, p0d5);
1243 : ymm = _mm256_min_ps(_mm256_max_ps(ymm, p0d5), ymm_max);
1244 :
1245 : __m256i ymm_i = _mm256_cvttps_epi32(ymm);
1246 :
1247 : ymm_i = _mm256_packus_epi32(ymm_i, ymm_i); // Pack int32 to uint16
1248 : ymm_i = _mm256_permute4x64_epi64(ymm_i, 0 | (2 << 2)); // AVX2
1249 :
1250 : _mm_storeu_si128(reinterpret_cast<__m128i *>(pValueOut),
1251 : _mm256_castsi256_si128(ymm_i));
1252 : }
1253 : #else
1254 : template <>
1255 7773521 : inline void GDALCopy8Words(const float *pValueIn, GUInt16 *const pValueOut)
1256 : {
1257 7773521 : __m128 xmm = _mm_loadu_ps(pValueIn);
1258 15547002 : __m128 xmm1 = _mm_loadu_ps(pValueIn + 4);
1259 :
1260 7773521 : const __m128 p0d5 = _mm_set1_ps(0.5f);
1261 7773521 : const __m128 xmm_max = _mm_set1_ps(65535);
1262 7773521 : xmm = _mm_add_ps(xmm, p0d5);
1263 7773521 : xmm1 = _mm_add_ps(xmm1, p0d5);
1264 15540102 : xmm = _mm_min_ps(_mm_max_ps(xmm, p0d5), xmm_max);
1265 15537302 : xmm1 = _mm_min_ps(_mm_max_ps(xmm1, p0d5), xmm_max);
1266 :
1267 7769721 : __m128i xmm_i = _mm_cvttps_epi32(xmm);
1268 7766511 : __m128i xmm1_i = _mm_cvttps_epi32(xmm1);
1269 :
1270 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
1271 : xmm_i = _mm_packus_epi32(xmm_i, xmm1_i); // Pack int32 to uint16
1272 : #else
1273 : // Translate to int16 range because _mm_packus_epi32 is SSE4.1 only
1274 15533002 : xmm_i = _mm_add_epi32(xmm_i, _mm_set1_epi32(-32768));
1275 15533002 : xmm1_i = _mm_add_epi32(xmm1_i, _mm_set1_epi32(-32768));
1276 7774561 : xmm_i = _mm_packs_epi32(xmm_i, xmm1_i); // Pack int32 to int16
1277 : // Translate back to uint16 range (actually -32768==32768 in int16)
1278 15549102 : xmm_i = _mm_add_epi16(xmm_i, _mm_set1_epi16(-32768));
1279 : #endif
1280 : _mm_storeu_si128(reinterpret_cast<__m128i *>(pValueOut), xmm_i);
1281 7774561 : }
1282 : #endif
1283 :
1284 : #endif // defined(__x86_64) || defined(_M_X64)
1285 :
1286 : #endif // GDAL_PRIV_TEMPLATES_HPP_INCLUDED
|
