LCOV - code coverage report
Current view: top level - gcore - gdal_priv_templates.hpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 389 413 94.2 %
Date: 2025-09-10 17:48:50 Functions: 363 365 99.5 %

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

Generated by: LCOV version 1.14