From ce7a57ae17a7c30f1bd7108ed9e220fe149fd55e Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Tue, 4 Feb 2025 03:52:44 -0500 Subject: [PATCH] Initial implementation of float16 support (RFC100) (#11180) --- .github/workflows/icc/build.sh | 9 +- CMakeLists.txt | 17 + MIGRATION_GUIDE.TXT | 8 +- alg/gdalchecksum.cpp | 3 +- alg/gdalmediancut.cpp | 5 +- alg/gdalpansharpen.cpp | 22 +- alg/gdalrasterize.cpp | 9 + alg/gdalwarper.cpp | 7 +- alg/gdalwarpkernel.cpp | 174 +++-- apps/gdaldem_lib.cpp | 9 +- apps/gdalinfo_lib.cpp | 6 + apps/gdalmdiminfo_lib.cpp | 6 + autotest/cpp/CMakeLists.txt | 2 + autotest/cpp/test_gdal.cpp | 581 +++++++++------- autotest/cpp/test_gdal_minmax_element.cpp | 42 +- autotest/cpp/testcopywords.cpp | 246 ++++++- autotest/cpp/testfloat16.cpp | 187 ++++++ autotest/gcore/numpy_rw.py | 2 +- autotest/gdrivers/vrtderived.py | 40 +- autotest/gdrivers/zarr_driver.py | 67 +- ci/travis/csa_common/install.sh | 2 +- ci/travis/csa_common/script.sh | 2 +- cmake/helpers/configure.cmake | 16 + cmake/template/cpl_config.h.in | 6 + doc/source/development/dev_practices.rst | 1 + doc/source/user/raster_data_model.rst | 6 +- frmts/envisat/records.h | 2 + frmts/gti/data/gdaltileindex.xsd | 2 + frmts/gtiff/gt_overview.cpp | 14 + frmts/gtiff/gtiffdataset_write.cpp | 1 + frmts/hdf5/hdf5multidim.cpp | 8 + frmts/netcdf/netcdfdataset.cpp | 1 + frmts/tiledb/tiledbdense.cpp | 12 +- frmts/tiledb/tiledbmultidimarray.cpp | 5 + frmts/tiledb/tiledbsparse.cpp | 1 + frmts/vrt/data/gdalvrt.xsd | 3 + frmts/vrt/pixelfunctions.cpp | 4 + frmts/vrt/vrtderivedrasterband.cpp | 25 +- frmts/vrt/vrtdriver.cpp | 4 +- frmts/vrt/vrtprocesseddatasetfunctions.cpp | 5 + frmts/vrt/vrtrasterband.cpp | 13 +- frmts/zarr/zarr_array.cpp | 2 +- frmts/zarr/zarr_v2_array.cpp | 13 +- frmts/zarr/zarr_v2_group.cpp | 22 +- frmts/zarr/zarr_v3_array.cpp | 20 +- frmts/zarr/zarr_v3_group.cpp | 16 +- frmts/zarr/zarrdrivercore.cpp | 7 +- fuzzers/build.sh | 4 +- gcore/gdal.h | 4 +- gcore/gdal_misc.cpp | 375 +++++++---- gcore/gdal_priv.h | 36 +- gcore/gdal_priv_templates.hpp | 533 +++++++++++---- gcore/gdal_typetraits.h | 46 +- gcore/gdalcachedpixelaccessor.h | 5 + gcore/gdalmultidim.cpp | 45 +- gcore/gdalmultidim_rat.cpp | 1 + gcore/gdalnodatamaskband.cpp | 8 + gcore/gdalnodatavaluesmaskband.cpp | 2 + gcore/gdalrasterband.cpp | 472 ++++++++----- gcore/overview.cpp | 9 +- gcore/rasterio.cpp | 39 +- gcore/rawdataset.cpp | 5 + .../gpkg/gdalgeopackagerasterband.cpp | 1 + port/CMakeLists.txt | 1 + port/cpl_float.cpp | 18 + port/cpl_float.h | 635 ++++++++++++++++++ port/cpl_json_streaming_writer.cpp | 20 + port/cpl_json_streaming_writer.h | 5 +- port/cpl_port.h | 91 +-- scripts/cppcheck.sh | 4 +- swig/include/gdal.i | 4 +- swig/include/gdal_array.i | 19 +- swig/include/gdalconst.i | 2 + swig/include/java/gdal_java.i | 24 +- swig/include/python/gdal_python.i | 16 +- swig/include/python/typemaps_python.i | 5 + .../gdal-utils/osgeo_utils/gdal_calc.py | 7 + .../gdal-utils/osgeo_utils/samples/fft.py | 4 + .../gdal-utils/osgeo_utils/samples/rel.py | 4 + 79 files changed, 3132 insertions(+), 967 deletions(-) create mode 100644 autotest/cpp/testfloat16.cpp diff --git a/.github/workflows/icc/build.sh b/.github/workflows/icc/build.sh index f353e49c479a..52ccde170ad6 100755 --- a/.github/workflows/icc/build.sh +++ b/.github/workflows/icc/build.sh @@ -5,11 +5,16 @@ set -eu # for precompiled headers ccache --set-config sloppiness=pch_defines,time_macros,include_file_mtime,include_file_ctime +# Set C and C++ compiler flags to disable `_Float16`. This is +# necessary because the system C and C++ compilers don't support it, +# and Python's `build_ext` will use the system compiler to build GDAL +# Python extensions. cmake ${GDAL_SOURCE_DIR:=..} \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_C_COMPILER=icx \ -DCMAKE_CXX_COMPILER=icx \ - "-DUSE_PRECOMPILED_HEADERS=ON" \ + -DCMAKE_C_FLAGS=-DGDAL_DISABLE_FLOAT16 \ + -DCMAKE_CXX_FLAGS=-DGDAL_DISABLE_FLOAT16 \ + -DUSE_PRECOMPILED_HEADERS=ON \ -DUSE_CCACHE=ON make -j$(nproc) - diff --git a/CMakeLists.txt b/CMakeLists.txt index f0c98ca89db8..d8ea752ca587 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -111,6 +111,23 @@ else() cmake_pop_check_state() endif () + +# Check whether std::float16_t is available and is working +include(CheckCXXSourceCompiles) +check_cxx_source_compiles( + " + #include + #include + int main() { + std::float16_t x = 0; + using std::nextafter; + std::float16_t y = nextafter(x, x); + return y == 0 ? 0 : 1; + } + " + HAVE_STD_FLOAT16_T +) + # option(CLANG_TIDY_ENABLED "Run clang-tidy with the compiler." OFF) set(CLANG_TIDY_CHECKS diff --git a/MIGRATION_GUIDE.TXT b/MIGRATION_GUIDE.TXT index 509ff1c29d86..c4ad53a63d5e 100644 --- a/MIGRATION_GUIDE.TXT +++ b/MIGRATION_GUIDE.TXT @@ -1,5 +1,11 @@ MIGRATION GUIDE FROM GDAL 3.10 to GDAL 3.11 ------------------------------------------- +------------------------------------------- + +- GDAL drivers may now return raster bands with the new data types + GDT_Float16 or GDT_CFloat16. Code that use the GDAL API must be + ready to react to the new data type, possibly by doing RasterIO() + requests with eBufType==GDT_Float32, if they can't deal natively + with Float16 values. - If only a specific GDAL Minor version is to be supported, this must now be specified in the find_package call in CMake via a version range specification. diff --git a/alg/gdalchecksum.cpp b/alg/gdalchecksum.cpp index 39ee2ba21314..517b7a010419 100644 --- a/alg/gdalchecksum.cpp +++ b/alg/gdalchecksum.cpp @@ -58,7 +58,8 @@ int CPL_STDCALL GDALChecksumImage(GDALRasterBandH hBand, int nXOff, int nYOff, const GDALDataType eDataType = GDALGetRasterDataType(hBand); const bool bComplex = CPL_TO_BOOL(GDALDataTypeIsComplex(eDataType)); const bool bIsFloatingPoint = - (eDataType == GDT_Float32 || eDataType == GDT_Float64 || + (eDataType == GDT_Float16 || eDataType == GDT_Float32 || + eDataType == GDT_Float64 || eDataType == GDT_CFloat16 || eDataType == GDT_CFloat32 || eDataType == GDT_CFloat64); const auto IntFromDouble = [](double dfVal) diff --git a/alg/gdalmediancut.cpp b/alg/gdalmediancut.cpp index 4ffa7931cfb7..cb9415e92473 100644 --- a/alg/gdalmediancut.cpp +++ b/alg/gdalmediancut.cpp @@ -32,6 +32,7 @@ #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_progress.h" #include "cpl_vsi.h" #include "gdal.h" @@ -328,7 +329,7 @@ int GDALComputeMedianCutPCTInternal( /* STEP 1: create empty boxes. */ /* ==================================================================== */ if (static_cast(nXSize) > - std::numeric_limits::max() / static_cast(nYSize)) + cpl::NumericLimits::max() / static_cast(nYSize)) { CPLError(CE_Warning, CPLE_AppDefined, "GDALComputeMedianCutPCTInternal() not called " @@ -339,7 +340,7 @@ int GDALComputeMedianCutPCTInternal( if (nBits == 8 && pabyRedBand != nullptr && pabyGreenBand != nullptr && pabyBlueBand != nullptr && static_cast(nXSize) <= - std::numeric_limits::max() / static_cast(nYSize)) + cpl::NumericLimits::max() / static_cast(nYSize)) { nPixels = static_cast(nXSize) * static_cast(nYSize); } diff --git a/alg/gdalpansharpen.cpp b/alg/gdalpansharpen.cpp index 1bce71b58211..77e8665284c5 100644 --- a/alg/gdalpansharpen.cpp +++ b/alg/gdalpansharpen.cpp @@ -26,6 +26,7 @@ #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_multiproc.h" #include "cpl_vsi.h" #include "../frmts/mem/memdataset.h" @@ -800,7 +801,7 @@ void GDALPansharpenOperation::WeightedBroveyPositiveWeights( } if (nMaxValue == 0) - nMaxValue = std::numeric_limits::max(); + nMaxValue = cpl::NumericLimits::max(); size_t j; if (psOptions->nInputSpectralBands == 3 && psOptions->nOutPansharpenedBands == 3 && @@ -1009,6 +1010,12 @@ CPLErr GDALPansharpenOperation::WeightedBrovey( nBandValues, nMaxValue); break; + case GDT_Float16: + WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer, + static_cast(pDataBuf), nValues, + nBandValues, nMaxValue); + break; + case GDT_Float32: WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer, static_cast(pDataBuf), nValues, nBandValues, @@ -1090,6 +1097,12 @@ CPLErr GDALPansharpenOperation::WeightedBrovey( static_cast(pDataBuf), nValues, nBandValues, 0); break; + case GDT_Float16: + WeightedBrovey3( + pPanBuffer, pUpsampledSpectralBuffer, + static_cast(pDataBuf), nValues, nBandValues, 0); + break; + case GDT_Float32: WeightedBrovey3( pPanBuffer, pUpsampledSpectralBuffer, @@ -1769,6 +1782,13 @@ CPLErr GDALPansharpenOperation::PansharpenChunk( pDataBuf, eBufDataType, nValues, nBandValues); break; + case GDT_Float16: + eErr = WeightedBrovey( + static_cast(pPanBuffer), + static_cast(pUpsampledSpectralBuffer), + pDataBuf, eBufDataType, nValues, nBandValues); + break; + case GDT_Float32: eErr = WeightedBrovey( static_cast(pPanBuffer), diff --git a/alg/gdalrasterize.cpp b/alg/gdalrasterize.cpp index 46f6470413a9..46d7a3c0646a 100644 --- a/alg/gdalrasterize.cpp +++ b/alg/gdalrasterize.cpp @@ -244,6 +244,10 @@ static void gvBurnScanline(void *pCBData, int nY, int nXStart, int nXEnd, gvBurnScanlineBasic(psInfo, nY, nXStart, nXEnd, dfVariant); break; + case GDT_Float16: + gvBurnScanlineBasic(psInfo, nY, nXStart, nXEnd, + dfVariant); + break; case GDT_Float32: gvBurnScanlineBasic(psInfo, nY, nXStart, nXEnd, dfVariant); break; @@ -252,6 +256,7 @@ static void gvBurnScanline(void *pCBData, int nY, int nXStart, int nXEnd, break; case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_Unknown: @@ -371,6 +376,9 @@ static void gvBurnPoint(void *pCBData, int nY, int nX, double dfVariant) case GDT_UInt64: gvBurnPointBasic(psInfo, nY, nX, dfVariant); break; + case GDT_Float16: + gvBurnPointBasic(psInfo, nY, nX, dfVariant); + break; case GDT_Float32: gvBurnPointBasic(psInfo, nY, nX, dfVariant); break; @@ -379,6 +387,7 @@ static void gvBurnPoint(void *pCBData, int nY, int nX, double dfVariant) break; case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_Unknown: diff --git a/alg/gdalwarper.cpp b/alg/gdalwarper.cpp index 8cddbca37b20..51ed9898b43c 100644 --- a/alg/gdalwarper.cpp +++ b/alg/gdalwarper.cpp @@ -23,6 +23,7 @@ #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_mask.h" #include "cpl_minixml.h" #include "cpl_progress.h" @@ -310,8 +311,8 @@ static CPLErr GDALWarpNoDataMaskerT(const double *padfNoData, size_t nPixels, int *pbOutAllValid) { // Nothing to do if value is out of range. - if (padfNoData[0] < std::numeric_limits::min() || - padfNoData[0] > std::numeric_limits::max() + 0.000001 || + if (padfNoData[0] < cpl::NumericLimits::min() || + padfNoData[0] > cpl::NumericLimits::max() + 0.000001 || padfNoData[1] != 0.0) { *pbOutAllValid = TRUE; @@ -2099,7 +2100,7 @@ GDALWarpOptions *CPL_STDCALL GDALDeserializeWarpOptions(CPLXMLNode *psTree) CPLString().Printf( "%.16g", -std::numeric_limits::max()) == pszValueIn) { - return -std::numeric_limits::max(); + return std::numeric_limits::lowest(); } else if (eDataType == GDT_Float32 && CPLString().Printf("%.16g", diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp index 527c15580e51..7b8a6a096363 100644 --- a/alg/gdalwarpkernel.cpp +++ b/alg/gdalwarpkernel.cpp @@ -32,6 +32,7 @@ #include "cpl_atomic_ops.h" #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_mask.h" #include "cpl_multiproc.h" #include "cpl_progress.h" @@ -1482,7 +1483,7 @@ template struct sGWKRoundValueT /* unsigned */ template static T GWKRoundValueT(double dfValue) { - return sGWKRoundValueT::is_signed>::eval(dfValue); + return sGWKRoundValueT::is_signed>::eval(dfValue); } template <> float GWKRoundValueT(double dfValue) @@ -1503,10 +1504,10 @@ template <> double GWKRoundValueT(double dfValue) template static CPL_INLINE T GWKClampValueT(double dfValue) { - if (dfValue < std::numeric_limits::min()) - return std::numeric_limits::min(); - else if (dfValue > std::numeric_limits::max()) - return std::numeric_limits::max(); + if (dfValue < cpl::NumericLimits::min()) + return cpl::NumericLimits::min(); + else if (dfValue > cpl::NumericLimits::max()) + return cpl::NumericLimits::max(); else return GWKRoundValueT(dfValue); } @@ -1537,28 +1538,30 @@ inline void AvoidNoData(const GDALWarpKernel *poWK, int iBand, if (poWK->padfDstNoDataReal != nullptr && poWK->padfDstNoDataReal[iBand] == static_cast(pDst[iDstOffset])) { - if constexpr (std::numeric_limits::is_integer) + if constexpr (cpl::NumericLimits::is_integer) { if (pDst[iDstOffset] == - static_cast(std::numeric_limits::lowest())) + static_cast(cpl::NumericLimits::lowest())) { pDst[iDstOffset] = - static_cast(std::numeric_limits::lowest() + 1); + static_cast(cpl::NumericLimits::lowest() + 1); } else pDst[iDstOffset]--; } else { - if (pDst[iDstOffset] == std::numeric_limits::max()) + if (pDst[iDstOffset] == cpl::NumericLimits::max()) { + using std::nextafter; pDst[iDstOffset] = - std::nextafter(pDst[iDstOffset], static_cast(0)); + nextafter(pDst[iDstOffset], static_cast(0)); } else { - pDst[iDstOffset] = std::nextafter( - pDst[iDstOffset], std::numeric_limits::max()); + using std::nextafter; + pDst[iDstOffset] = + nextafter(pDst[iDstOffset], cpl::NumericLimits::max()); } } @@ -1653,16 +1656,17 @@ inline void ClampRoundAndAvoidNoData(const GDALWarpKernel *poWK, int iBand, GByte *pabyDst = poWK->papabyDstImage[iBand]; T *pDst = reinterpret_cast(pabyDst); - if constexpr (std::numeric_limits::is_integer) + if constexpr (cpl::NumericLimits::is_integer) { - if (dfReal < static_cast(std::numeric_limits::lowest())) - pDst[iDstOffset] = static_cast(std::numeric_limits::lowest()); - else if (dfReal > static_cast(std::numeric_limits::max())) - pDst[iDstOffset] = static_cast(std::numeric_limits::max()); + using std::floor; + if (dfReal < static_cast(cpl::NumericLimits::lowest())) + pDst[iDstOffset] = static_cast(cpl::NumericLimits::lowest()); + else if (dfReal > static_cast(cpl::NumericLimits::max())) + pDst[iDstOffset] = static_cast(cpl::NumericLimits::max()); + else if constexpr (cpl::NumericLimits::is_signed) + pDst[iDstOffset] = static_cast(floor(dfReal + 0.5)); else - pDst[iDstOffset] = (std::numeric_limits::is_signed) - ? static_cast(floor(dfReal + 0.5)) - : static_cast(dfReal + 0.5); + pDst[iDstOffset] = static_cast(dfReal + 0.5); } else { @@ -1753,6 +1757,11 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand, dfDstImag = 0.0; break; + case GDT_Float16: + dfDstReal = reinterpret_cast(pabyDst)[iDstOffset]; + dfDstImag = 0.0; + break; + case GDT_Float32: dfDstReal = reinterpret_cast(pabyDst)[iDstOffset]; dfDstImag = 0.0; @@ -1775,6 +1784,13 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand, reinterpret_cast(pabyDst)[iDstOffset * 2 + 1]; break; + case GDT_CFloat16: + dfDstReal = + reinterpret_cast(pabyDst)[iDstOffset * 2]; + dfDstImag = + reinterpret_cast(pabyDst)[iDstOffset * 2 + 1]; + break; + case GDT_CFloat32: dfDstReal = reinterpret_cast(pabyDst)[iDstOffset * 2]; dfDstImag = @@ -1847,6 +1863,10 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand, dfReal); break; + case GDT_Float16: + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); + break; + case GDT_Float32: ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; @@ -1858,23 +1878,21 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand, case GDT_CInt16: { typedef GInt16 T; - if (dfReal < static_cast(std::numeric_limits::min())) + if (dfReal < static_cast(cpl::NumericLimits::min())) reinterpret_cast(pabyDst)[iDstOffset * 2] = - std::numeric_limits::min(); - else if (dfReal > - static_cast(std::numeric_limits::max())) + cpl::NumericLimits::min(); + else if (dfReal > static_cast(cpl::NumericLimits::max())) reinterpret_cast(pabyDst)[iDstOffset * 2] = - std::numeric_limits::max(); + cpl::NumericLimits::max(); else reinterpret_cast(pabyDst)[iDstOffset * 2] = static_cast(floor(dfReal + 0.5)); - if (dfImag < static_cast(std::numeric_limits::min())) + if (dfImag < static_cast(cpl::NumericLimits::min())) reinterpret_cast(pabyDst)[iDstOffset * 2 + 1] = - std::numeric_limits::min(); - else if (dfImag > - static_cast(std::numeric_limits::max())) + cpl::NumericLimits::min(); + else if (dfImag > static_cast(cpl::NumericLimits::max())) reinterpret_cast(pabyDst)[iDstOffset * 2 + 1] = - std::numeric_limits::max(); + cpl::NumericLimits::max(); else reinterpret_cast(pabyDst)[iDstOffset * 2 + 1] = static_cast(floor(dfImag + 0.5)); @@ -1884,29 +1902,34 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand, case GDT_CInt32: { typedef GInt32 T; - if (dfReal < static_cast(std::numeric_limits::min())) + if (dfReal < static_cast(cpl::NumericLimits::min())) reinterpret_cast(pabyDst)[iDstOffset * 2] = - std::numeric_limits::min(); - else if (dfReal > - static_cast(std::numeric_limits::max())) + cpl::NumericLimits::min(); + else if (dfReal > static_cast(cpl::NumericLimits::max())) reinterpret_cast(pabyDst)[iDstOffset * 2] = - std::numeric_limits::max(); + cpl::NumericLimits::max(); else reinterpret_cast(pabyDst)[iDstOffset * 2] = static_cast(floor(dfReal + 0.5)); - if (dfImag < static_cast(std::numeric_limits::min())) + if (dfImag < static_cast(cpl::NumericLimits::min())) reinterpret_cast(pabyDst)[iDstOffset * 2 + 1] = - std::numeric_limits::min(); - else if (dfImag > - static_cast(std::numeric_limits::max())) + cpl::NumericLimits::min(); + else if (dfImag > static_cast(cpl::NumericLimits::max())) reinterpret_cast(pabyDst)[iDstOffset * 2 + 1] = - std::numeric_limits::max(); + cpl::NumericLimits::max(); else reinterpret_cast(pabyDst)[iDstOffset * 2 + 1] = static_cast(floor(dfImag + 0.5)); break; } + case GDT_CFloat16: + reinterpret_cast(pabyDst)[iDstOffset * 2] = + static_cast(dfReal); + reinterpret_cast(pabyDst)[iDstOffset * 2 + 1] = + static_cast(dfImag); + break; + case GDT_CFloat32: reinterpret_cast(pabyDst)[iDstOffset * 2] = static_cast(dfReal); @@ -1999,6 +2022,10 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand, reinterpret_cast(pabyDst)[iDstOffset]); break; + case GDT_Float16: + dfDstReal = reinterpret_cast(pabyDst)[iDstOffset]; + break; + case GDT_Float32: dfDstReal = reinterpret_cast(pabyDst)[iDstOffset]; break; @@ -2009,6 +2036,7 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand, case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_Unknown: @@ -2068,6 +2096,10 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand, dfReal); break; + case GDT_Float16: + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); + break; + case GDT_Float32: ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; @@ -2078,6 +2110,7 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand, case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: return false; @@ -2160,6 +2193,11 @@ static bool GWKGetPixelValue(const GDALWarpKernel *poWK, int iBand, *pdfImag = 0.0; break; + case GDT_Float16: + *pdfReal = reinterpret_cast(pabySrc)[iSrcOffset]; + *pdfImag = 0.0; + break; + case GDT_Float32: *pdfReal = reinterpret_cast(pabySrc)[iSrcOffset]; *pdfImag = 0.0; @@ -2180,6 +2218,12 @@ static bool GWKGetPixelValue(const GDALWarpKernel *poWK, int iBand, *pdfImag = reinterpret_cast(pabySrc)[iSrcOffset * 2 + 1]; break; + case GDT_CFloat16: + *pdfReal = reinterpret_cast(pabySrc)[iSrcOffset * 2]; + *pdfImag = + reinterpret_cast(pabySrc)[iSrcOffset * 2 + 1]; + break; + case GDT_CFloat32: *pdfReal = reinterpret_cast(pabySrc)[iSrcOffset * 2]; *pdfImag = reinterpret_cast(pabySrc)[iSrcOffset * 2 + 1]; @@ -2260,6 +2304,10 @@ static bool GWKGetPixelValueReal(const GDALWarpKernel *poWK, int iBand, reinterpret_cast(pabySrc)[iSrcOffset]); break; + case GDT_Float16: + *pdfReal = reinterpret_cast(pabySrc)[iSrcOffset]; + break; + case GDT_Float32: *pdfReal = reinterpret_cast(pabySrc)[iSrcOffset]; break; @@ -2270,6 +2318,7 @@ static bool GWKGetPixelValueReal(const GDALWarpKernel *poWK, int iBand, case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_Unknown: @@ -2466,6 +2515,19 @@ static bool GWKGetPixelRow(const GDALWarpKernel *poWK, int iBand, break; } + case GDT_Float16: + { + GFloat16 *pSrc = + reinterpret_cast(poWK->papabySrcImage[iBand]); + pSrc += iSrcOffset; + for (int i = 0; i < nSrcLen; i += 2) + { + adfReal[i] = pSrc[i]; + adfReal[i + 1] = pSrc[i + 1]; + } + break; + } + case GDT_Float32: { float *pSrc = @@ -2524,6 +2586,22 @@ static bool GWKGetPixelRow(const GDALWarpKernel *poWK, int iBand, break; } + case GDT_CFloat16: + { + GFloat16 *pSrc = + reinterpret_cast(poWK->papabySrcImage[iBand]); + pSrc += 2 * iSrcOffset; + for (int i = 0; i < nSrcLen; i += 2) + { + adfReal[i] = pSrc[2 * i]; + padfImag[i] = pSrc[2 * i + 1]; + + adfReal[i + 1] = pSrc[2 * i + 2]; + padfImag[i + 1] = pSrc[2 * i + 3]; + } + break; + } + case GDT_CFloat32: { float *pSrc = @@ -4439,7 +4517,7 @@ static void GWKComputeWeights(GDALResampleAlg eResample, int iMin, int iMax, int iC = 0; // Used after for. // Not zero, but as close as possible to it, to avoid potential division by // zero at end of function - double dfAccumulatorWeightHorizontal = std::numeric_limits::min(); + double dfAccumulatorWeightHorizontal = cpl::NumericLimits::min(); for (; i + 2 < iMax; i += 4, iC += 4) { padfWeightsHorizontal[iC] = (i - dfDeltaX) * dfXScale; @@ -4462,7 +4540,7 @@ static void GWKComputeWeights(GDALResampleAlg eResample, int iMin, int iMax, int jC = 0; // Used after for. // Not zero, but as close as possible to it, to avoid potential division by // zero at end of function - double dfAccumulatorWeightVertical = std::numeric_limits::min(); + double dfAccumulatorWeightVertical = cpl::NumericLimits::min(); for (; j + 2 < jMax; j += 4, jC += 4) { padfWeightsVertical[jC] = (j - dfDeltaY) * dfYScale; @@ -7816,7 +7894,7 @@ static void GWKAverageOrModeThread(void *pData) // poWK->eResample == GRA_Max. { bool bFoundValid = false; - double dfTotalReal = std::numeric_limits::lowest(); + double dfTotalReal = cpl::NumericLimits::lowest(); // This code adapted from nAlgo 1 method, GRA_Average. for (int iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++) { @@ -7876,7 +7954,7 @@ static void GWKAverageOrModeThread(void *pData) // poWK->eResample == GRA_Min. { bool bFoundValid = false; - double dfTotalReal = std::numeric_limits::max(); + double dfTotalReal = cpl::NumericLimits::max(); // This code adapted from nAlgo 1 method, GRA_Average. for (int iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++) { @@ -8225,8 +8303,8 @@ static void getConvexPolyIntersection(const XYPoly &poly1, const XYPoly &poly2, return; // Find lowest-left point in intersection set - double lowest_x = std::numeric_limits::max(); - double lowest_y = std::numeric_limits::max(); + double lowest_x = cpl::NumericLimits::max(); + double lowest_y = cpl::NumericLimits::max(); for (const auto &pair : intersection) { const double x = pair.first; @@ -8265,13 +8343,13 @@ static void getConvexPolyIntersection(const XYPoly &poly1, const XYPoly &poly2, double tan_p1; if (p1x_diff == 0.0) - tan_p1 = p1y_diff == 0.0 ? 0.0 : std::numeric_limits::max(); + tan_p1 = p1y_diff == 0.0 ? 0.0 : cpl::NumericLimits::max(); else tan_p1 = p1y_diff / p1x_diff; double tan_p2; if (p2x_diff == 0.0) - tan_p2 = p2y_diff == 0.0 ? 0.0 : std::numeric_limits::max(); + tan_p2 = p2y_diff == 0.0 ? 0.0 : cpl::NumericLimits::max(); else tan_p2 = p2y_diff / p2x_diff; diff --git a/apps/gdaldem_lib.cpp b/apps/gdaldem_lib.cpp index 6fec3cb4e356..9378306f2b5c 100644 --- a/apps/gdaldem_lib.cpp +++ b/apps/gdaldem_lib.cpp @@ -99,6 +99,7 @@ #include #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_progress.h" #include "cpl_string.h" #include "cpl_vsi.h" @@ -329,7 +330,7 @@ static CPLErr GDALGeneric3x3Processing( int bIsSrcNoDataNan = FALSE; T fSrcNoDataValue = 0; - if (std::numeric_limits::is_integer) + if (cpl::NumericLimits::is_integer) { eReadDT = GDT_Int32; if (bSrcHasNoData) @@ -398,7 +399,7 @@ static CPLErr GDALGeneric3x3Processing( return CE_Failure; } - if (std::numeric_limits::is_integer && bSrcHasNoData) + if (cpl::NumericLimits::is_integer && bSrcHasNoData) { abLineHasNoDataValue[i] = false; for (int iX = 0; iX < nXSize; iX++) @@ -484,7 +485,7 @@ static CPLErr GDALGeneric3x3Processing( // In case none of the 3 lines have nodata values, then no need to // check it in ComputeVal() bool bOneOfThreeLinesHasNoData = CPL_TO_BOOL(bSrcHasNoData); - if (std::numeric_limits::is_integer && bSrcHasNoData) + if (cpl::NumericLimits::is_integer && bSrcHasNoData) { bool bLastLineHasNoDataValue = false; int iX = 0; @@ -2654,7 +2655,7 @@ GDALGeneric3x3RasterBand::GDALGeneric3x3RasterBand( const double dfNoDataValue = GDALGetRasterNoDataValue(poDSIn->hSrcBand, &bSrcHasNoData); - if (std::numeric_limits::is_integer) + if (cpl::NumericLimits::is_integer) { eReadDT = GDT_Int32; if (bSrcHasNoData) diff --git a/apps/gdalinfo_lib.cpp b/apps/gdalinfo_lib.cpp index 502a748a7281..3fe5e60fc4f8 100644 --- a/apps/gdalinfo_lib.cpp +++ b/apps/gdalinfo_lib.cpp @@ -1071,6 +1071,9 @@ char *GDALInfo(GDALDatasetH hDataset, const GDALInfoOptions *psOptions) case GDT_Int64: stacDataType = "int64"; break; + case GDT_Float16: + stacDataType = "float16"; + break; case GDT_Float32: stacDataType = "float32"; break; @@ -1083,6 +1086,9 @@ char *GDALInfo(GDALDatasetH hDataset, const GDALInfoOptions *psOptions) case GDT_CInt32: stacDataType = "cint32"; break; + case GDT_CFloat16: + stacDataType = "cfloat16"; + break; case GDT_CFloat32: stacDataType = "cfloat32"; break; diff --git a/apps/gdalmdiminfo_lib.cpp b/apps/gdalmdiminfo_lib.cpp index cf59e16d4dd4..903bf9f48a82 100644 --- a/apps/gdalmdiminfo_lib.cpp +++ b/apps/gdalmdiminfo_lib.cpp @@ -163,6 +163,9 @@ static void DumpValue(CPLJSonStreamingWriter &serializer, const GByte *bytes, case GDT_UInt64: DumpValue(serializer, bytes); break; + case GDT_Float16: + DumpValue(serializer, bytes); + break; case GDT_Float32: DumpValue(serializer, bytes); break; @@ -175,6 +178,9 @@ static void DumpValue(CPLJSonStreamingWriter &serializer, const GByte *bytes, case GDT_CInt32: DumpComplexValue(serializer, bytes); break; + case GDT_CFloat16: + DumpComplexValue(serializer, bytes); + break; case GDT_CFloat32: DumpComplexValue(serializer, bytes); break; diff --git a/autotest/cpp/CMakeLists.txt b/autotest/cpp/CMakeLists.txt index 67a8feea8cc5..9c630e16c2d4 100644 --- a/autotest/cpp/CMakeLists.txt +++ b/autotest/cpp/CMakeLists.txt @@ -299,6 +299,7 @@ set(QUICKTEST_LIST test-block-cache-4 test-block-cache-5 test-block-cache-6 + test-float16 test-copy-words test-closed-on-destroy-DM test-threaded-condition @@ -312,6 +313,7 @@ set(QUICKTEST_LIST test-deferred-plugin ) +gdal_gtest_target(testfloat16 test-float16 testfloat16.cpp) gdal_gtest_target(testcopywords test-copy-words testcopywords.cpp) gdal_gtest_target(testclosedondestroydm test-closed-on-destroy-DM testclosedondestroydm.cpp) gdal_gtest_target(testthreadcond test-threaded-condition testthreadcond.cpp) diff --git a/autotest/cpp/test_gdal.cpp b/autotest/cpp/test_gdal.cpp index e9fde5c5a900..957b7cdb2ff0 100644 --- a/autotest/cpp/test_gdal.cpp +++ b/autotest/cpp/test_gdal.cpp @@ -159,13 +159,17 @@ TEST_F(test_gdal, GDALDataTypeUnion_special_cases) EXPECT_EQ(GDALDataTypeUnion(GDT_Int16, GDT_UInt32), GDT_Int64); EXPECT_EQ(GDALDataTypeUnion(GDT_UInt32, GDT_Int16), GDT_Int64); EXPECT_EQ(GDALDataTypeUnion(GDT_Int64, GDT_UInt64), GDT_Float64); + EXPECT_EQ(GDALDataTypeUnion(GDT_Int64, GDT_Float16), GDT_Float64); EXPECT_EQ(GDALDataTypeUnion(GDT_Int64, GDT_Float32), GDT_Float64); EXPECT_EQ(GDALDataTypeUnion(GDT_Int64, GDT_Float64), GDT_Float64); + EXPECT_EQ(GDALDataTypeUnion(GDT_UInt64, GDT_Float16), GDT_Float64); EXPECT_EQ(GDALDataTypeUnion(GDT_UInt64, GDT_Float32), GDT_Float64); EXPECT_EQ(GDALDataTypeUnion(GDT_UInt64, GDT_Float64), GDT_Float64); EXPECT_EQ(GDALDataTypeUnion(GDT_UInt32, GDT_CInt16), GDT_CFloat64); + EXPECT_EQ(GDALDataTypeUnion(GDT_Float16, GDT_CInt32), GDT_CFloat64); EXPECT_EQ(GDALDataTypeUnion(GDT_Float32, GDT_CInt32), GDT_CFloat64); EXPECT_EQ(GDALDataTypeUnion(GDT_CInt16, GDT_UInt32), GDT_CFloat64); + EXPECT_EQ(GDALDataTypeUnion(GDT_CInt16, GDT_CFloat16), GDT_CFloat32); EXPECT_EQ(GDALDataTypeUnion(GDT_CInt16, GDT_CFloat32), GDT_CFloat32); EXPECT_EQ(GDALDataTypeUnion(GDT_CInt32, GDT_Byte), GDT_CInt32); EXPECT_EQ(GDALDataTypeUnion(GDT_CInt32, GDT_UInt16), GDT_CInt32); @@ -175,6 +179,14 @@ TEST_F(test_gdal, GDALDataTypeUnion_special_cases) EXPECT_EQ(GDALDataTypeUnion(GDT_CInt32, GDT_Float32), GDT_CFloat64); EXPECT_EQ(GDALDataTypeUnion(GDT_CInt32, GDT_CInt16), GDT_CInt32); EXPECT_EQ(GDALDataTypeUnion(GDT_CInt32, GDT_CFloat32), GDT_CFloat64); + EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat16, GDT_Byte), GDT_CFloat32); + EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat16, GDT_UInt16), GDT_CFloat32); + EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat16, GDT_Int16), GDT_CFloat32); + EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat16, GDT_UInt32), GDT_CFloat64); + EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat16, GDT_Int32), GDT_CFloat64); + EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat16, GDT_Float32), GDT_CFloat32); + EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat16, GDT_CInt16), GDT_CFloat32); + EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat16, GDT_CInt32), GDT_CFloat64); EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat32, GDT_Byte), GDT_CFloat32); EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat32, GDT_UInt16), GDT_CFloat32); EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat32, GDT_Int16), GDT_CFloat32); @@ -184,65 +196,41 @@ TEST_F(test_gdal, GDALDataTypeUnion_special_cases) EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat32, GDT_CInt16), GDT_CFloat32); EXPECT_EQ(GDALDataTypeUnion(GDT_CFloat32, GDT_CInt32), GDT_CFloat64); - EXPECT_EQ(GDALFindDataType(0, false /* signed */, false /* floating */, - false /* complex */), - GDT_Byte); - EXPECT_EQ(GDALFindDataType(0, true /* signed */, false /* floating */, - false /* complex */), - GDT_Int8); - EXPECT_EQ(GDALFindDataType(0, false /* signed */, false /* floating */, - true /* complex */), - GDT_CInt32); - EXPECT_EQ(GDALFindDataType(0, true /* signed */, false /* floating */, - true /* complex */), - GDT_CInt16); - EXPECT_EQ(GDALFindDataType(0, false /* signed */, true /* floating */, - false /* complex */), - GDT_Float32); - EXPECT_EQ(GDALFindDataType(0, true /* signed */, true /* floating */, - false /* complex */), - GDT_Float32); - EXPECT_EQ(GDALFindDataType(0, false /* signed */, true /* floating */, - true /* complex */), - GDT_CFloat32); - EXPECT_EQ(GDALFindDataType(0, true /* signed */, true /* floating */, - true /* complex */), - GDT_CFloat32); - - EXPECT_EQ(GDALFindDataType(8, false /* signed */, false /* floating */, - false /* complex */), - GDT_Byte); - EXPECT_EQ(GDALFindDataType(8, true /* signed */, false /* floating */, - false /* complex */), - GDT_Int8); - - EXPECT_EQ(GDALFindDataType(16, false /* signed */, false /* floating */, - false /* complex */), - GDT_UInt16); - EXPECT_EQ(GDALFindDataType(16, true /* signed */, false /* floating */, - false /* complex */), - GDT_Int16); - - EXPECT_EQ(GDALFindDataType(32, false /* signed */, false /* floating */, - false /* complex */), - GDT_UInt32); - EXPECT_EQ(GDALFindDataType(32, true /* signed */, false /* floating */, - false /* complex */), - GDT_Int32); + // Define brief abbreviations to make calls to `GDALFindDataType` + // more legible + const bool u = false, s = true; // signed + const bool i = false, f = true; // floating + const bool r = false, c = true; // complex - EXPECT_EQ(GDALFindDataType(64, false /* signed */, true /* floating */, - false /* complex */), - GDT_Float64); - EXPECT_EQ(GDALFindDataType(64, false /* signed */, true /* floating */, - true /* complex */), - GDT_CFloat64); + EXPECT_EQ(GDALFindDataType(0, u, i, r), GDT_Byte); + EXPECT_EQ(GDALFindDataType(0, s, i, r), GDT_Int8); + EXPECT_EQ(GDALFindDataType(0, u, i, c), GDT_CInt32); + EXPECT_EQ(GDALFindDataType(0, s, i, c), GDT_CInt16); + EXPECT_EQ(GDALFindDataType(0, u, f, r), GDT_Float32); + EXPECT_EQ(GDALFindDataType(0, s, f, r), GDT_Float32); + EXPECT_EQ(GDALFindDataType(0, u, f, c), GDT_CFloat32); + EXPECT_EQ(GDALFindDataType(0, s, f, c), GDT_CFloat32); - EXPECT_EQ(GDALFindDataType(64, false /* signed */, false /* floating */, - false /* complex */), - GDT_UInt64); - EXPECT_EQ(GDALFindDataType(64, true /* signed */, false /* floating */, - false /* complex */), - GDT_Int64); + EXPECT_EQ(GDALFindDataType(8, u, i, r), GDT_Byte); + EXPECT_EQ(GDALFindDataType(8, s, i, r), GDT_Int8); + + EXPECT_EQ(GDALFindDataType(16, u, f, r), GDT_Float32); + EXPECT_EQ(GDALFindDataType(16, u, f, c), GDT_CFloat32); + + EXPECT_EQ(GDALFindDataType(16, u, i, r), GDT_UInt16); + EXPECT_EQ(GDALFindDataType(16, s, i, r), GDT_Int16); + + EXPECT_EQ(GDALFindDataType(32, u, f, r), GDT_Float32); + EXPECT_EQ(GDALFindDataType(32, u, f, c), GDT_CFloat32); + + EXPECT_EQ(GDALFindDataType(32, u, i, r), GDT_UInt32); + EXPECT_EQ(GDALFindDataType(32, s, i, r), GDT_Int32); + + EXPECT_EQ(GDALFindDataType(64, u, f, r), GDT_Float64); + EXPECT_EQ(GDALFindDataType(64, u, f, c), GDT_CFloat64); + + EXPECT_EQ(GDALFindDataType(64, u, i, r), GDT_UInt64); + EXPECT_EQ(GDALFindDataType(64, s, i, r), GDT_Int64); EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -128, false), GDT_Int16); EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32768, false), GDT_Int16); @@ -285,30 +273,37 @@ TEST_F(test_gdal, GDALDataTypeUnion_special_cases) GDALDataTypeUnionWithValue(GDT_UInt64, 18446744073709555712.0, false), GDT_Float64); + EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float16, -999, false), + GDT_Float32); + EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float16, -99999, false), + GDT_Float32); + EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float16, -99999.9876, false), + GDT_Float64); + EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999, false), GDT_Float32); EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999.9876, false), GDT_Float64); EXPECT_EQ(GDALDataTypeUnionWithValue( - GDT_Float32, std::numeric_limits::quiet_NaN(), false), + GDT_Float32, cpl::NumericLimits::quiet_NaN(), false), GDT_Float32); EXPECT_EQ(GDALDataTypeUnionWithValue( - GDT_Float32, -std::numeric_limits::infinity(), false), + GDT_Float32, -cpl::NumericLimits::infinity(), false), GDT_Float32); EXPECT_EQ(GDALDataTypeUnionWithValue( - GDT_Float32, -std::numeric_limits::infinity(), false), + GDT_Float32, -cpl::NumericLimits::infinity(), false), GDT_Float32); EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float64, -99999.9876, false), GDT_Float64); EXPECT_EQ(GDALDataTypeUnionWithValue( - GDT_Float64, std::numeric_limits::quiet_NaN(), false), + GDT_Float64, cpl::NumericLimits::quiet_NaN(), false), GDT_Float64); EXPECT_EQ(GDALDataTypeUnionWithValue( - GDT_Float64, -std::numeric_limits::infinity(), false), + GDT_Float64, -cpl::NumericLimits::infinity(), false), GDT_Float64); EXPECT_EQ(GDALDataTypeUnionWithValue( - GDT_Float64, -std::numeric_limits::infinity(), false), + GDT_Float64, -cpl::NumericLimits::infinity(), false), GDT_Float64); EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Unknown, 0, false), GDT_Byte); @@ -419,6 +414,40 @@ TEST_F(test_gdal, GDALAdjustValueToDataType) &bRounded) == 10000000000.0 && !bClamped && !bRounded); + EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float16, 0.0, &bClamped, + &bRounded) == 0.0 && + !bClamped && !bRounded); + EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float16, 1e-10, &bClamped, + &bRounded) == 0.0 && + !bClamped && !bRounded); + EXPECT_TRUE( + GDALAdjustValueToDataType(GDT_Float16, 1.23, &bClamped, &bRounded) == + static_cast(static_cast(1.23f)) && + !bClamped && !bRounded); + EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float16, -1e300, &bClamped, + &bRounded) == -65504 && + bClamped && !bRounded); + EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float16, 1e30, &bClamped, + &bRounded) == 65504 && + bClamped && !bRounded); + EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float16, + cpl::NumericLimits::infinity(), + &bClamped, &bRounded) == + cpl::NumericLimits::infinity() && + !bClamped && !bRounded); + EXPECT_TRUE(GDALAdjustValueToDataType( + GDT_Float16, -cpl::NumericLimits::infinity(), + &bClamped, + &bRounded) == -cpl::NumericLimits::infinity() && + !bClamped && !bRounded); + { + double dfNan = cpl::NumericLimits::quiet_NaN(); + double dfGot = + GDALAdjustValueToDataType(GDT_Float16, dfNan, &bClamped, &bRounded); + EXPECT_TRUE(memcmp(&dfNan, &dfGot, sizeof(double)) == 0 && !bClamped && + !bRounded); + } + EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float32, 0.0, &bClamped, &bRounded) == 0.0 && !bClamped && !bRounded); @@ -431,24 +460,24 @@ TEST_F(test_gdal, GDALAdjustValueToDataType) !bClamped && !bRounded); EXPECT_TRUE( GDALAdjustValueToDataType(GDT_Float32, -1e300, &bClamped, &bRounded) == - -std::numeric_limits::max() && + -cpl::NumericLimits::max() && bClamped && !bRounded); EXPECT_TRUE( GDALAdjustValueToDataType(GDT_Float32, 1e300, &bClamped, &bRounded) == - std::numeric_limits::max() && + cpl::NumericLimits::max() && bClamped && !bRounded); - EXPECT_TRUE(GDALAdjustValueToDataType( - GDT_Float32, std::numeric_limits::infinity(), - &bClamped, - &bRounded) == std::numeric_limits::infinity() && + EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float32, + cpl::NumericLimits::infinity(), + &bClamped, &bRounded) == + cpl::NumericLimits::infinity() && !bClamped && !bRounded); EXPECT_TRUE(GDALAdjustValueToDataType( - GDT_Float32, -std::numeric_limits::infinity(), + GDT_Float32, -cpl::NumericLimits::infinity(), &bClamped, - &bRounded) == -std::numeric_limits::infinity() && + &bRounded) == -cpl::NumericLimits::infinity() && !bClamped && !bRounded); { - double dfNan = std::numeric_limits::quiet_NaN(); + double dfNan = cpl::NumericLimits::quiet_NaN(); double dfGot = GDALAdjustValueToDataType(GDT_Float32, dfNan, &bClamped, &bRounded); EXPECT_TRUE(memcmp(&dfNan, &dfGot, sizeof(double)) == 0 && !bClamped && @@ -467,18 +496,18 @@ TEST_F(test_gdal, GDALAdjustValueToDataType) EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float64, 1e40, &bClamped, &bRounded) == 1e40 && !bClamped && !bRounded); - EXPECT_TRUE(GDALAdjustValueToDataType( - GDT_Float64, std::numeric_limits::infinity(), - &bClamped, - &bRounded) == std::numeric_limits::infinity() && + EXPECT_TRUE(GDALAdjustValueToDataType(GDT_Float64, + cpl::NumericLimits::infinity(), + &bClamped, &bRounded) == + cpl::NumericLimits::infinity() && !bClamped && !bRounded); EXPECT_TRUE(GDALAdjustValueToDataType( - GDT_Float64, -std::numeric_limits::infinity(), + GDT_Float64, -cpl::NumericLimits::infinity(), &bClamped, - &bRounded) == -std::numeric_limits::infinity() && + &bRounded) == -cpl::NumericLimits::infinity() && !bClamped && !bRounded); { - double dfNan = std::numeric_limits::quiet_NaN(); + double dfNan = cpl::NumericLimits::quiet_NaN(); double dfGot = GDALAdjustValueToDataType(GDT_Float64, dfNan, &bClamped, &bRounded); EXPECT_TRUE(memcmp(&dfNan, &dfGot, sizeof(double)) == 0 && !bClamped && @@ -716,44 +745,44 @@ TEST_F(test_gdal, ARE_REAL_EQUAL) EXPECT_TRUE(!ARE_REAL_EQUAL(0.1, 0.0)); EXPECT_TRUE(ARE_REAL_EQUAL(1.0, 1.0)); EXPECT_TRUE(!ARE_REAL_EQUAL(1.0, 0.99)); - EXPECT_TRUE(ARE_REAL_EQUAL(-std::numeric_limits::min(), - -std::numeric_limits::min())); - EXPECT_TRUE(ARE_REAL_EQUAL(std::numeric_limits::min(), - std::numeric_limits::min())); - EXPECT_TRUE(!ARE_REAL_EQUAL(std::numeric_limits::min(), 0.0)); - EXPECT_TRUE(ARE_REAL_EQUAL(-std::numeric_limits::max(), - -std::numeric_limits::max())); - EXPECT_TRUE(ARE_REAL_EQUAL(std::numeric_limits::max(), - std::numeric_limits::max())); - EXPECT_TRUE(ARE_REAL_EQUAL(-std::numeric_limits::infinity(), - -std::numeric_limits::infinity())); - EXPECT_TRUE(ARE_REAL_EQUAL(std::numeric_limits::infinity(), - std::numeric_limits::infinity())); - EXPECT_TRUE(!ARE_REAL_EQUAL(std::numeric_limits::infinity(), - std::numeric_limits::max())); - EXPECT_TRUE(ARE_REAL_EQUAL(-std::numeric_limits::min(), - -std::numeric_limits::min())); + EXPECT_TRUE(ARE_REAL_EQUAL(-cpl::NumericLimits::min(), + -cpl::NumericLimits::min())); + EXPECT_TRUE(ARE_REAL_EQUAL(cpl::NumericLimits::min(), + cpl::NumericLimits::min())); + EXPECT_TRUE(!ARE_REAL_EQUAL(cpl::NumericLimits::min(), 0.0)); + EXPECT_TRUE(ARE_REAL_EQUAL(-cpl::NumericLimits::max(), + -cpl::NumericLimits::max())); + EXPECT_TRUE(ARE_REAL_EQUAL(cpl::NumericLimits::max(), + cpl::NumericLimits::max())); + EXPECT_TRUE(ARE_REAL_EQUAL(-cpl::NumericLimits::infinity(), + -cpl::NumericLimits::infinity())); + EXPECT_TRUE(ARE_REAL_EQUAL(cpl::NumericLimits::infinity(), + cpl::NumericLimits::infinity())); + EXPECT_TRUE(!ARE_REAL_EQUAL(cpl::NumericLimits::infinity(), + cpl::NumericLimits::max())); + EXPECT_TRUE(ARE_REAL_EQUAL(-cpl::NumericLimits::min(), + -cpl::NumericLimits::min())); EXPECT_TRUE(ARE_REAL_EQUAL(0.0f, 0.0f)); EXPECT_TRUE(!ARE_REAL_EQUAL(0.0f, 0.1f)); EXPECT_TRUE(!ARE_REAL_EQUAL(0.1f, 0.0f)); EXPECT_TRUE(ARE_REAL_EQUAL(1.0f, 1.0f)); EXPECT_TRUE(!ARE_REAL_EQUAL(1.0f, 0.99f)); - EXPECT_TRUE(ARE_REAL_EQUAL(-std::numeric_limits::min(), - -std::numeric_limits::min())); - EXPECT_TRUE(ARE_REAL_EQUAL(std::numeric_limits::min(), - std::numeric_limits::min())); - EXPECT_TRUE(!ARE_REAL_EQUAL(std::numeric_limits::min(), 0.0f)); - EXPECT_TRUE(ARE_REAL_EQUAL(-std::numeric_limits::max(), - -std::numeric_limits::max())); - EXPECT_TRUE(ARE_REAL_EQUAL(std::numeric_limits::max(), - std::numeric_limits::max())); - EXPECT_TRUE(ARE_REAL_EQUAL(-std::numeric_limits::infinity(), - -std::numeric_limits::infinity())); - EXPECT_TRUE(ARE_REAL_EQUAL(std::numeric_limits::infinity(), - std::numeric_limits::infinity())); - EXPECT_TRUE(!ARE_REAL_EQUAL(std::numeric_limits::infinity(), - std::numeric_limits::max())); + EXPECT_TRUE(ARE_REAL_EQUAL(-cpl::NumericLimits::min(), + -cpl::NumericLimits::min())); + EXPECT_TRUE(ARE_REAL_EQUAL(cpl::NumericLimits::min(), + cpl::NumericLimits::min())); + EXPECT_TRUE(!ARE_REAL_EQUAL(cpl::NumericLimits::min(), 0.0f)); + EXPECT_TRUE(ARE_REAL_EQUAL(-cpl::NumericLimits::max(), + -cpl::NumericLimits::max())); + EXPECT_TRUE(ARE_REAL_EQUAL(cpl::NumericLimits::max(), + cpl::NumericLimits::max())); + EXPECT_TRUE(ARE_REAL_EQUAL(-cpl::NumericLimits::infinity(), + -cpl::NumericLimits::infinity())); + EXPECT_TRUE(ARE_REAL_EQUAL(cpl::NumericLimits::infinity(), + cpl::NumericLimits::infinity())); + EXPECT_TRUE(!ARE_REAL_EQUAL(cpl::NumericLimits::infinity(), + cpl::NumericLimits::max())); } // Test GDALIsValueInRange() @@ -764,7 +793,7 @@ TEST_F(test_gdal, GDALIsValueInRange) EXPECT_TRUE(!GDALIsValueInRange(-1)); EXPECT_TRUE(!GDALIsValueInRange(256)); EXPECT_TRUE( - !GDALIsValueInRange(std::numeric_limits::quiet_NaN())); + !GDALIsValueInRange(cpl::NumericLimits::quiet_NaN())); EXPECT_TRUE(GDALIsValueInRange(-128)); EXPECT_TRUE(GDALIsValueInRange(127)); @@ -787,27 +816,25 @@ TEST_F(test_gdal, GDALIsValueInRange) EXPECT_TRUE(!GDALIsValueInRange(18446744073709551616.0)); EXPECT_TRUE(!GDALIsValueInRange(-0.5)); - EXPECT_TRUE(GDALIsValueInRange(-std::numeric_limits::max())); - EXPECT_TRUE(GDALIsValueInRange(std::numeric_limits::max())); + EXPECT_TRUE(GDALIsValueInRange(-cpl::NumericLimits::max())); + EXPECT_TRUE(GDALIsValueInRange(cpl::NumericLimits::max())); EXPECT_TRUE( - GDALIsValueInRange(-std::numeric_limits::infinity())); + GDALIsValueInRange(-cpl::NumericLimits::infinity())); EXPECT_TRUE( - GDALIsValueInRange(std::numeric_limits::infinity())); + GDALIsValueInRange(cpl::NumericLimits::infinity())); EXPECT_TRUE( - !GDALIsValueInRange(std::numeric_limits::quiet_NaN())); - EXPECT_TRUE( - !GDALIsValueInRange(-std::numeric_limits::max())); - EXPECT_TRUE(!GDALIsValueInRange(std::numeric_limits::max())); + !GDALIsValueInRange(cpl::NumericLimits::quiet_NaN())); + EXPECT_TRUE(!GDALIsValueInRange(-cpl::NumericLimits::max())); + EXPECT_TRUE(!GDALIsValueInRange(cpl::NumericLimits::max())); EXPECT_TRUE( - GDALIsValueInRange(-std::numeric_limits::infinity())); - EXPECT_TRUE( - GDALIsValueInRange(std::numeric_limits::infinity())); + GDALIsValueInRange(-cpl::NumericLimits::infinity())); EXPECT_TRUE( - GDALIsValueInRange(-std::numeric_limits::max())); - EXPECT_TRUE(GDALIsValueInRange(std::numeric_limits::max())); + GDALIsValueInRange(cpl::NumericLimits::infinity())); + EXPECT_TRUE(GDALIsValueInRange(-cpl::NumericLimits::max())); + EXPECT_TRUE(GDALIsValueInRange(cpl::NumericLimits::max())); EXPECT_TRUE( - !GDALIsValueInRange(std::numeric_limits::quiet_NaN())); + !GDALIsValueInRange(cpl::NumericLimits::quiet_NaN())); } #ifdef _MSC_VER @@ -827,7 +854,7 @@ TEST_F(test_gdal, GDALIsValueExactAs) EXPECT_TRUE(!GDALIsValueExactAs(255.5)); EXPECT_TRUE(!GDALIsValueExactAs(256)); EXPECT_TRUE( - !GDALIsValueExactAs(std::numeric_limits::quiet_NaN())); + !GDALIsValueExactAs(cpl::NumericLimits::quiet_NaN())); // -(1 << 63) EXPECT_TRUE(GDALIsValueExactAs(-9223372036854775808.0)); @@ -845,27 +872,25 @@ TEST_F(test_gdal, GDALIsValueExactAs) EXPECT_TRUE(!GDALIsValueExactAs(18446744073709551616.0)); EXPECT_TRUE(!GDALIsValueExactAs(-0.5)); - EXPECT_TRUE(GDALIsValueExactAs(-std::numeric_limits::max())); - EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits::max())); - EXPECT_TRUE( - GDALIsValueExactAs(-std::numeric_limits::infinity())); + EXPECT_TRUE(GDALIsValueExactAs(-cpl::NumericLimits::max())); + EXPECT_TRUE(GDALIsValueExactAs(cpl::NumericLimits::max())); EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::infinity())); + GDALIsValueExactAs(-cpl::NumericLimits::infinity())); EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::quiet_NaN())); + GDALIsValueExactAs(cpl::NumericLimits::infinity())); EXPECT_TRUE( - !GDALIsValueExactAs(-std::numeric_limits::max())); - EXPECT_TRUE(!GDALIsValueExactAs(std::numeric_limits::max())); + GDALIsValueExactAs(cpl::NumericLimits::quiet_NaN())); + EXPECT_TRUE(!GDALIsValueExactAs(-cpl::NumericLimits::max())); + EXPECT_TRUE(!GDALIsValueExactAs(cpl::NumericLimits::max())); EXPECT_TRUE( - GDALIsValueExactAs(-std::numeric_limits::infinity())); + GDALIsValueExactAs(-cpl::NumericLimits::infinity())); EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::infinity())); + GDALIsValueExactAs(cpl::NumericLimits::infinity())); + EXPECT_TRUE(GDALIsValueExactAs(-cpl::NumericLimits::max())); + EXPECT_TRUE(GDALIsValueExactAs(cpl::NumericLimits::max())); EXPECT_TRUE( - GDALIsValueExactAs(-std::numeric_limits::max())); - EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits::max())); - EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::quiet_NaN())); + GDALIsValueExactAs(cpl::NumericLimits::quiet_NaN())); } // Test GDALIsValueExactAs() @@ -895,74 +920,73 @@ TEST_F(test_gdal, GDALIsValueExactAs_C_func) EXPECT_FALSE(GDALIsValueExactAs(32768, GDT_Int16)); EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int16)); - EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits::lowest(), - GDT_UInt32)); EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::max(), GDT_UInt32)); + GDALIsValueExactAs(cpl::NumericLimits::lowest(), GDT_UInt32)); + EXPECT_TRUE( + GDALIsValueExactAs(cpl::NumericLimits::max(), GDT_UInt32)); EXPECT_FALSE(GDALIsValueExactAs( - std::numeric_limits::lowest() - 1.0, GDT_UInt32)); - EXPECT_FALSE(GDALIsValueExactAs(std::numeric_limits::max() + 1.0, + cpl::NumericLimits::lowest() - 1.0, GDT_UInt32)); + EXPECT_FALSE(GDALIsValueExactAs(cpl::NumericLimits::max() + 1.0, GDT_UInt32)); EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_UInt32)); EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::lowest(), GDT_Int32)); + GDALIsValueExactAs(cpl::NumericLimits::lowest(), GDT_Int32)); EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::max(), GDT_Int32)); - EXPECT_FALSE(GDALIsValueExactAs( - std::numeric_limits::lowest() - 1.0, GDT_Int32)); - EXPECT_FALSE(GDALIsValueExactAs(std::numeric_limits::max() + 1.0, + GDALIsValueExactAs(cpl::NumericLimits::max(), GDT_Int32)); + EXPECT_FALSE(GDALIsValueExactAs(cpl::NumericLimits::lowest() - 1.0, + GDT_Int32)); + EXPECT_FALSE(GDALIsValueExactAs(cpl::NumericLimits::max() + 1.0, GDT_Int32)); EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int32)); EXPECT_TRUE(GDALIsValueExactAs( - static_cast(std::numeric_limits::lowest()), + static_cast(cpl::NumericLimits::lowest()), GDT_UInt64)); // (1 << 64) - 2048 EXPECT_TRUE(GDALIsValueExactAs(18446744073709549568.0, GDT_UInt64)); EXPECT_FALSE(GDALIsValueExactAs( - static_cast(std::numeric_limits::lowest()) - 1.0, + static_cast(cpl::NumericLimits::lowest()) - 1.0, GDT_UInt64)); // (1 << 64) EXPECT_FALSE(GDALIsValueExactAs(18446744073709551616.0, GDT_UInt64)); EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_UInt64)); EXPECT_TRUE(GDALIsValueExactAs( - static_cast(std::numeric_limits::lowest()), - GDT_Int64)); + static_cast(cpl::NumericLimits::lowest()), GDT_Int64)); // (1 << 63) - 1024 EXPECT_TRUE(GDALIsValueExactAs(9223372036854774784.0, GDT_Int64)); EXPECT_FALSE(GDALIsValueExactAs( - static_cast(std::numeric_limits::lowest()) - 2048.0, + static_cast(cpl::NumericLimits::lowest()) - 2048.0, GDT_Int64)); // (1 << 63) - 512 EXPECT_FALSE(GDALIsValueExactAs(9223372036854775296.0, GDT_Int64)); EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int64)); EXPECT_TRUE( - GDALIsValueExactAs(-std::numeric_limits::max(), GDT_Float32)); + GDALIsValueExactAs(-cpl::NumericLimits::max(), GDT_Float32)); EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::max(), GDT_Float32)); - EXPECT_TRUE(GDALIsValueExactAs(-std::numeric_limits::infinity(), + GDALIsValueExactAs(cpl::NumericLimits::max(), GDT_Float32)); + EXPECT_TRUE(GDALIsValueExactAs(-cpl::NumericLimits::infinity(), GDT_Float32)); - EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits::infinity(), - GDT_Float32)); - EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits::quiet_NaN(), + EXPECT_TRUE( + GDALIsValueExactAs(cpl::NumericLimits::infinity(), GDT_Float32)); + EXPECT_TRUE(GDALIsValueExactAs(cpl::NumericLimits::quiet_NaN(), GDT_Float32)); EXPECT_TRUE( - !GDALIsValueExactAs(-std::numeric_limits::max(), GDT_Float32)); + !GDALIsValueExactAs(-cpl::NumericLimits::max(), GDT_Float32)); EXPECT_TRUE( - !GDALIsValueExactAs(std::numeric_limits::max(), GDT_Float32)); + !GDALIsValueExactAs(cpl::NumericLimits::max(), GDT_Float32)); - EXPECT_TRUE(GDALIsValueExactAs(-std::numeric_limits::infinity(), + EXPECT_TRUE(GDALIsValueExactAs(-cpl::NumericLimits::infinity(), GDT_Float64)); - EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits::infinity(), + EXPECT_TRUE(GDALIsValueExactAs(cpl::NumericLimits::infinity(), GDT_Float64)); EXPECT_TRUE( - GDALIsValueExactAs(-std::numeric_limits::max(), GDT_Float64)); + GDALIsValueExactAs(-cpl::NumericLimits::max(), GDT_Float64)); EXPECT_TRUE( - GDALIsValueExactAs(std::numeric_limits::max(), GDT_Float64)); - EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits::quiet_NaN(), + GDALIsValueExactAs(cpl::NumericLimits::max(), GDT_Float64)); + EXPECT_TRUE(GDALIsValueExactAs(cpl::NumericLimits::quiet_NaN(), GDT_Float64)); EXPECT_TRUE(GDALIsValueExactAs(0, GDT_CInt16)); @@ -1636,8 +1660,8 @@ TEST_F(test_gdal, GDALExtendedDataType) // Another error case of ProcessPerChunk { - const auto M64 = std::numeric_limits::max(); - const auto Msize_t = std::numeric_limits::max(); + const auto M64 = cpl::NumericLimits::max(); + const auto Msize_t = cpl::NumericLimits::max(); myArray ar(GDT_UInt16, {M64, M64, M64}, {32, 256, 128}); // Product of myCustomChunkSize[] > Msize_t @@ -1728,7 +1752,7 @@ TEST_F(test_gdal, GDALExtendedDataType) } } { - const auto M = std::numeric_limits::max(); + const auto M = cpl::NumericLimits::max(); myArray ar(GDT_UInt16, {M, M, M}, {M, M, M / 2}); { auto cs = ar.GetProcessingChunkSize(0); @@ -1744,15 +1768,15 @@ TEST_F(test_gdal, GDALExtendedDataType) } #if SIZEOF_VOIDP == 8 { - const auto M = std::numeric_limits::max(); + const auto M = cpl::NumericLimits::max(); myArray ar(GDT_UInt16, {M, M, M}, {M, M, M / 4}); { auto cs = - ar.GetProcessingChunkSize(std::numeric_limits::max()); + ar.GetProcessingChunkSize(cpl::NumericLimits::max()); EXPECT_EQ(cs.size(), 3U); EXPECT_EQ(cs[0], 1U); EXPECT_EQ(cs[1], 1U); - EXPECT_EQ(cs[2], (std::numeric_limits::max() / 4) * 2); + EXPECT_EQ(cs[2], (cpl::NumericLimits::max() / 4) * 2); } } #endif @@ -2946,7 +2970,7 @@ TEST_F(test_gdal, GDALBufferHasOnlyNoData) EXPECT_TRUE(!GDALBufferHasOnlyNoData(&float32val, 1e50, 1, 1, 1, 1, 32, GSF_FLOATING_POINT)); - float float32nan = std::numeric_limits::quiet_NaN(); + float float32nan = cpl::NumericLimits::quiet_NaN(); EXPECT_TRUE(GDALBufferHasOnlyNoData(&float32nan, float32nan, 1, 1, 1, 1, 32, GSF_FLOATING_POINT)); EXPECT_TRUE(!GDALBufferHasOnlyNoData(&float32nan, 0.0, 1, 1, 1, 1, 32, @@ -2958,7 +2982,7 @@ TEST_F(test_gdal, GDALBufferHasOnlyNoData) EXPECT_TRUE(!GDALBufferHasOnlyNoData(&float64val, 0.0, 1, 1, 1, 1, 64, GSF_FLOATING_POINT)); - double float64nan = std::numeric_limits::quiet_NaN(); + double float64nan = cpl::NumericLimits::quiet_NaN(); EXPECT_TRUE(GDALBufferHasOnlyNoData(&float64nan, float64nan, 1, 1, 1, 1, 64, GSF_FLOATING_POINT)); EXPECT_TRUE(!GDALBufferHasOnlyNoData(&float64nan, 0.0, 1, 1, 1, 1, 64, @@ -2970,185 +2994,220 @@ TEST_F(test_gdal, GetRasterNoDataReplacementValue) { // Test GDT_Byte EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Byte, std::numeric_limits::lowest()), + GDT_Byte, cpl::NumericLimits::lowest()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Byte, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Byte, std::numeric_limits::lowest()), - std::numeric_limits::lowest() + 1); - EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Byte, std::numeric_limits::max()), - std::numeric_limits::max() - 1); + GDT_Byte, cpl::NumericLimits::lowest()), + cpl::NumericLimits::lowest() + 1); + EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Byte, + cpl::NumericLimits::max()), + cpl::NumericLimits::max() - 1); // Test GDT_Int8 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int8, std::numeric_limits::lowest()), + GDT_Int8, cpl::NumericLimits::lowest()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Int8, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int8, std::numeric_limits::lowest()), - std::numeric_limits::lowest() + 1); + GDT_Int8, cpl::NumericLimits::lowest()), + cpl::NumericLimits::lowest() + 1); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Int8, - std::numeric_limits::max()), - std::numeric_limits::max() - 1); + cpl::NumericLimits::max()), + cpl::NumericLimits::max() - 1); // Test GDT_UInt16 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_UInt16, std::numeric_limits::lowest()), + GDT_UInt16, cpl::NumericLimits::lowest()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_UInt16, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_UInt16, std::numeric_limits::lowest()), - std::numeric_limits::lowest() + 1); + GDT_UInt16, cpl::NumericLimits::lowest()), + cpl::NumericLimits::lowest() + 1); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_UInt16, std::numeric_limits::max()), - std::numeric_limits::max() - 1); + GDT_UInt16, cpl::NumericLimits::max()), + cpl::NumericLimits::max() - 1); // Test GDT_Int16 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int16, std::numeric_limits::lowest()), + GDT_Int16, cpl::NumericLimits::lowest()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Int16, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int16, std::numeric_limits::lowest()), - std::numeric_limits::lowest() + 1); - EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int16, std::numeric_limits::max()), - std::numeric_limits::max() - 1); + GDT_Int16, cpl::NumericLimits::lowest()), + cpl::NumericLimits::lowest() + 1); + EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Int16, + cpl::NumericLimits::max()), + cpl::NumericLimits::max() - 1); // Test GDT_UInt32 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_UInt32, std::numeric_limits::lowest()), + GDT_UInt32, cpl::NumericLimits::lowest()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_UInt32, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_UInt32, std::numeric_limits::lowest()), - std::numeric_limits::lowest() + 1); + GDT_UInt32, cpl::NumericLimits::lowest()), + cpl::NumericLimits::lowest() + 1); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_UInt32, std::numeric_limits::max()), - std::numeric_limits::max() - 1); + GDT_UInt32, cpl::NumericLimits::max()), + cpl::NumericLimits::max() - 1); // Test GDT_Int32 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int32, std::numeric_limits::lowest()), + GDT_Int32, cpl::NumericLimits::lowest()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Int32, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int32, std::numeric_limits::lowest()), - std::numeric_limits::lowest() + 1); - EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int32, std::numeric_limits::max()), - std::numeric_limits::max() - 1); + GDT_Int32, cpl::NumericLimits::lowest()), + cpl::NumericLimits::lowest() + 1); + EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Int32, + cpl::NumericLimits::max()), + cpl::NumericLimits::max() - 1); // Test GDT_UInt64 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_UInt64, std::numeric_limits::lowest()), + GDT_UInt64, cpl::NumericLimits::lowest()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_UInt64, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue( GDT_UInt64, - static_cast(std::numeric_limits::lowest())), - static_cast(std::numeric_limits::lowest()) + 1); + static_cast(cpl::NumericLimits::lowest())), + static_cast(cpl::NumericLimits::lowest()) + 1); // uin64_t max is not representable in double so we expect the next value to be returned + using std::nextafter; EXPECT_EQ( GDALGetNoDataReplacementValue( GDT_UInt64, - static_cast(std::numeric_limits::max())), - std::nextafter( - static_cast(std::numeric_limits::max()), 0) - + static_cast(cpl::NumericLimits::max())), + nextafter(static_cast(cpl::NumericLimits::max()), 0) - 1); // Test GDT_Int64 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int64, std::numeric_limits::lowest()), + GDT_Int64, cpl::NumericLimits::lowest()), 0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Int64, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0); // in64_t max is not representable in double so we expect the next value to be returned EXPECT_EQ(GDALGetNoDataReplacementValue( GDT_Int64, - static_cast(std::numeric_limits::lowest())), - static_cast(std::numeric_limits::lowest()) + 1); - EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Int64, - static_cast(std::numeric_limits::max())), - std::nextafter( - static_cast(std::numeric_limits::max()), 0) - - 1); + static_cast(cpl::NumericLimits::lowest())), + static_cast(cpl::NumericLimits::lowest()) + 1); + EXPECT_EQ( + GDALGetNoDataReplacementValue( + GDT_Int64, static_cast(cpl::NumericLimits::max())), + nextafter(static_cast(cpl::NumericLimits::max()), 0) - + 1); // Test floating point types + // NOTE: Google Test's output for GFloat16 values is very wrong. + // It seems to round GFloat16 values to integers before outputting + // them. Do not trust the screen output when there is an error. + // However, the tests themselves are reliable. + + // out of range for float16 + EXPECT_EQ(GDALGetNoDataReplacementValue( + GDT_Float16, cpl::NumericLimits::lowest()), + 0.0); + EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Float16, + cpl::NumericLimits::max()), + 0.0); + EXPECT_EQ(GDALGetNoDataReplacementValue( + GDT_Float16, cpl::NumericLimits::infinity()), + 0.0); + EXPECT_EQ(GDALGetNoDataReplacementValue( + GDT_Float16, -cpl::NumericLimits::infinity()), + 0.0); + + // in range for float 16 + EXPECT_EQ( + static_cast(GDALGetNoDataReplacementValue(GDT_Float16, -1.0)), + nextafter(GFloat16(-1.0), GFloat16(0.0f))); + EXPECT_EQ( + static_cast(GDALGetNoDataReplacementValue(GDT_Float16, 1.1)), + nextafter(GFloat16(1.1), GFloat16(2.0f))); + EXPECT_EQ( + GDALGetNoDataReplacementValue(GDT_Float16, + cpl::NumericLimits::lowest()), + nextafter(cpl::NumericLimits::lowest(), GFloat16(0.0f))); + + EXPECT_EQ(GDALGetNoDataReplacementValue( + GDT_Float16, cpl::NumericLimits::max()), + static_cast(nextafter(cpl::NumericLimits::max(), + GFloat16(0.0f)))); + // out of range for float32 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Float32, std::numeric_limits::lowest()), + GDT_Float32, cpl::NumericLimits::lowest()), 0.0); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Float32, - std::numeric_limits::max()), + cpl::NumericLimits::max()), 0.0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Float32, std::numeric_limits::infinity()), + GDT_Float32, cpl::NumericLimits::infinity()), 0.0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Float32, -std::numeric_limits::infinity()), + GDT_Float32, -cpl::NumericLimits::infinity()), 0.0); // in range for float 32 EXPECT_EQ( static_cast(GDALGetNoDataReplacementValue(GDT_Float32, -1.0)), - std::nextafter(float(-1.0), 0.0f)); + nextafter(float(-1.0), 0.0f)); EXPECT_EQ( static_cast(GDALGetNoDataReplacementValue(GDT_Float32, 1.1)), - std::nextafter(float(1.1), 2.0f)); + nextafter(float(1.1), 2.0f)); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Float32, std::numeric_limits::lowest()), - std::nextafter(std::numeric_limits::lowest(), 0.0f)); + GDT_Float32, cpl::NumericLimits::lowest()), + nextafter(cpl::NumericLimits::lowest(), 0.0f)); - EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Float32, - std::numeric_limits::max()), - static_cast( - std::nextafter(std::numeric_limits::max(), 0.0f))); + EXPECT_EQ( + GDALGetNoDataReplacementValue(GDT_Float32, + cpl::NumericLimits::max()), + static_cast(nextafter(cpl::NumericLimits::max(), 0.0f))); // in range for float64 EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Float64, std::numeric_limits::lowest()), - std::nextafter(std::numeric_limits::lowest(), 0.0)); + GDT_Float64, cpl::NumericLimits::lowest()), + nextafter(cpl::NumericLimits::lowest(), 0.0)); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Float64, - std::numeric_limits::max()), - std::nextafter(std::numeric_limits::max(), 0.0)); + cpl::NumericLimits::max()), + nextafter(cpl::NumericLimits::max(), 0.0)); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Float64, std::numeric_limits::lowest()), - std::nextafter(std::numeric_limits::lowest(), 0.0)); + GDT_Float64, cpl::NumericLimits::lowest()), + nextafter(cpl::NumericLimits::lowest(), 0.0)); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Float64, - std::numeric_limits::max()), - std::nextafter(std::numeric_limits::max(), 0.0)); + cpl::NumericLimits::max()), + nextafter(cpl::NumericLimits::max(), 0.0)); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Float64, double(-1.0)), - std::nextafter(double(-1.0), 0.0)); + nextafter(double(-1.0), 0.0)); EXPECT_EQ(GDALGetNoDataReplacementValue(GDT_Float64, double(1.1)), - std::nextafter(double(1.1), 2.0)); + nextafter(double(1.1), 2.0)); // test infinity EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Float64, std::numeric_limits::infinity()), + GDT_Float64, cpl::NumericLimits::infinity()), 0.0); EXPECT_EQ(GDALGetNoDataReplacementValue( - GDT_Float64, -std::numeric_limits::infinity()), + GDT_Float64, -cpl::NumericLimits::infinity()), 0.0); } @@ -4754,12 +4813,12 @@ TEST_F(test_gdal, ReadRaster) std::vector res; EXPECT_EQ(poDS->GetRasterBand(1)->ReadRaster(res), CE_None); const auto expected_res = - std::vector{-std::numeric_limits::infinity(), + std::vector{-cpl::NumericLimits::infinity(), -1.0f, 1.0f, 128.0f, 32768.0f, - std::numeric_limits::infinity()}; + cpl::NumericLimits::infinity()}; EXPECT_EQ(res, expected_res); std::fill(res.begin(), res.end(), 0.0f); @@ -4772,12 +4831,12 @@ TEST_F(test_gdal, ReadRaster) std::vector> res; EXPECT_EQ(poDS->GetRasterBand(1)->ReadRaster(res), CE_None); const auto expected_res = std::vector>{ - -std::numeric_limits::infinity(), + -cpl::NumericLimits::infinity(), -1.0f, 1.0f, 128.0f, 32768.0f, - std::numeric_limits::infinity()}; + cpl::NumericLimits::infinity()}; EXPECT_EQ(res, expected_res); std::fill(res.begin(), res.end(), 0.0f); diff --git a/autotest/cpp/test_gdal_minmax_element.cpp b/autotest/cpp/test_gdal_minmax_element.cpp index b6d681dda76d..2c904b8158f1 100644 --- a/autotest/cpp/test_gdal_minmax_element.cpp +++ b/autotest/cpp/test_gdal_minmax_element.cpp @@ -698,8 +698,8 @@ TEST_F(test_gdal_minmax_element, float32) } { T nodata = 2.0f; - std::vector v{std::numeric_limits::quiet_NaN(), - std::numeric_limits::quiet_NaN(), nodata, max_v, + std::vector v{cpl::NumericLimits::quiet_NaN(), + cpl::NumericLimits::quiet_NaN(), nodata, max_v, min_v}; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata); EXPECT_EQ(v[idx_min], min_v); @@ -707,9 +707,9 @@ TEST_F(test_gdal_minmax_element, float32) EXPECT_EQ(v[idx_max], max_v); } { - T nodata = std::numeric_limits::quiet_NaN(); - std::vector v{std::numeric_limits::quiet_NaN(), - std::numeric_limits::quiet_NaN(), nodata, max_v, + T nodata = cpl::NumericLimits::quiet_NaN(); + std::vector v{cpl::NumericLimits::quiet_NaN(), + cpl::NumericLimits::quiet_NaN(), nodata, max_v, min_v}; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata); EXPECT_EQ(v[idx_min], min_v); @@ -717,26 +717,26 @@ TEST_F(test_gdal_minmax_element, float32) EXPECT_EQ(v[idx_max], max_v); } { - std::vector v{std::numeric_limits::quiet_NaN(), - std::numeric_limits::quiet_NaN(), + std::vector v{cpl::NumericLimits::quiet_NaN(), + cpl::NumericLimits::quiet_NaN(), max_v, - std::numeric_limits::quiet_NaN(), + cpl::NumericLimits::quiet_NaN(), min_v, - std::numeric_limits::quiet_NaN()}; + cpl::NumericLimits::quiet_NaN()}; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0); EXPECT_EQ(v[idx_min], min_v); auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0); EXPECT_EQ(v[idx_max], max_v); } { - std::vector v{max_v, std::numeric_limits::quiet_NaN(), min_v}; + std::vector v{max_v, cpl::NumericLimits::quiet_NaN(), min_v}; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0); EXPECT_EQ(v[idx_min], min_v); auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0); EXPECT_EQ(v[idx_max], max_v); } { - std::vector v(257, std::numeric_limits::quiet_NaN()); + std::vector v(257, cpl::NumericLimits::quiet_NaN()); v[125] = static_cast(min_v + 0.1f); v[126] = min_v; v[127] = static_cast(min_v + 0.1f); @@ -758,7 +758,7 @@ TEST_F(test_gdal_minmax_element, float32) EXPECT_EQ(v[idx_max], max_v); } { - std::vector v(255, std::numeric_limits::quiet_NaN()); + std::vector v(255, cpl::NumericLimits::quiet_NaN()); v[v.size() - 2] = min_v; v.back() = max_v; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0); @@ -818,8 +818,8 @@ TEST_F(test_gdal_minmax_element, float64) } { T nodata = 2.0; - std::vector v{std::numeric_limits::quiet_NaN(), - std::numeric_limits::quiet_NaN(), nodata, max_v, + std::vector v{cpl::NumericLimits::quiet_NaN(), + cpl::NumericLimits::quiet_NaN(), nodata, max_v, min_v}; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, true, nodata); EXPECT_EQ(v[idx_min], min_v); @@ -827,26 +827,26 @@ TEST_F(test_gdal_minmax_element, float64) EXPECT_EQ(v[idx_max], max_v); } { - std::vector v{max_v, std::numeric_limits::quiet_NaN(), min_v}; + std::vector v{max_v, cpl::NumericLimits::quiet_NaN(), min_v}; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0); EXPECT_EQ(v[idx_min], min_v); auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0); EXPECT_EQ(v[idx_max], max_v); } { - std::vector v{std::numeric_limits::quiet_NaN(), - std::numeric_limits::quiet_NaN(), + std::vector v{cpl::NumericLimits::quiet_NaN(), + cpl::NumericLimits::quiet_NaN(), max_v, - std::numeric_limits::quiet_NaN(), + cpl::NumericLimits::quiet_NaN(), min_v, - std::numeric_limits::quiet_NaN()}; + cpl::NumericLimits::quiet_NaN()}; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0); EXPECT_EQ(v[idx_min], min_v); auto idx_max = gdal::max_element(v.data(), v.size(), eDT, false, 0); EXPECT_EQ(v[idx_max], max_v); } { - std::vector v(33, std::numeric_limits::quiet_NaN()); + std::vector v(33, cpl::NumericLimits::quiet_NaN()); v[5] = min_v; v[15] = max_v; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0); @@ -855,7 +855,7 @@ TEST_F(test_gdal_minmax_element, float64) EXPECT_EQ(v[idx_max], max_v); } { - std::vector v(255, std::numeric_limits::quiet_NaN()); + std::vector v(255, cpl::NumericLimits::quiet_NaN()); v[v.size() - 2] = min_v; v.back() = max_v; auto idx_min = gdal::min_element(v.data(), v.size(), eDT, false, 0); diff --git a/autotest/cpp/testcopywords.cpp b/autotest/cpp/testcopywords.cpp index 2fbdcdd3f50b..c843476358f1 100644 --- a/autotest/cpp/testcopywords.cpp +++ b/autotest/cpp/testcopywords.cpp @@ -11,6 +11,7 @@ ****************************************************************************/ #include "cpl_conv.h" +#include "cpl_float.h" #include "gdal.h" #include @@ -158,6 +159,9 @@ class TestCopyWords : public ::testing::Test else if (outtype == GDT_UInt64) Test( intype, inval, invali, outtype, outval, outvali, numLine); + else if (outtype == GDT_Float16) + Test(intype, inval, invali, outtype, + outval, outvali, numLine); else if (outtype == GDT_Float32) Test(intype, inval, invali, outtype, outval, outvali, numLine); @@ -170,6 +174,9 @@ class TestCopyWords : public ::testing::Test else if (outtype == GDT_CInt32) Test(intype, inval, invali, outtype, outval, outvali, numLine); + else if (outtype == GDT_CFloat16) + Test(intype, inval, invali, outtype, + outval, outvali, numLine); else if (outtype == GDT_CFloat32) Test(intype, inval, invali, outtype, outval, outvali, numLine); @@ -207,6 +214,9 @@ class TestCopyWords : public ::testing::Test else if (intype == GDT_UInt64) FromR_2(intype, inval, invali, outtype, outval, outvali, numLine); + else if (intype == GDT_Float16) + FromR_2(intype, inval, invali, outtype, + outval, outvali, numLine); else if (intype == GDT_Float32) FromR_2(intype, inval, invali, outtype, outval, outvali, numLine); @@ -219,6 +229,9 @@ class TestCopyWords : public ::testing::Test else if (intype == GDT_CInt32) FromR_2(intype, inval, invali, outtype, outval, outvali, numLine); + else if (intype == GDT_CFloat16) + FromR_2(intype, inval, invali, outtype, + outval, outvali, numLine); else if (intype == GDT_CFloat32) FromR_2(intype, inval, invali, outtype, outval, outvali, numLine); @@ -241,8 +254,8 @@ class TestCopyWords : public ::testing::Test #define IS_UNSIGNED(x) \ (x == GDT_Byte || x == GDT_UInt16 || x == GDT_UInt32 || x == GDT_UInt64) #define IS_FLOAT(x) \ - (x == GDT_Float32 || x == GDT_Float64 || x == GDT_CFloat32 || \ - x == GDT_CFloat64) + (x == GDT_Float16 || x == GDT_Float32 || x == GDT_Float64 || \ + x == GDT_CFloat16 || x == GDT_CFloat32 || x == GDT_CFloat64) #define CST_3000000000 (((GIntBig)3000) * 1000 * 1000) #define CST_5000000000 (((GIntBig)5000) * 1000 * 1000) @@ -290,10 +303,12 @@ TEST_F(TestCopyWords, GDT_Int8) FROM_R(GDT_Int8, -128, GDT_UInt32, 0); /* clamp */ FROM_R(GDT_Int8, -128, GDT_Int64, -128); FROM_R(GDT_Int8, -128, GDT_UInt64, 0); /* clamp */ + FROM_R(GDT_Int8, -128, GDT_Float16, -128); FROM_R(GDT_Int8, -128, GDT_Float32, -128); FROM_R(GDT_Int8, -128, GDT_Float64, -128); FROM_R(GDT_Int8, -128, GDT_CInt16, -128); FROM_R(GDT_Int8, -128, GDT_CInt32, -128); + FROM_R(GDT_Int8, -128, GDT_CFloat16, -128); FROM_R(GDT_Int8, -128, GDT_CFloat32, -128); FROM_R(GDT_Int8, -128, GDT_CFloat64, -128); for (GDALDataType outtype = GDT_Byte; outtype < GDT_TypeCount; @@ -507,6 +522,144 @@ TEST_F(TestCopyWords, GDT_UInt64) FROM_R(GDT_UInt64, nVal, GDT_CFloat64, nVal); } +TEST_F(TestCopyWords, GDT_Float16only) +{ + GDALDataType intype = GDT_Float16; + for (GDALDataType outtype = GDT_Byte; outtype < GDT_TypeCount; + outtype = (GDALDataType)(outtype + 1)) + { + if (IS_FLOAT(outtype)) + { + FROM_R_F(intype, 127.1, outtype, 127.1); + FROM_R_F(intype, -127.1, outtype, -127.1); + } + else + { + FROM_R_F(intype, 125.1, outtype, 125); + FROM_R_F(intype, 125.9, outtype, 126); + + FROM_R_F(intype, 0.4, outtype, 0); + FROM_R_F(intype, 0.5, outtype, + 1); /* We could argue how to do this rounding */ + FROM_R_F(intype, 0.6, outtype, 1); + FROM_R_F(intype, 126.5, outtype, + 127); /* We could argue how to do this rounding */ + + if (!IS_UNSIGNED(outtype)) + { + FROM_R_F(intype, -125.9, outtype, -126); + FROM_R_F(intype, -127.1, outtype, -127); + + FROM_R_F(intype, -0.4, outtype, 0); + FROM_R_F(intype, -0.5, outtype, + -1); /* We could argue how to do this rounding */ + FROM_R_F(intype, -0.6, outtype, -1); + FROM_R_F(intype, -127.5, outtype, + -128); /* We could argue how to do this rounding */ + } + } + } + FROM_R(intype, -30000, GDT_Byte, 0); + FROM_R(intype, -32768, GDT_Byte, 0); + FROM_R(intype, -1, GDT_Byte, 0); + FROM_R(intype, 256, GDT_Byte, 255); + FROM_R(intype, 30000, GDT_Byte, 255); + FROM_R(intype, -330000, GDT_Int16, -32768); + FROM_R(intype, -33000, GDT_Int16, -32768); + FROM_R(intype, 33000, GDT_Int16, 32767); + FROM_R(intype, -33000, GDT_UInt16, 0); + FROM_R(intype, -1, GDT_UInt16, 0); + FROM_R(intype, 60000, GDT_UInt16, 60000); + FROM_R(intype, -33000, GDT_Int32, -32992); + FROM_R(intype, 33000, GDT_Int32, 32992); + FROM_R(intype, -1, GDT_UInt32, 0); + FROM_R(intype, 60000, GDT_UInt32, 60000U); + FROM_R(intype, 33000, GDT_Float32, 32992); + FROM_R(intype, -33000, GDT_Float32, -32992); + FROM_R(intype, 33000, GDT_Float64, 32992); + FROM_R(intype, -33000, GDT_Float64, -32992); + FROM_R(intype, -33000, GDT_CInt16, -32768); + FROM_R(intype, 33000, GDT_CInt16, 32767); + FROM_R(intype, -33000, GDT_CInt32, -32992); + FROM_R(intype, 33000, GDT_CInt32, 32992); + FROM_R(intype, 33000, GDT_CFloat32, 32992); + FROM_R(intype, -33000, GDT_CFloat32, -32992); + FROM_R(intype, 33000, GDT_CFloat64, 32992); + FROM_R(intype, -33000, GDT_CFloat64, -32992); + + // Float16 to Int64 + { + GFloat16 in_value = cpl::NumericLimits::quiet_NaN(); + int64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_Int64, 0, 1); + EXPECT_EQ(out_value, 0); + } + + { + GFloat16 in_value = -cpl::NumericLimits::infinity(); + int64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_Int64, 0, 1); + EXPECT_EQ(out_value, INT64_MIN); + } + + { + GFloat16 in_value = -cpl::NumericLimits::max(); + int64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_Int64, 0, 1); + EXPECT_EQ(out_value, -65504); + } + + { + GFloat16 in_value = cpl::NumericLimits::max(); + int64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_Int64, 0, 1); + EXPECT_EQ(out_value, 65504); + } + + { + GFloat16 in_value = cpl::NumericLimits::infinity(); + int64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_Int64, 0, 1); + EXPECT_EQ(out_value, INT64_MAX); + } + + // Float16 to UInt64 + { + GFloat16 in_value = cpl::NumericLimits::quiet_NaN(); + uint64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_UInt64, 0, 1); + EXPECT_EQ(out_value, 0); + } + + { + GFloat16 in_value = -cpl::NumericLimits::infinity(); + uint64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_UInt64, 0, 1); + EXPECT_EQ(out_value, 0); + } + + { + GFloat16 in_value = -cpl::NumericLimits::max(); + uint64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_UInt64, 0, 1); + EXPECT_EQ(out_value, 0); + } + + { + GFloat16 in_value = cpl::NumericLimits::max(); + uint64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_UInt64, 0, 1); + EXPECT_EQ(out_value, 65504); + } + + { + GFloat16 in_value = cpl::NumericLimits::infinity(); + uint64_t out_value = 0; + GDALCopyWords(&in_value, GDT_Float16, 0, &out_value, GDT_UInt64, 0, 1); + EXPECT_EQ(out_value, UINT64_MAX); + } +} + TEST_F(TestCopyWords, GDT_Float32and64) { /* GDT_Float32 and GDT_Float64 */ @@ -581,35 +734,35 @@ TEST_F(TestCopyWords, GDT_Float32and64) // Float32 to Int64 { - float in_value = std::numeric_limits::quiet_NaN(); + float in_value = cpl::NumericLimits::quiet_NaN(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, 0); } { - float in_value = -std::numeric_limits::infinity(); + float in_value = -cpl::NumericLimits::infinity(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, INT64_MIN); } { - float in_value = -std::numeric_limits::max(); + float in_value = -cpl::NumericLimits::max(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, INT64_MIN); } { - float in_value = std::numeric_limits::max(); + float in_value = cpl::NumericLimits::max(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, INT64_MAX); } { - float in_value = std::numeric_limits::infinity(); + float in_value = cpl::NumericLimits::infinity(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, INT64_MAX); @@ -617,35 +770,35 @@ TEST_F(TestCopyWords, GDT_Float32and64) // Float64 to Int64 { - double in_value = std::numeric_limits::quiet_NaN(); + double in_value = cpl::NumericLimits::quiet_NaN(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, 0); } { - double in_value = -std::numeric_limits::infinity(); + double in_value = -cpl::NumericLimits::infinity(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, INT64_MIN); } { - double in_value = -std::numeric_limits::max(); + double in_value = -cpl::NumericLimits::max(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, INT64_MIN); } { - double in_value = std::numeric_limits::max(); + double in_value = cpl::NumericLimits::max(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, INT64_MAX); } { - double in_value = std::numeric_limits::infinity(); + double in_value = cpl::NumericLimits::infinity(); int64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_Int64, 0, 1); EXPECT_EQ(out_value, INT64_MAX); @@ -653,35 +806,35 @@ TEST_F(TestCopyWords, GDT_Float32and64) // Float32 to UInt64 { - float in_value = std::numeric_limits::quiet_NaN(); + float in_value = cpl::NumericLimits::quiet_NaN(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, 0); } { - float in_value = -std::numeric_limits::infinity(); + float in_value = -cpl::NumericLimits::infinity(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, 0); } { - float in_value = -std::numeric_limits::max(); + float in_value = -cpl::NumericLimits::max(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, 0); } { - float in_value = std::numeric_limits::max(); + float in_value = cpl::NumericLimits::max(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, UINT64_MAX); } { - float in_value = std::numeric_limits::infinity(); + float in_value = cpl::NumericLimits::infinity(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float32, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, UINT64_MAX); @@ -689,35 +842,35 @@ TEST_F(TestCopyWords, GDT_Float32and64) // Float64 to UInt64 { - double in_value = -std::numeric_limits::quiet_NaN(); + double in_value = -cpl::NumericLimits::quiet_NaN(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, 0); } { - double in_value = -std::numeric_limits::infinity(); + double in_value = -cpl::NumericLimits::infinity(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, 0); } { - double in_value = -std::numeric_limits::max(); + double in_value = -cpl::NumericLimits::max(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, 0); } { - double in_value = std::numeric_limits::max(); + double in_value = cpl::NumericLimits::max(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, UINT64_MAX); } { - double in_value = std::numeric_limits::infinity(); + double in_value = cpl::NumericLimits::infinity(); uint64_t out_value = 0; GDALCopyWords(&in_value, GDT_Float64, 0, &out_value, GDT_UInt64, 0, 1); EXPECT_EQ(out_value, UINT64_MAX); @@ -842,6 +995,49 @@ TEST_F(TestCopyWords, GDT_CFloat32and64) } } +TEST_F(TestCopyWords, GDT_CFloat16only) +{ + GDALDataType intype = GDT_CFloat16; + for (GDALDataType outtype = GDT_Byte; outtype < GDT_TypeCount; + outtype = (GDALDataType)(outtype + 1)) + { + if (IS_FLOAT(outtype)) + { + FROM_C_F(intype, 127.1, 127.9, outtype, 127.1, 127.9); + FROM_C_F(intype, -127.1, -127.9, outtype, -127.1, -127.9); + } + else + { + FROM_C_F(intype, 126.1, 150.9, outtype, 126, 151); + FROM_C_F(intype, 126.9, 150.1, outtype, 127, 150); + if (!IS_UNSIGNED(outtype)) + { + FROM_C_F(intype, -125.9, -127.1, outtype, -126, -127); + } + } + } + FROM_C(intype, -1, 256, GDT_Byte, 0, 0); + FROM_C(intype, 256, -1, GDT_Byte, 255, 0); + FROM_C(intype, -33000, 33000, GDT_Int16, -32768, 0); + FROM_C(intype, 33000, -33000, GDT_Int16, 32767, 0); + FROM_C(intype, -1, 66000, GDT_UInt16, 0, 0); + FROM_C(intype, 66000, -1, GDT_UInt16, 65535, 0); + FROM_C(intype, -33000, -33000, GDT_Int32, -32992, 0); + FROM_C(intype, 33000, 33000, GDT_Int32, 32992, 0); + FROM_C(intype, -1, 33000, GDT_UInt32, 0, 0); + FROM_C(intype, 33000, -1, GDT_UInt32, 32992, 0); + FROM_C(intype, 33000, -1, GDT_Float32, 32992, 0); + FROM_C(intype, 33000, -1, GDT_Float64, 32992, 0); + FROM_C(intype, -33000, -1, GDT_Float32, -32992, 0); + FROM_C(intype, -33000, -1, GDT_Float64, -32992, 0); + FROM_C(intype, -33000, 33000, GDT_CInt16, -32768, 32767); + FROM_C(intype, 33000, -33000, GDT_CInt16, 32767, -32768); + FROM_C(intype, -33000, -33000, GDT_CInt32, -32992, -32992); + FROM_C(intype, 33000, 33000, GDT_CInt32, 32992, 32992); + FROM_C(intype, 33000, -33000, GDT_CFloat32, 32992, -32992); + FROM_C(intype, 33000, -33000, GDT_CFloat64, 32992, -32992); +} + template void CheckPackedGeneric(GDALDataType eIn, GDALDataType eOut) { @@ -945,6 +1141,9 @@ template void CheckPacked(GDALDataType eIn, GDALDataType eOut) case GDT_Int64: CheckPacked(eIn, eOut); break; + case GDT_Float16: + CheckPacked(eIn, eOut); + break; case GDT_Float32: CheckPacked(eIn, eOut); break; @@ -984,6 +1183,9 @@ static void CheckPacked(GDALDataType eIn, GDALDataType eOut) case GDT_Int64: CheckPacked(eIn, eOut); break; + case GDT_Float16: + CheckPacked(eIn, eOut); + break; case GDT_Float32: CheckPacked(eIn, eOut); break; diff --git a/autotest/cpp/testfloat16.cpp b/autotest/cpp/testfloat16.cpp new file mode 100644 index 000000000000..7d141a32c1b4 --- /dev/null +++ b/autotest/cpp/testfloat16.cpp @@ -0,0 +1,187 @@ +/****************************************************************************** + * $Id$ + * + * Project: GDAL Core + * Purpose: Test GFloat16. + * Author: Erik Schnetter + * + ****************************************************************************** + * Copyright (c) 2024, Even Rouault + * + * SPDX-License-Identifier: MIT + ****************************************************************************/ + +#include "cpl_float.h" + +#include "gtest_include.h" + +#include +#include + +namespace +{ + +TEST(TestFloat16, conversions) +{ + for (int i = -2048; i <= +2048; ++i) + { + EXPECT_EQ(GFloat16(i), static_cast(i)); + if (i >= -128 && i <= 127) + { + EXPECT_EQ(GFloat16(static_cast(i)), + static_cast(i)); + } + EXPECT_EQ(GFloat16(static_cast(i)), static_cast(i)); + EXPECT_EQ(GFloat16(static_cast(i)), static_cast(i)); + EXPECT_EQ(GFloat16(static_cast(i)), static_cast(i)); + EXPECT_EQ(GFloat16(static_cast(i)), static_cast(i)); + if (i >= 0) + { + if (i <= 255) + { + EXPECT_EQ(GFloat16(static_cast(i)), + static_cast(i)); + } + EXPECT_EQ(GFloat16(static_cast(i)), + static_cast(i)); + EXPECT_EQ(GFloat16(static_cast(i)), + static_cast(i)); + EXPECT_EQ(GFloat16(static_cast(i)), + static_cast(i)); + EXPECT_EQ(GFloat16(static_cast(i)), + static_cast(i)); + } + EXPECT_EQ(GFloat16(i), GFloat16(i)); + EXPECT_EQ(GFloat16(i), static_cast(i)); + } + + EXPECT_EQ(GFloat16(65504), 65504.0); + EXPECT_EQ(GFloat16(-65504), -65504.0); + // Work around the Windows compiler reporting "error C2124: divide + // or mod by zero". See also + // . + volatile double zero = 0.0; + EXPECT_EQ(GFloat16(1.0 / zero), 1.0 / zero); + EXPECT_EQ(GFloat16(-1.0 / zero), -1.0 / zero); + EXPECT_EQ(GFloat16(0.0), -0.0); + EXPECT_EQ(GFloat16(-0.0), 0.0); +} + +TEST(TestFloat16, arithmetic) +{ + for (int i = -100; i <= +100; ++i) + { + double x = i; + + EXPECT_EQ(+GFloat16(x), +x); + EXPECT_EQ(-GFloat16(x), -x); + } + + for (int i = -100; i <= +100; ++i) + { + for (int j = -100; j <= +100; ++j) + { + double x = i; + double y = j; + + EXPECT_EQ(GFloat16(x) + GFloat16(y), x + y); + EXPECT_EQ(GFloat16(x) - GFloat16(y), x - y); + using std::fabs; + EXPECT_NEAR(GFloat16(x) * GFloat16(y), x * y, fabs(x * y / 1024)); + if (j != 0) + { + EXPECT_NEAR(GFloat16(x) / GFloat16(y), x / y, + fabs(x / y / 1024)); + } + } + } +} + +TEST(TestFloat16, comparisons) +{ + for (int i = -100; i <= +100; ++i) + { + for (int j = -100; j <= +100; ++j) + { + double x = i; + double y = j; + + EXPECT_EQ(GFloat16(x) == GFloat16(y), x == y); + EXPECT_EQ(GFloat16(x) != GFloat16(y), x != y); + EXPECT_EQ(GFloat16(x) < GFloat16(y), x < y); + EXPECT_EQ(GFloat16(x) > GFloat16(y), x > y); + EXPECT_EQ(GFloat16(x) <= GFloat16(y), x <= y); + EXPECT_EQ(GFloat16(x) >= GFloat16(y), x >= y); + } + } +} + +TEST(TestFloat16, math) +{ + for (int i = -100; i <= +100; ++i) + { + double x = i; + + using std::isfinite; + EXPECT_EQ(isfinite(GFloat16(x)), isfinite(x)); + using std::isinf; + EXPECT_EQ(isinf(GFloat16(x)), isinf(x)); + using std::isnan; + EXPECT_EQ(isnan(GFloat16(x)), isnan(x)); + using std::abs; + EXPECT_EQ(abs(GFloat16(x)), abs(x)); + using std::cbrt; + using std::fabs; + EXPECT_NEAR(cbrt(GFloat16(x)), cbrt(x), fabs(cbrt(x) / 1024)); + using std::ceil; + EXPECT_EQ(ceil(GFloat16(x)), ceil(x)); + using std::fabs; + EXPECT_EQ(fabs(GFloat16(x)), fabs(x)); + using std::floor; + EXPECT_EQ(floor(GFloat16(x)), floor(x)); + using std::round; + EXPECT_EQ(round(GFloat16(x)), round(x)); + if (x >= 0) + { + using std::sqrt; + EXPECT_NEAR(sqrt(GFloat16(x)), sqrt(x), fabs(sqrt(x) / 1024)); + } + } + + for (int i = -100; i <= +100; ++i) + { + for (int j = -100; j <= +100; ++j) + { + double x = i; + double y = j; + + using std::fmax; + EXPECT_EQ(fmax(GFloat16(x), GFloat16(y)), GFloat16(fmax(x, y))); + using std::fmin; + EXPECT_EQ(fmin(GFloat16(x), GFloat16(y)), GFloat16(fmin(x, y))); + using std::hypot; + EXPECT_EQ(hypot(GFloat16(x), GFloat16(y)), GFloat16(hypot(x, y))); + using std::max; + EXPECT_EQ(max(GFloat16(x), GFloat16(y)), GFloat16(max(x, y))); + using std::min; + EXPECT_EQ(min(GFloat16(x), GFloat16(y)), GFloat16(min(x, y))); + using std::pow; + EXPECT_EQ(pow(GFloat16(x), GFloat16(y)), GFloat16(pow(x, y))); + using std::fabs; + using std::isfinite; + GFloat16 r1 = GFloat16(pow(GFloat16(x), j)); + GFloat16 r2 = GFloat16(pow(x, j)); + if (!isfinite(r1)) + { + EXPECT_EQ(r1, r2); + } + else + { + GFloat16 tol = (1 + fabs(r2)) / 1024; + EXPECT_NEAR(r1, r2, tol); + } + } + } +} + +} // namespace diff --git a/autotest/gcore/numpy_rw.py b/autotest/gcore/numpy_rw.py index 666146fbe8dc..2b2dd099d4d3 100755 --- a/autotest/gcore/numpy_rw.py +++ b/autotest/gcore/numpy_rw.py @@ -658,7 +658,7 @@ def test_numpy_rw_16(): assert ds is None # Unsupported data type - array = numpy.empty([1, 1], numpy.float16) + array = numpy.empty([1, 1], numpy.bool_) with gdal.quiet_errors(): ds = gdal_array.OpenArray(array) assert ds is None diff --git a/autotest/gdrivers/vrtderived.py b/autotest/gdrivers/vrtderived.py index 0547ca0e33a0..48cb238dbfdb 100755 --- a/autotest/gdrivers/vrtderived.py +++ b/autotest/gdrivers/vrtderived.py @@ -751,14 +751,19 @@ def test_vrtderived_12(): for dt in [ "Byte", + "Int8", "UInt16", "Int16", "UInt32", "Int32", + "UInt64", + "Int64", + "Float16", "Float32", "Float64", "CInt16", "CInt32", + "CFloat16", "CFloat32", "CFloat64", ]: @@ -777,18 +782,30 @@ def test_vrtderived_12(): "GDAL_VRT_ENABLE_PYTHON", "YES" ), gdaltest.error_handler(): cs = ds.GetRasterBand(1).Checksum() - # CInt16/CInt32 do not map to native numpy types - if dt == "CInt16" or dt == "CInt32": - expected_cs = -1 # error + # CInt16/CInt32/CFloat16 do not map to native numpy types + if dt == "CInt16" or dt == "CInt32" or dt == "CFloat16": + expected_cs = [-1] # error + elif dt == "Float16": + # Might or might not be supported by GDAL + expected_cs = [-1, 100] else: - expected_cs = 100 - if cs != expected_cs: + expected_cs = [100] + if cs not in expected_cs: print(dt) print(gdal.GetLastErrorMsg()) - pytest.fail("invalid checksum") + if len(expected_cs) == 1: + pytest.fail( + "invalid checksum, datatype %s, have %d, expected %d" + % (dt, cs, expected_cs[0]) + ) + else: + pytest.fail( + "invalid checksum, datatype %s, have %d, expected one of [%d, %d]" + % (dt, cs, expected_cs[0], expected_cs[1]) + ) # Same for SourceTransferType - for dt in ["CInt16", "CInt32"]: + for dt in ["CInt16", "CInt32", "CFloat16"]: ds = gdal.Open( """ @@ -1050,7 +1067,14 @@ def identity(in_ar, out_ar, *args, **kwargs): with gdal.config_option("GDAL_VRT_ENABLE_PYTHON", "YES"): with gdal.Open(vrt_xml) as vrt_ds: arr = vrt_ds.ReadAsArray() - if dtype not in {gdal.GDT_CInt16, gdal.GDT_CInt32}: + # The complex int/float types are not available in numpy. + # Float16 may or may not be supported by GDAL. + if dtype not in { + gdal.GDT_CInt16, + gdal.GDT_CInt32, + gdal.GDT_Float16, + gdal.GDT_CFloat16, + }: assert arr[0, 0] == 1 assert vrt_ds.GetRasterBand(1).DataType == dtype diff --git a/autotest/gdrivers/zarr_driver.py b/autotest/gdrivers/zarr_driver.py index caed01d26f2e..e3292285a128 100644 --- a/autotest/gdrivers/zarr_driver.py +++ b/autotest/gdrivers/zarr_driver.py @@ -42,8 +42,10 @@ def module_disable_exceptions(): gdal.GDT_UInt32: "I", gdal.GDT_Int64: "q", gdal.GDT_UInt64: "Q", + gdal.GDT_Float16: "e", gdal.GDT_Float32: "f", gdal.GDT_Float64: "d", + gdal.GDT_CFloat16: "e", gdal.GDT_CFloat32: "f", gdal.GDT_CFloat64: "d", } @@ -87,6 +89,13 @@ def module_disable_exceptions(): (1 << 64) - 1, ], # not really legit to have the fill_value as a str, but libjson-c can't support numeric values in int64::max(), uint64::max() range. [">u8", gdal.GDT_UInt64, None, None], + # We would like to test these, but SWIG does not support float16 (yet?) + # ["f2", gdal.GDT_Float16, None, None], + # ["f4", gdal.GDT_Float32, None, None], ["c4", gdal.GDT_CFloat16, None, None], ["c8", gdal.GDT_CFloat32, None, None], ["u8", gdal.GDT_Float64, None, None], + # TODO: Test reading/writing GDT_Float16 via float32 Python data + # ["f2", gdal.GDT_Float16, None, None], + # ["f4", gdal.GDT_Float32, None, None], ["c4", gdal.GDT_CFloat16, None, None], ["c8", gdal.GDT_CFloat32, None, None], [" tmp.txt + (sed 's/fileLocation/artifactLocation/g' < $f) |jq '.runs[].results[] | (if .locations[].physicalLocation.artifactLocation.uri | (contains("/usr/include") or contains("degrib") or contains("libpng") or contains("libjpeg") or contains("EHapi") or contains("GDapi") or contains("SWapi") or contains("osr_cs_wkt_parser") or contains("ods_formula_parser") or contains("swq_parser") or contains("libjson") or contains("flatbuffers") or contains("cpl_float.h") or contains("cpl_minizip_zip.cpp") or contains("gdal_rpc.cpp") or contains("gdal_interpolateatpoint.cpp") or contains("internal_libqhull") ) then empty else { "uri": .locations[].physicalLocation.artifactLocation.uri, "msg": .message.text, "location": .codeFlows[-1].threadFlows[-1].locations[-1] } end)' > tmp.txt if [ -s tmp.txt ]; then echo "Errors from $f: " cat $f diff --git a/cmake/helpers/configure.cmake b/cmake/helpers/configure.cmake index 55b1df0f6d6e..7161215fbafa 100644 --- a/cmake/helpers/configure.cmake +++ b/cmake/helpers/configure.cmake @@ -44,6 +44,22 @@ check_type_size("long int" SIZEOF_LONG_INT) check_type_size("void*" SIZEOF_VOIDP) check_type_size("size_t" SIZEOF_SIZE_T) +# Check whether the type `_Float16` exists, and whether type +# conversions work (there might be linker problems) +check_cxx_source_compiles( + " + #ifdef GDAL_DISABLE_FLOAT16 + Explicitly disable _Float16 support + #endif + int main() { + _Float16 h = 1; + float f = h; + double d = h; + return f != 0 && d != 0 ? 0 : 1; + } + " + HAVE__FLOAT16) + check_function_exists(ctime_r HAVE_CTIME_R) check_function_exists(gmtime_r HAVE_GMTIME_R) check_function_exists(localtime_r HAVE_LOCALTIME_R) diff --git a/cmake/template/cpl_config.h.in b/cmake/template/cpl_config.h.in index cfd56cb60cb4..6baf7610e81a 100644 --- a/cmake/template/cpl_config.h.in +++ b/cmake/template/cpl_config.h.in @@ -24,6 +24,12 @@ /* The size of `size_t', as computed by sizeof. */ #cmakedefine SIZEOF_SIZE_T @SIZEOF_SIZE_T@ +/* Whether `std::float16_t` is available (and working). */ +#cmakedefine HAVE_STD_FLOAT16_T @HAVE_STD_FLOAT16_T@ + +/* Whether `_Float16' is supported. */ +#cmakedefine HAVE__FLOAT16 @HAVE__FLOAT16@ + /* Define to 1, if you have LARGEFILE64_SOURCE */ #cmakedefine VSI_NEED_LARGEFILE64_SOURCE 1 diff --git a/doc/source/development/dev_practices.rst b/doc/source/development/dev_practices.rst index 3121e138bce9..ed2fbb722327 100644 --- a/doc/source/development/dev_practices.rst +++ b/doc/source/development/dev_practices.rst @@ -56,6 +56,7 @@ used in GDAL/OGR. - *i*: integer number used as a zero based array or loop index. - *f*: floating point value (single precision) - *h*: an opaque handle (such as GDALDatasetH). +- *hf*: floating point value (half precision) - *n*: integer number (size unspecified) - *o*: C++ object - *os*: CPLString or std::string diff --git a/doc/source/user/raster_data_model.rst b/doc/source/user/raster_data_model.rst index fbc727ce0e1e..b2fbaaae5c6c 100644 --- a/doc/source/user/raster_data_model.rst +++ b/doc/source/user/raster_data_model.rst @@ -217,11 +217,13 @@ A raster band is represented in GDAL with the :cpp:class:`GDALRasterBand` class. A raster band has the following properties: - A width and height in pixels and lines. This is the same as that defined for the dataset, if this is a full resolution band. -- A datatype (GDALDataType). One of Byte, Int8, UInt16, Int16, UInt32, Int32, UInt64, Int64, Float32, Float64, and the complex types CInt16, CInt32, CFloat32, and CFloat64. +- A datatype (GDALDataType). One of Byte, Int8, UInt16, Int16, UInt32, Int32, UInt64, Int64, Float16, Float32, Float64, and the complex types CInt16, CInt32, CFloat16, CFloat32, and CFloat64. UInt64 and Int64 data types have been added in GDAL 3.5. Beyond reading and write pixel values, their support is limited. Some algorithms might use 64-bit floating-point internally (warping), as well as some methods returning only double values (GetMinimum(), GetMaximum(), etc.), or even 32-bit floating point (overview, RasterIO resampling). Hence the range where exact values are preserved can be [0, 2^53] (or less if 32-bit floating-point is used). - Int8 data type has been added in GDAL 3.7. + The Int8 data type has been added in GDAL 3.7. + + The Float16 and CFloat16 data types have been added in GDAL 3.11. If this data type is not supported by the hardware, then a software emulation is used. Not all drivers support Float16 yet. - A block size. This is a preferred (efficient) access chunk size. For tiled images this will be one tile. For scanline oriented images this will normally be one scanline. - A list of name/value pair metadata in the same format as the dataset, but of information that is potentially specific to this band. diff --git a/frmts/envisat/records.h b/frmts/envisat/records.h index a59e0db67fd9..7fb6dc8d4794 100644 --- a/frmts/envisat/records.h +++ b/frmts/envisat/records.h @@ -32,10 +32,12 @@ extern "C" /*! Sixteen bit signed integer */ EDT_Int16 = GDT_Int16, /*! Thirty two bit unsigned integer */ EDT_UInt32 = GDT_UInt32, /*! Thirty two bit signed integer */ EDT_Int32 = GDT_Int32, + /*! Sixteen bit floating point */ EDT_Float16 = GDT_Float16, /*! Thirty two bit floating point */ EDT_Float32 = GDT_Float32, /*! Sixty four bit floating point */ EDT_Float64 = GDT_Float64, /*! Complex Int16 */ EDT_CInt16 = GDT_CInt16, /*! Complex Int32 */ EDT_CInt32 = GDT_CInt32, + /*! Complex Float16 */ EDT_CFloat16 = GDT_CFloat16, /*! Complex Float32 */ EDT_CFloat32 = GDT_CFloat32, /*! Complex Float64 */ EDT_CFloat64 = GDT_CFloat64, /*! Modified Julian Dated */ EDT_MJD = GDT_TypeCount + 1, diff --git a/frmts/gti/data/gdaltileindex.xsd b/frmts/gti/data/gdaltileindex.xsd index 9289018c5f7e..90393065789f 100644 --- a/frmts/gti/data/gdaltileindex.xsd +++ b/frmts/gti/data/gdaltileindex.xsd @@ -134,10 +134,12 @@ + + diff --git a/frmts/gtiff/gt_overview.cpp b/frmts/gtiff/gt_overview.cpp index 001273ee9614..8d813861bcb4 100644 --- a/frmts/gtiff/gt_overview.cpp +++ b/frmts/gtiff/gt_overview.cpp @@ -368,6 +368,13 @@ CPLErr GTIFFBuildOverviewsEx(const char *pszFilename, int nBands, nBandFormat = SAMPLEFORMAT_INT; break; + case GDT_Float16: + // Convert Float16 to float. + // TODO: At some point we should support Float16. + nBandBits = 32; + nBandFormat = SAMPLEFORMAT_IEEEFP; + break; + case GDT_Float32: nBandBits = 32; nBandFormat = SAMPLEFORMAT_IEEEFP; @@ -388,6 +395,13 @@ CPLErr GTIFFBuildOverviewsEx(const char *pszFilename, int nBands, nBandFormat = SAMPLEFORMAT_COMPLEXINT; break; + case GDT_CFloat16: + // Convert Float16 to float. + // TODO: At some point we should support Float16. + nBandBits = 64; + nBandFormat = SAMPLEFORMAT_COMPLEXIEEEFP; + break; + case GDT_CFloat32: nBandBits = 64; nBandFormat = SAMPLEFORMAT_COMPLEXIEEEFP; diff --git a/frmts/gtiff/gtiffdataset_write.cpp b/frmts/gtiff/gtiffdataset_write.cpp index a1123263941e..1829db745a35 100644 --- a/frmts/gtiff/gtiffdataset_write.cpp +++ b/frmts/gtiff/gtiffdataset_write.cpp @@ -30,6 +30,7 @@ #include "cpl_error.h" #include "cpl_error_internal.h" // CPLErrorHandlerAccumulatorStruct +#include "cpl_float.h" #include "cpl_md5.h" #include "cpl_vsi.h" #include "cpl_vsi_virtual.h" diff --git a/frmts/hdf5/hdf5multidim.cpp b/frmts/hdf5/hdf5multidim.cpp index 9706d920e46e..6def9810a91a 100644 --- a/frmts/hdf5/hdf5multidim.cpp +++ b/frmts/hdf5/hdf5multidim.cpp @@ -2206,6 +2206,13 @@ GetHDF5DataTypeFromGDALDataType(const GDALExtendedDataType &dt, hid_t hNativeDT, case GDT_Int64: hBufferType = H5Tcopy(H5T_NATIVE_INT64); break; + case GDT_Float16: +#ifdef HDF5_HAVE_FLOAT16 + hBufferType = H5Tcopy(H5T_NATIVE_FLOAT16); + break; +#else + return H5I_INVALID_HID; +#endif case GDT_Float32: hBufferType = H5Tcopy(H5T_NATIVE_FLOAT); break; @@ -2214,6 +2221,7 @@ GetHDF5DataTypeFromGDALDataType(const GDALExtendedDataType &dt, hid_t hNativeDT, break; case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: if (bufferDataType != dt) diff --git a/frmts/netcdf/netcdfdataset.cpp b/frmts/netcdf/netcdfdataset.cpp index fa460ac4fd1c..89ab7f5e1704 100644 --- a/frmts/netcdf/netcdfdataset.cpp +++ b/frmts/netcdf/netcdfdataset.cpp @@ -47,6 +47,7 @@ #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_json.h" #include "cpl_minixml.h" #include "cpl_multiproc.h" diff --git a/frmts/tiledb/tiledbdense.cpp b/frmts/tiledb/tiledbdense.cpp index f3079e20085c..8111cf25f4d4 100644 --- a/frmts/tiledb/tiledbdense.cpp +++ b/frmts/tiledb/tiledbdense.cpp @@ -514,6 +514,11 @@ double TileDBRasterBand::GetNoDataValue(int *pbHasNoData) dfNoData = static_cast( *static_cast(value)); break; + case GDT_Float16: + case GDT_CFloat16: + // tileDB does not support float16 + CPLAssert(false); + break; case GDT_Float32: case GDT_CFloat32: dfNoData = *static_cast(value); @@ -600,9 +605,14 @@ CPLErr TileDBRasterBand::SetNoDataValue(double dfNoData) case GDT_Int64: bIsValid = IsValidNoData(dfNoData); break; + case GDT_Float16: + case GDT_CFloat16: + // tileDB does not support float16 + bIsValid = CPLIsNan(dfNoData) || IsValidNoData(dfNoData); + break; case GDT_Float32: case GDT_CFloat32: - bIsValid = std::isnan(dfNoData) || IsValidNoData(dfNoData); + bIsValid = CPLIsNan(dfNoData) || IsValidNoData(dfNoData); break; case GDT_Float64: case GDT_CFloat64: diff --git a/frmts/tiledb/tiledbmultidimarray.cpp b/frmts/tiledb/tiledbmultidimarray.cpp index 09afa8e868b9..fd8d69d3fce5 100644 --- a/frmts/tiledb/tiledbmultidimarray.cpp +++ b/frmts/tiledb/tiledbmultidimarray.cpp @@ -1281,6 +1281,11 @@ FillBlockSize(const std::vector> &aoDimensions, case GDT_Int64: tiledb_dt = TILEDB_INT64; break; + case GDT_CFloat16: + case GDT_Float16: + // tileDB does not support float16 + tiledb_dt = TILEDB_FLOAT32; + break; case GDT_CFloat32: case GDT_Float32: tiledb_dt = TILEDB_FLOAT32; diff --git a/frmts/tiledb/tiledbsparse.cpp b/frmts/tiledb/tiledbsparse.cpp index 9dfbc0a35fc3..4ec7384f6f3f 100644 --- a/frmts/tiledb/tiledbsparse.cpp +++ b/frmts/tiledb/tiledbsparse.cpp @@ -12,6 +12,7 @@ #include "tiledbheaders.h" +#include "cpl_float.h" #include "cpl_json.h" #include "cpl_time.h" #include "ogr_p.h" diff --git a/frmts/vrt/data/gdalvrt.xsd b/frmts/vrt/data/gdalvrt.xsd index 6cc250615091..a34a332608f2 100644 --- a/frmts/vrt/data/gdalvrt.xsd +++ b/frmts/vrt/data/gdalvrt.xsd @@ -425,16 +425,19 @@ + + + diff --git a/frmts/vrt/pixelfunctions.cpp b/frmts/vrt/pixelfunctions.cpp index da919aef1a0f..1e2412d85257 100644 --- a/frmts/vrt/pixelfunctions.cpp +++ b/frmts/vrt/pixelfunctions.cpp @@ -44,6 +44,8 @@ inline double GetSrcVal(const void *pSource, GDALDataType eSrcType, T ii) case GDT_Int64: return static_cast( static_cast(pSource)[ii]); + case GDT_Float16: + return static_cast(pSource)[ii]; case GDT_Float32: return static_cast(pSource)[ii]; case GDT_Float64: @@ -52,6 +54,8 @@ inline double GetSrcVal(const void *pSource, GDALDataType eSrcType, T ii) return static_cast(pSource)[2 * ii]; case GDT_CInt32: return static_cast(pSource)[2 * ii]; + case GDT_CFloat16: + return static_cast(pSource)[2 * ii]; case GDT_CFloat32: return static_cast(pSource)[2 * ii]; case GDT_CFloat64: diff --git a/frmts/vrt/vrtderivedrasterband.cpp b/frmts/vrt/vrtderivedrasterband.cpp index dfc7816c520c..3dd206191a07 100644 --- a/frmts/vrt/vrtderivedrasterband.cpp +++ b/frmts/vrt/vrtderivedrasterband.cpp @@ -89,6 +89,9 @@ static PyObject *GDALCreateNumpyArray(PyObject *pCreateArray, void *pBuffer, case GDT_UInt64: pszDataType = "uint64"; break; + case GDT_Float16: + pszDataType = "float16"; + break; case GDT_Float32: pszDataType = "float32"; break; @@ -99,6 +102,9 @@ static PyObject *GDALCreateNumpyArray(PyObject *pCreateArray, void *pBuffer, case GDT_CInt32: CPLAssert(FALSE); break; + case GDT_CFloat16: + CPLAssert(FALSE); + break; case GDT_CFloat32: pszDataType = "complex64"; break; @@ -1209,18 +1215,21 @@ CPLErr VRTDerivedRasterBand::IRasterIO( { eErr = CE_Failure; - // numpy doesn't have native cint16/cint32 - if (eSrcType == GDT_CInt16 || eSrcType == GDT_CInt32) + // numpy doesn't have native cint16/cint32/cfloat16 + if (eSrcType == GDT_CInt16 || eSrcType == GDT_CInt32 || + eSrcType == GDT_CFloat16) { - CPLError( - CE_Failure, CPLE_AppDefined, - "CInt16/CInt32 data type not supported for SourceTransferType"); + CPLError(CE_Failure, CPLE_AppDefined, + "CInt16/CInt32/CFloat16 data type not supported for " + "SourceTransferType"); goto end; } - if (eDataType == GDT_CInt16 || eDataType == GDT_CInt32) + if (eDataType == GDT_CInt16 || eDataType == GDT_CInt32 || + eDataType == GDT_CFloat16) { - CPLError(CE_Failure, CPLE_AppDefined, - "CInt16/CInt32 data type not supported for data type"); + CPLError( + CE_Failure, CPLE_AppDefined, + "CInt16/CInt32/CFloat16 data type not supported for data type"); goto end; } diff --git a/frmts/vrt/vrtdriver.cpp b/frmts/vrt/vrtdriver.cpp index b262b412ec40..36f922b03a26 100644 --- a/frmts/vrt/vrtdriver.cpp +++ b/frmts/vrt/vrtdriver.cpp @@ -512,8 +512,8 @@ void GDALRegister_VRT() poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Byte Int8 Int16 UInt16 Int32 UInt32 Int64 UInt64 " - "Float32 Float64 " - "CInt16 CInt32 CFloat32 CFloat64"); + "Float16 Float32 Float64 " + "CInt16 CInt32 CFloat16 CFloat32 CFloat64"); poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST, "\n" diff --git a/frmts/vrt/vrtprocesseddatasetfunctions.cpp b/frmts/vrt/vrtprocesseddatasetfunctions.cpp index 1f8465b979b2..a93fe48abcd1 100644 --- a/frmts/vrt/vrtprocesseddatasetfunctions.cpp +++ b/frmts/vrt/vrtprocesseddatasetfunctions.cpp @@ -39,6 +39,11 @@ static inline double GetDstValue(double dfVal, double dfDstNoData, { return dfReplacementDstNodata; } + else if (eIntendedDstDT == GDT_Float16 && + static_cast(dfVal) == static_cast(dfDstNoData)) + { + return dfReplacementDstNodata; + } else if (eIntendedDstDT == GDT_Float32 && static_cast(dfVal) == static_cast(dfDstNoData)) { diff --git a/frmts/vrt/vrtrasterband.cpp b/frmts/vrt/vrtrasterband.cpp index e7934238c6b0..60ed41307301 100644 --- a/frmts/vrt/vrtrasterband.cpp +++ b/frmts/vrt/vrtrasterband.cpp @@ -606,6 +606,16 @@ CPLString VRTSerializeNoData(double dfVal, GDALDataType eDataType, { return "nan"; } + else if (eDataType == GDT_Float16 && dfVal == -6.55e4) + { + // To avoid rounding out of the range of GFloat16 + return "-6.55e4"; + } + else if (eDataType == GDT_Float16 && dfVal == 6.55e4) + { + // To avoid rounding out of the range of GFloat16 + return "6.55e4"; + } else if (eDataType == GDT_Float32 && dfVal == -std::numeric_limits::max()) { @@ -885,7 +895,8 @@ bool VRTRasterBand::IsNoDataValueInDataTypeRange() const if (!m_bNoDataValueSet) return true; if (!std::isfinite(m_dfNoDataValue)) - return eDataType == GDT_Float32 || eDataType == GDT_Float64; + return eDataType == GDT_Float16 || eDataType == GDT_Float32 || + eDataType == GDT_Float64; GByte abyTempBuffer[2 * sizeof(double)]; CPLAssert(GDALGetDataTypeSizeBytes(eDataType) <= static_cast(sizeof(abyTempBuffer))); diff --git a/frmts/zarr/zarr_array.cpp b/frmts/zarr/zarr_array.cpp index f2e8ec04d352..586c6b1dfef3 100644 --- a/frmts/zarr/zarr_array.cpp +++ b/frmts/zarr/zarr_array.cpp @@ -1496,7 +1496,7 @@ bool ZarrArray::IRead(const GUInt64 *arrayStartIdx, const size_t *count, } /************************************************************************/ -/* ZarrArray::IRead() */ +/* ZarrArray::IWrite() */ /************************************************************************/ bool ZarrArray::IWrite(const GUInt64 *arrayStartIdx, const size_t *count, diff --git a/frmts/zarr/zarr_v2_array.cpp b/frmts/zarr/zarr_v2_array.cpp index ee5f1c161c4d..e62af1fe9b9c 100644 --- a/frmts/zarr/zarr_v2_array.cpp +++ b/frmts/zarr/zarr_v2_array.cpp @@ -1140,9 +1140,11 @@ static GDALExtendedDataType ParseDtype(const CPLJSONObject &obj, } else if (chType == 'f' && nBytes == 2) { + // elt.nativeType = DtypeElt::NativeType::IEEEFP; + // elt.gdalTypeIsApproxOfNative = true; + // eDT = GDT_Float32; elt.nativeType = DtypeElt::NativeType::IEEEFP; - elt.gdalTypeIsApproxOfNative = true; - eDT = GDT_Float32; + eDT = GDT_Float16; } else if (chType == 'f' && nBytes == 4) { @@ -1698,6 +1700,13 @@ ZarrV2Group::LoadArray(const std::string &osArrayName, CPLError(CE_Failure, CPLE_AppDefined, "Invalid fill_value"); return nullptr; } + if (oType.GetNumericDataType() == GDT_Float16) + { + const GFloat16 hfNoDataValue = + static_cast(dfNoDataValue); + abyNoData.resize(sizeof(hfNoDataValue)); + memcpy(&abyNoData[0], &hfNoDataValue, sizeof(hfNoDataValue)); + } if (oType.GetNumericDataType() == GDT_Float32) { const float fNoDataValue = static_cast(dfNoDataValue); diff --git a/frmts/zarr/zarr_v2_group.cpp b/frmts/zarr/zarr_v2_group.cpp index 5364cad43fd6..d3828a8a8491 100644 --- a/frmts/zarr/zarr_v2_group.cpp +++ b/frmts/zarr/zarr_v2_group.cpp @@ -735,6 +735,12 @@ static CPLJSONObject FillDTypeElts(const GDALExtendedDataType &oDataType, dtype.Set(dummy, " ZarrV2Group::CreateMDArray( case GDT_Int64: oFilter.Add("dtype", " ZarrV2Group::CreateMDArray( case GDT_CInt32: oFilter.Add("dtype", "(dfNoDataValue); + abyNoData.resize(sizeof(hfNoDataValue)); + memcpy(&abyNoData[0], &hfNoDataValue, sizeof(hfNoDataValue)); + } else if (oType.GetNumericDataType() == GDT_Float32) { const float fNoDataValue = static_cast(dfNoDataValue); diff --git a/frmts/zarr/zarr_v3_group.cpp b/frmts/zarr/zarr_v3_group.cpp index 132e8dd839be..4a340778a3be 100644 --- a/frmts/zarr/zarr_v3_group.cpp +++ b/frmts/zarr/zarr_v3_group.cpp @@ -417,6 +417,12 @@ static CPLJSONObject FillDTypeElts(const GDALExtendedDataType &oDataType, dtype.Set(dummy, "int64"); break; } + case GDT_Float16: + { + elt.nativeType = DtypeElt::NativeType::IEEEFP; + dtype.Set(dummy, "float16"); + break; + } case GDT_Float32: { elt.nativeType = DtypeElt::NativeType::IEEEFP; @@ -436,6 +442,12 @@ static CPLJSONObject FillDTypeElts(const GDALExtendedDataType &oDataType, bUnsupported = true; break; } + case GDT_CFloat16: + { + elt.nativeType = DtypeElt::NativeType::COMPLEX_IEEEFP; + dtype.Set(dummy, "complex32"); + break; + } case GDT_CFloat32: { elt.nativeType = DtypeElt::NativeType::COMPLEX_IEEEFP; @@ -450,8 +462,8 @@ static CPLJSONObject FillDTypeElts(const GDALExtendedDataType &oDataType, } case GDT_TypeCount: { - static_assert(GDT_TypeCount == GDT_Int8 + 1, - "GDT_TypeCount == GDT_Int8 + 1"); + static_assert(GDT_TypeCount == GDT_CFloat16 + 1, + "GDT_TypeCount == GDT_CFloat16 + 1"); break; } } diff --git a/frmts/zarr/zarrdrivercore.cpp b/frmts/zarr/zarrdrivercore.cpp index f9f6b896240c..17f6ff4e01e5 100644 --- a/frmts/zarr/zarrdrivercore.cpp +++ b/frmts/zarr/zarrdrivercore.cpp @@ -68,9 +68,10 @@ void ZARRDriverSetCommonMetadata(GDALDriver *poDriver) poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES"); poDriver->SetMetadataItem(GDAL_DCAP_MULTIDIM_RASTER, "YES"); poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Zarr"); - poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, - "Byte Int16 UInt16 Int32 UInt32 Int64 UInt64 " - "Float32 Float64 CFloat32 CFloat64"); + poDriver->SetMetadataItem( + GDAL_DMD_CREATIONDATATYPES, + "Int8 Byte Int16 UInt16 Int32 UInt32 Int64 UInt64 " + "Float16 Float32 Float64 CFLoat16 CFloat32 CFloat64"); poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES"); poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES"); diff --git a/fuzzers/build.sh b/fuzzers/build.sh index e60e29195217..6b042d10bea2 100755 --- a/fuzzers/build.sh +++ b/fuzzers/build.sh @@ -89,7 +89,7 @@ if test "${CIFUZZ:-}" = "True"; then -L$SRC_DIR/build -lgdal \ -lproj \ -Wl,-Bstatic -lzstd -lwebp -llzma -lexpat -lsqlite3 -lgif -ljpeg -lz \ - -Wl,-Bdynamic -ldl -lpthread + -Wl,-Bdynamic -ldl -lpthread -lclang_rt.builtins echo "Building ogr_fuzzer" $CXX $CXXFLAGS \ @@ -103,7 +103,7 @@ if test "${CIFUZZ:-}" = "True"; then -L$SRC_DIR/build -lgdal \ -L$SRC/install/lib -lproj \ -Wl,-Bstatic -lzstd -lwebp -llzma -lexpat -lsqlite3 -lgif -ljpeg -lz \ - -Wl,-Bdynamic -ldl -lpthread + -Wl,-Bdynamic -ldl -lpthread -lclang_rt.builtins echo "Building gdal_fuzzer_seed_corpus.zip" cd $(dirname $0)/../autotest/gcore/data diff --git a/gcore/gdal.h b/gcore/gdal.h index 02adb37c43f6..0c5d9e01c86b 100644 --- a/gcore/gdal.h +++ b/gcore/gdal.h @@ -55,14 +55,16 @@ typedef enum /*! Thirty two bit signed integer */ GDT_Int32 = 5, /*! 64 bit unsigned integer (GDAL >= 3.5)*/ GDT_UInt64 = 12, /*! 64 bit signed integer (GDAL >= 3.5)*/ GDT_Int64 = 13, + /*! Sixteen bit floating point */ GDT_Float16 = 15, /*! Thirty two bit floating point */ GDT_Float32 = 6, /*! Sixty four bit floating point */ GDT_Float64 = 7, /*! Complex Int16 */ GDT_CInt16 = 8, /*! Complex Int32 */ GDT_CInt32 = 9, /* TODO?(#6879): GDT_CInt64 */ + /*! Complex Float16 */ GDT_CFloat16 = 16, /*! Complex Float32 */ GDT_CFloat32 = 10, /*! Complex Float64 */ GDT_CFloat64 = 11, - GDT_TypeCount = 15 /* maximum type # + 1 */ + GDT_TypeCount = 17 /* maximum type # + 1 */ } GDALDataType; int CPL_DLL CPL_STDCALL GDALGetDataTypeSize(GDALDataType); // Deprecated. diff --git a/gcore/gdal_misc.cpp b/gcore/gdal_misc.cpp index 09e65c2ca824..98353426cebe 100644 --- a/gcore/gdal_misc.cpp +++ b/gcore/gdal_misc.cpp @@ -30,6 +30,7 @@ #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_json.h" #include "cpl_minixml.h" #include "cpl_multiproc.h" @@ -87,7 +88,9 @@ static int GetDataTypeElementSizeBits(GDALDataType eDataType) case GDT_UInt16: case GDT_Int16: + case GDT_Float16: case GDT_CInt16: + case GDT_CFloat16: return 16; case GDT_UInt32: @@ -224,10 +227,12 @@ GDALDataType CPL_STDCALL GDALDataTypeUnionWithValue(GDALDataType eDT, return eDT; break; } + // Do not return `GDT_Float16` because that type is not supported everywhere + case GDT_Float16: case GDT_Float32: { if (GDALIsValueExactAs(dfValue)) - return eDT; + return GDT_Float32; break; } case GDT_Float64: @@ -237,6 +242,7 @@ GDALDataType CPL_STDCALL GDALDataTypeUnionWithValue(GDALDataType eDT, case GDT_Unknown: case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_TypeCount: @@ -255,34 +261,34 @@ static int GetMinBitsForValue(double dValue) { if (round(dValue) == dValue) { - if (dValue <= std::numeric_limits::max() && - dValue >= std::numeric_limits::min()) + if (dValue <= cpl::NumericLimits::max() && + dValue >= cpl::NumericLimits::lowest()) return 8; - if (dValue <= std::numeric_limits::max() && - dValue >= std::numeric_limits::min()) + if (dValue <= cpl::NumericLimits::max() && + dValue >= cpl::NumericLimits::lowest()) return 8; - if (dValue <= std::numeric_limits::max() && - dValue >= std::numeric_limits::min()) + if (dValue <= cpl::NumericLimits::max() && + dValue >= cpl::NumericLimits::lowest()) return 16; - if (dValue <= std::numeric_limits::max() && - dValue >= std::numeric_limits::min()) + if (dValue <= cpl::NumericLimits::max() && + dValue >= cpl::NumericLimits::lowest()) return 16; - if (dValue <= std::numeric_limits::max() && - dValue >= std::numeric_limits::min()) + if (dValue <= cpl::NumericLimits::max() && + dValue >= cpl::NumericLimits::lowest()) return 32; - if (dValue <= std::numeric_limits::max() && - dValue >= std::numeric_limits::min()) + if (dValue <= cpl::NumericLimits::max() && + dValue >= cpl::NumericLimits::lowest()) return 32; - if (dValue <= static_cast( - std::numeric_limits::max()) && - dValue >= - static_cast(std::numeric_limits::min())) + if (dValue <= + static_cast(cpl::NumericLimits::max()) && + dValue >= static_cast( + cpl::NumericLimits::lowest())) return 64; } else if (static_cast(dValue) == dValue) @@ -312,53 +318,80 @@ static int GetMinBitsForValue(double dValue) GDALDataType CPL_STDCALL GDALFindDataType(int nBits, int bSigned, int bFloating, int bComplex) { - if (bComplex) - { - nBits = std::max(nBits, !bSigned ? 32 : 16); - } // we don't have complex unsigned data types, so for a complex uint16, - // promote to complex int32 - if (bFloating) - { - nBits = std::max(nBits, 32); - } - - if (nBits <= 8) + if (!bFloating) { - return bSigned ? GDT_Int8 : GDT_Byte; - } - - if (nBits <= 16) - { - if (bComplex) - return GDT_CInt16; - if (bSigned) - return GDT_Int16; - return GDT_UInt16; + if (!bComplex) + { + if (!bSigned) + { + if (nBits <= 8) + return GDT_Byte; + if (nBits <= 16) + return GDT_UInt16; + if (nBits <= 32) + return GDT_UInt32; + if (nBits <= 64) + return GDT_UInt64; + return GDT_Float64; + } + else // bSigned + { + if (nBits <= 8) + return GDT_Int8; + if (nBits <= 16) + return GDT_Int16; + if (nBits <= 32) + return GDT_Int32; + if (nBits <= 64) + return GDT_Int64; + return GDT_Float64; + } + } + else // bComplex + { + if (!bSigned) + { + // We don't have complex unsigned data types, so + // return a large-enough complex signed type + + // Do not choose CInt16 for backward compatibility + // if (nBits <= 15) + // return GDT_CInt16; + if (nBits <= 31) + return GDT_CInt32; + return GDT_CFloat64; + } + else // bSigned + { + if (nBits <= 16) + return GDT_CInt16; + if (nBits <= 32) + return GDT_CInt32; + return GDT_CFloat64; + } + } } - - if (nBits <= 32) + else // bFloating { - if (bFloating) + if (!bComplex) + { + // Do not choose Float16 since is not supported everywhere + // if (nBits <= 16) + // return GDT_Float16; + if (nBits <= 32) + return GDT_Float32; + return GDT_Float64; + } + else // bComplex { - if (bComplex) + // Do not choose Float16 since is not supported everywhere + // if (nBits <= 16) + // return GDT_CFloat16; + if (nBits <= 32) return GDT_CFloat32; - return GDT_Float32; + return GDT_CFloat64; } - - if (bComplex) - return GDT_CInt32; - if (bSigned) - return GDT_Int32; - return GDT_UInt32; } - - if (nBits == 64 && !bFloating && !bComplex) - return bSigned ? GDT_Int64 : GDT_UInt64; - - if (bComplex) - return GDT_CFloat64; - - return GDT_Float64; } /************************************************************************/ @@ -379,9 +412,9 @@ GDALDataType CPL_STDCALL GDALFindDataTypeForValue(double dValue, int bComplex) const bool bFloating = round(dValue) != dValue || dValue > - static_cast(std::numeric_limits::max()) || + static_cast(cpl::NumericLimits::max()) || dValue < - static_cast(std::numeric_limits::lowest()); + static_cast(cpl::NumericLimits::lowest()); const bool bSigned = bFloating || dValue < 0; const int nBits = GetMinBitsForValue(dValue); @@ -413,12 +446,14 @@ int CPL_STDCALL GDALGetDataTypeSizeBytes(GDALDataType eDataType) case GDT_UInt16: case GDT_Int16: + case GDT_Float16: return 2; case GDT_UInt32: case GDT_Int32: case GDT_Float32: case GDT_CInt16: + case GDT_CFloat16: return 4; case GDT_Float64: @@ -499,6 +534,7 @@ int CPL_STDCALL GDALDataTypeIsComplex(GDALDataType eDataType) { case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: return TRUE; @@ -511,6 +547,7 @@ int CPL_STDCALL GDALDataTypeIsComplex(GDALDataType eDataType) case GDT_UInt32: case GDT_Int64: case GDT_UInt64: + case GDT_Float16: case GDT_Float32: case GDT_Float64: return FALSE; @@ -529,8 +566,8 @@ int CPL_STDCALL GDALDataTypeIsComplex(GDALDataType eDataType) /** * \brief Is data type floating? (might be complex) * - * @return TRUE if the passed type is floating (one of GDT_Float32, GDT_Float64, - * GDT_CFloat32, GDT_CFloat64) + * @return TRUE if the passed type is floating (one of GDT_Float32, GDT_Float16, + * GDT_Float64, GDT_CFloat16, GDT_CFloat32, GDT_CFloat64) * @since GDAL 2.3 */ @@ -538,8 +575,10 @@ int CPL_STDCALL GDALDataTypeIsFloating(GDALDataType eDataType) { switch (eDataType) { + case GDT_Float16: case GDT_Float32: case GDT_Float64: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: return TRUE; @@ -592,8 +631,10 @@ int CPL_STDCALL GDALDataTypeIsInteger(GDALDataType eDataType) case GDT_Int64: return TRUE; + case GDT_Float16: case GDT_Float32: case GDT_Float64: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: return FALSE; @@ -630,10 +671,12 @@ int CPL_STDCALL GDALDataTypeIsSigned(GDALDataType eDataType) case GDT_Int16: case GDT_Int32: case GDT_Int64: + case GDT_Float16: case GDT_Float32: case GDT_Float64: case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: return TRUE; @@ -693,6 +736,15 @@ int CPL_STDCALL GDALDataTypeIsConversionLossy(GDALDataType eTypeFrom, return FALSE; } + if (eTypeTo == GDT_Float16 && + (eTypeFrom == GDT_Int16 || eTypeFrom == GDT_UInt16 || + eTypeFrom == GDT_Int32 || eTypeFrom == GDT_UInt32 || + eTypeFrom == GDT_Int64 || eTypeFrom == GDT_UInt64 || + eTypeFrom == GDT_Float32 || eTypeFrom == GDT_Float64)) + { + return TRUE; + } + if (eTypeTo == GDT_Float32 && (eTypeFrom == GDT_Int32 || eTypeFrom == GDT_UInt32 || eTypeFrom == GDT_Int64 || eTypeFrom == GDT_UInt64 || @@ -760,6 +812,9 @@ const char *CPL_STDCALL GDALGetDataTypeName(GDALDataType eDataType) case GDT_Int64: return "Int64"; + case GDT_Float16: + return "Float16"; + case GDT_Float32: return "Float32"; @@ -772,6 +827,9 @@ const char *CPL_STDCALL GDALGetDataTypeName(GDALDataType eDataType) case GDT_CInt32: return "CInt32"; + case GDT_CFloat16: + return "CFloat16"; + case GDT_CFloat32: return "CFloat32"; @@ -825,17 +883,15 @@ template static inline void ClampAndRound(double &dfValue, bool &bClamped, bool &bRounded) { - // TODO(schwehr): Rework this template. ::min() versus ::lowest. - - if (dfValue < static_cast(std::numeric_limits::min())) + if (dfValue < static_cast(cpl::NumericLimits::lowest())) { bClamped = true; - dfValue = static_cast(std::numeric_limits::min()); + dfValue = static_cast(cpl::NumericLimits::lowest()); } - else if (dfValue > static_cast(std::numeric_limits::max())) + else if (dfValue > static_cast(cpl::NumericLimits::max())) { bClamped = true; - dfValue = static_cast(std::numeric_limits::max()); + dfValue = static_cast(cpl::NumericLimits::max()); } else if (dfValue != static_cast(static_cast(dfValue))) { @@ -892,28 +948,53 @@ double GDALAdjustValueToDataType(GDALDataType eDT, double dfValue, case GDT_UInt64: ClampAndRound(dfValue, bClamped, bRounded); break; - case GDT_Float32: + case GDT_Float16: { if (!std::isfinite(dfValue)) break; - // TODO(schwehr): ::min() versus ::lowest. - // Use ClampAndRound after it has been fixed. - if (dfValue < -std::numeric_limits::max()) + // TODO: Use ClampAndRound + if (dfValue < cpl::NumericLimits::lowest()) { bClamped = TRUE; dfValue = - static_cast(-std::numeric_limits::max()); + static_cast(cpl::NumericLimits::lowest()); } - else if (dfValue > std::numeric_limits::max()) + else if (dfValue > cpl::NumericLimits::max()) { bClamped = TRUE; dfValue = - static_cast(std::numeric_limits::max()); + static_cast(cpl::NumericLimits::max()); } else { - // Intentionally loose precision. + // Intentionally lose precision. + // TODO(schwehr): Is the double cast really necessary? + // If so, why? What will fail? + dfValue = static_cast(static_cast(dfValue)); + } + break; + } + case GDT_Float32: + { + if (!std::isfinite(dfValue)) + break; + + // TODO: Use ClampAndRound + if (dfValue < cpl::NumericLimits::lowest()) + { + bClamped = TRUE; + dfValue = + static_cast(cpl::NumericLimits::lowest()); + } + else if (dfValue > cpl::NumericLimits::max()) + { + bClamped = TRUE; + dfValue = static_cast(cpl::NumericLimits::max()); + } + else + { + // Intentionally lose precision. // TODO(schwehr): Is the double cast really necessary? // If so, why? What will fail? dfValue = static_cast(static_cast(dfValue)); @@ -923,6 +1004,7 @@ double GDALAdjustValueToDataType(GDALDataType eDT, double dfValue, case GDT_Float64: case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_Unknown: @@ -973,6 +1055,8 @@ bool GDALIsValueExactAs(double dfValue, GDALDataType eDT) return GDALIsValueExactAs(dfValue); case GDT_Int64: return GDALIsValueExactAs(dfValue); + case GDT_Float16: + return GDALIsValueExactAs(dfValue); case GDT_Float32: return GDALIsValueExactAs(dfValue); case GDT_Float64: @@ -980,6 +1064,7 @@ bool GDALIsValueExactAs(double dfValue, GDALDataType eDT) case GDT_Unknown: case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_TypeCount: @@ -1011,6 +1096,8 @@ GDALDataType CPL_STDCALL GDALGetNonComplexDataType(GDALDataType eDataType) return GDT_Int16; case GDT_CInt32: return GDT_Int32; + case GDT_CFloat16: + return GDT_Float16; case GDT_CFloat32: return GDT_Float32; case GDT_CFloat64: @@ -1024,6 +1111,7 @@ GDALDataType CPL_STDCALL GDALGetNonComplexDataType(GDALDataType eDataType) case GDT_Int16: case GDT_Int32: case GDT_Int64: + case GDT_Float16: case GDT_Float32: case GDT_Float64: break; @@ -1549,6 +1637,10 @@ int CPL_STDCALL GDALGetRandomRasterSample(GDALRasterBandH hBand, int nSamples, reinterpret_cast( pDataRef)[iOffset]); break; + case GDT_Float16: + dfValue = reinterpret_cast( + pDataRef)[iOffset]; + break; case GDT_Float32: dfValue = reinterpret_cast(pDataRef)[iOffset]; @@ -1576,6 +1668,17 @@ int CPL_STDCALL GDALGetRandomRasterSample(GDALRasterBandH hBand, int nSamples, dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); break; } + case GDT_CFloat16: + { + const double dfReal = + reinterpret_cast( + pDataRef)[iOffset * 2]; + const double dfImag = + reinterpret_cast( + pDataRef)[iOffset * 2 + 1]; + dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; + } case GDT_CFloat32: { const double dfReal = reinterpret_cast( @@ -4223,7 +4326,7 @@ int CPL_STDCALL GDALGeneralCmdLineProcessor(int nArgc, char ***ppapszArgv, { std::cout << "Hit to Continue." << std::endl; std::cin.clear(); - std::cin.ignore(std::numeric_limits::max(), '\n'); + std::cin.ignore(cpl::NumericLimits::max(), '\n'); } /* -------------------------------------------------------------------- @@ -5057,7 +5160,7 @@ bool GDALCanReliablyUseSiblingFileList(const char *pszFilename) double GDALAdjustNoDataCloseToFloatMax(double dfVal) { - const auto kMaxFloat = std::numeric_limits::max(); + const auto kMaxFloat = cpl::NumericLimits::max(); if (std::fabs(dfVal - -kMaxFloat) < 1e-10 * kMaxFloat) return -kMaxFloat; if (std::fabs(dfVal - kMaxFloat) < 1e-10 * kMaxFloat) @@ -5109,7 +5212,7 @@ void GDALCopyNoDataValue(GDALRasterBand *poDstBand, GDALRasterBand *poSrcBand) else if (eDstDataType == GDT_Int64) { if (nNoData < - static_cast(std::numeric_limits::max())) + static_cast(cpl::NumericLimits::max())) { poDstBand->SetNoDataValueAsInt64( static_cast(nNoData)); @@ -5130,9 +5233,9 @@ void GDALCopyNoDataValue(GDALRasterBand *poDstBand, GDALRasterBand *poSrcBand) if (eDstDataType == GDT_Int64) { if (dfNoData >= static_cast( - std::numeric_limits::min()) && + cpl::NumericLimits::lowest()) && dfNoData <= static_cast( - std::numeric_limits::max()) && + cpl::NumericLimits::max()) && dfNoData == static_cast(static_cast(dfNoData))) { @@ -5143,9 +5246,9 @@ void GDALCopyNoDataValue(GDALRasterBand *poDstBand, GDALRasterBand *poSrcBand) else if (eDstDataType == GDT_UInt64) { if (dfNoData >= static_cast( - std::numeric_limits::min()) && + cpl::NumericLimits::lowest()) && dfNoData <= static_cast( - std::numeric_limits::max()) && + cpl::NumericLimits::max()) && dfNoData == static_cast(static_cast(dfNoData))) { @@ -5430,78 +5533,78 @@ double GDALGetNoDataReplacementValue(GDALDataType dt, double dfNoDataValue) if (dt == GDT_Byte) { if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } - if (dfNoDataValue == std::numeric_limits::max()) - dfReplacementVal = std::numeric_limits::max() - 1; + if (dfNoDataValue == cpl::NumericLimits::max()) + dfReplacementVal = cpl::NumericLimits::max() - 1; else dfReplacementVal = dfNoDataValue + 1; } else if (dt == GDT_Int8) { if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } - if (dfNoDataValue == std::numeric_limits::max()) - dfReplacementVal = std::numeric_limits::max() - 1; + if (dfNoDataValue == cpl::NumericLimits::max()) + dfReplacementVal = cpl::NumericLimits::max() - 1; else dfReplacementVal = dfNoDataValue + 1; } else if (dt == GDT_UInt16) { if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } - if (dfNoDataValue == std::numeric_limits::max()) - dfReplacementVal = std::numeric_limits::max() - 1; + if (dfNoDataValue == cpl::NumericLimits::max()) + dfReplacementVal = cpl::NumericLimits::max() - 1; else dfReplacementVal = dfNoDataValue + 1; } else if (dt == GDT_Int16) { if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } - if (dfNoDataValue == std::numeric_limits::max()) - dfReplacementVal = std::numeric_limits::max() - 1; + if (dfNoDataValue == cpl::NumericLimits::max()) + dfReplacementVal = cpl::NumericLimits::max() - 1; else dfReplacementVal = dfNoDataValue + 1; } else if (dt == GDT_UInt32) { if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } - if (dfNoDataValue == std::numeric_limits::max()) - dfReplacementVal = std::numeric_limits::max() - 1; + if (dfNoDataValue == cpl::NumericLimits::max()) + dfReplacementVal = cpl::NumericLimits::max() - 1; else dfReplacementVal = dfNoDataValue + 1; } else if (dt == GDT_Int32) { if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } - if (dfNoDataValue == std::numeric_limits::max()) - dfReplacementVal = std::numeric_limits::max() - 1; + if (dfNoDataValue == cpl::NumericLimits::max()) + dfReplacementVal = cpl::NumericLimits::max() - 1; else dfReplacementVal = dfNoDataValue + 1; } @@ -5511,18 +5614,18 @@ double GDALGetNoDataReplacementValue(GDALDataType dt, double dfNoDataValue) // so we take the next lower value representable as a double 18446744073709549567 static const double dfMaxUInt64Value{ std::nextafter( - static_cast(std::numeric_limits::max()), 0) - + static_cast(cpl::NumericLimits::max()), 0) - 1}; if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } if (dfNoDataValue >= - static_cast(std::numeric_limits::max())) + static_cast(cpl::NumericLimits::max())) dfReplacementVal = dfMaxUInt64Value; else dfReplacementVal = dfNoDataValue + 1; @@ -5533,61 +5636,83 @@ double GDALGetNoDataReplacementValue(GDALDataType dt, double dfNoDataValue) // so we take the next lower value representable as a double 9223372036854774784 static const double dfMaxInt64Value{ std::nextafter( - static_cast(std::numeric_limits::max()), 0) - + static_cast(cpl::NumericLimits::max()), 0) - 1}; if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } if (dfNoDataValue >= - static_cast(std::numeric_limits::max())) + static_cast(cpl::NumericLimits::max())) dfReplacementVal = dfMaxInt64Value; else dfReplacementVal = dfNoDataValue + 1; } + else if (dt == GDT_Float16) + { + + if (GDALClampDoubleValue(dfNoDataValue, + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) + { + return 0; + } + + if (dfNoDataValue == cpl::NumericLimits::max()) + { + using std::nextafter; + dfReplacementVal = + nextafter(static_cast(dfNoDataValue), GFloat16(0.0f)); + } + else + { + using std::nextafter; + dfReplacementVal = nextafter(static_cast(dfNoDataValue), + cpl::NumericLimits::max()); + } + } else if (dt == GDT_Float32) { if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } - if (dfNoDataValue == std::numeric_limits::max()) + if (dfNoDataValue == cpl::NumericLimits::max()) { dfReplacementVal = std::nextafter(static_cast(dfNoDataValue), 0.0f); } else { - dfReplacementVal = - std::nextafter(static_cast(dfNoDataValue), - std::numeric_limits::max()); + dfReplacementVal = std::nextafter(static_cast(dfNoDataValue), + cpl::NumericLimits::max()); } } else if (dt == GDT_Float64) { if (GDALClampDoubleValue(dfNoDataValue, - std::numeric_limits::lowest(), - std::numeric_limits::max())) + cpl::NumericLimits::lowest(), + cpl::NumericLimits::max())) { return 0; } - if (dfNoDataValue == std::numeric_limits::max()) + if (dfNoDataValue == cpl::NumericLimits::max()) { - dfReplacementVal = std::nextafter(dfNoDataValue, 0.0f); + dfReplacementVal = std::nextafter(dfNoDataValue, 0.0); } else { dfReplacementVal = std::nextafter( - dfNoDataValue, std::numeric_limits::max()); + dfNoDataValue, cpl::NumericLimits::max()); } } diff --git a/gcore/gdal_priv.h b/gcore/gdal_priv.h index 5d6c653ac131..2ae5bf76f764 100644 --- a/gcore/gdal_priv.h +++ b/gcore/gdal_priv.h @@ -46,6 +46,7 @@ class GDALRelationship; #include "gdalsubdatasetinfo.h" #include "cpl_vsi.h" #include "cpl_conv.h" +#include "cpl_float.h" #include "cpl_string.h" #include "cpl_minixml.h" #include "cpl_multiproc.h" @@ -3695,10 +3696,10 @@ class CPL_DLL GDALMDArray : virtual public GDALAbstractMDArray, Transpose(const std::vector &anMapNewAxisToOldAxis) const; std::shared_ptr GetUnscaled( - double dfOverriddenScale = std::numeric_limits::quiet_NaN(), - double dfOverriddenOffset = std::numeric_limits::quiet_NaN(), + double dfOverriddenScale = cpl::NumericLimits::quiet_NaN(), + double dfOverriddenOffset = cpl::NumericLimits::quiet_NaN(), double dfOverriddenDstNodata = - std::numeric_limits::quiet_NaN()) const; + cpl::NumericLimits::quiet_NaN()) const; virtual std::shared_ptr GetMask(CSLConstList papszOptions) const; @@ -4570,11 +4571,34 @@ GDALDataset *GDALCreateOverviewDataset(GDALDataset *poDS, int nOvrLevel, // Should cover particular cases of #3573, #4183, #4506, #6578 // Behavior is undefined if fVal1 or fVal2 are NaN (should be tested before // calling this function) -template inline bool ARE_REAL_EQUAL(T fVal1, T fVal2, int ulp = 2) + +// TODO: The expression `abs(fVal1 + fVal2)` looks strange; is this a bug? +// Should this be `abs(fVal1) + abs(fVal2)` instead? + +inline bool ARE_REAL_EQUAL(GFloat16 dfVal1, GFloat16 dfVal2, int ulp = 2) +{ + using std::abs; + return dfVal1 == dfVal2 || /* Should cover infinity */ + abs(dfVal1 - dfVal2) < cpl::NumericLimits::epsilon() * + abs(dfVal1 + dfVal2) * ulp; +} + +inline bool ARE_REAL_EQUAL(float fVal1, float fVal2, int ulp = 2) { + using std::abs; return fVal1 == fVal2 || /* Should cover infinity */ - std::abs(fVal1 - fVal2) < std::numeric_limits::epsilon() * - std::abs(fVal1 + fVal2) * ulp; + abs(fVal1 - fVal2) < + cpl::NumericLimits::epsilon() * abs(fVal1 + fVal2) * ulp; +} + +// We are using `cpl::NumericLimits::epsilon()` for backward +// compatibility +inline bool ARE_REAL_EQUAL(double dfVal1, double dfVal2, int ulp = 2) +{ + using std::abs; + return dfVal1 == dfVal2 || /* Should cover infinity */ + abs(dfVal1 - dfVal2) < cpl::NumericLimits::epsilon() * + abs(dfVal1 + dfVal2) * ulp; } double GDALAdjustNoDataCloseToFloatMax(double dfVal); diff --git a/gcore/gdal_priv_templates.hpp b/gcore/gdal_priv_templates.hpp index 0622ee2a76c9..410e315a2f02 100644 --- a/gcore/gdal_priv_templates.hpp +++ b/gcore/gdal_priv_templates.hpp @@ -22,6 +22,8 @@ #include #include +#include "cpl_float.h" + /************************************************************************/ /* GDALGetDataLimits() */ /************************************************************************/ @@ -37,37 +39,39 @@ template inline void GDALGetDataLimits(Tin &tMaxValue, Tin &tMinValue) { - tMaxValue = std::numeric_limits::max(); - tMinValue = std::numeric_limits::min(); + tMaxValue = cpl::NumericLimits::max(); + tMinValue = cpl::NumericLimits::lowest(); // Compute the actual minimum value of Tout in terms of Tin. - if constexpr (std::numeric_limits::is_signed && - std::numeric_limits::is_integer) + if constexpr (cpl::NumericLimits::is_signed && + cpl::NumericLimits::is_integer) { // the minimum value is less than zero - if constexpr (std::numeric_limits::digits < - std::numeric_limits::digits || - !std::numeric_limits::is_integer) + // cppcheck-suppress knownConditionTrueFalse + if constexpr (cpl::NumericLimits::digits < + cpl::NumericLimits::digits || + !cpl::NumericLimits::is_integer) { // Tout is smaller than Tin, so we need to clamp values in input // to the range of Tout's min/max values - if (std::numeric_limits::is_signed) + if constexpr (cpl::NumericLimits::is_signed) { - tMinValue = static_cast(std::numeric_limits::min()); + tMinValue = + static_cast(cpl::NumericLimits::lowest()); } - tMaxValue = static_cast(std::numeric_limits::max()); + tMaxValue = static_cast(cpl::NumericLimits::max()); } } - else if constexpr (std::numeric_limits::is_integer) + else if constexpr (cpl::NumericLimits::is_integer) { // the output is unsigned, so we just need to determine the max /* coverity[same_on_both_sides] */ - if constexpr (std::numeric_limits::digits <= - std::numeric_limits::digits) + if constexpr (cpl::NumericLimits::digits <= + cpl::NumericLimits::digits) { // Tout is smaller than Tin, so we need to clamp the input values // to the range of Tout's max - tMaxValue = static_cast(std::numeric_limits::max()); + tMaxValue = static_cast(cpl::NumericLimits::max()); } tMinValue = 0; } @@ -129,30 +133,37 @@ inline bool GDALClampDoubleValue(double &tValue, const T2 tMin, const T3 tMax) */ template inline bool GDALIsValueInRange(double dfValue) { - return dfValue >= static_cast(std::numeric_limits::lowest()) && - dfValue <= static_cast(std::numeric_limits::max()); + return dfValue >= static_cast(cpl::NumericLimits::lowest()) && + dfValue <= static_cast(cpl::NumericLimits::max()); } template <> inline bool GDALIsValueInRange(double dfValue) { - return !std::isnan(dfValue); + return !CPLIsNan(dfValue); } template <> inline bool GDALIsValueInRange(double dfValue) { - return std::isinf(dfValue) || + return CPLIsInf(dfValue) || (dfValue >= -std::numeric_limits::max() && dfValue <= std::numeric_limits::max()); } +template <> inline bool GDALIsValueInRange(double dfValue) +{ + return CPLIsInf(dfValue) || + (dfValue >= -cpl::NumericLimits::max() && + dfValue <= cpl::NumericLimits::max()); +} + template <> inline bool GDALIsValueInRange(double dfValue) { // Values in the range [INT64_MAX - 1023, INT64_MAX - 1] // get converted to a double that once cast to int64_t is // INT64_MAX + 1, hence the < strict comparison. return dfValue >= - static_cast(std::numeric_limits::min()) && - dfValue < static_cast(std::numeric_limits::max()); + static_cast(cpl::NumericLimits::lowest()) && + dfValue < static_cast(cpl::NumericLimits::max()); } template <> inline bool GDALIsValueInRange(double dfValue) @@ -161,7 +172,7 @@ template <> inline bool GDALIsValueInRange(double dfValue) // get converted to a double that once cast to uint64_t is // UINT64_MAX + 1, hence the < strict comparison. return dfValue >= 0 && - dfValue < static_cast(std::numeric_limits::max()); + dfValue < static_cast(cpl::NumericLimits::max()); } /************************************************************************/ @@ -186,11 +197,18 @@ template inline bool GDALIsValueExactAs(double dfValue) template <> inline bool GDALIsValueExactAs(double dfValue) { - return std::isnan(dfValue) || + return CPLIsNan(dfValue) || (GDALIsValueInRange(dfValue) && static_cast(static_cast(dfValue)) == dfValue); } +template <> inline bool GDALIsValueExactAs(double dfValue) +{ + return CPLIsNan(dfValue) || + (GDALIsValueInRange(dfValue) && + static_cast(static_cast(dfValue)) == dfValue); +} + template <> inline bool GDALIsValueExactAs(double) { return true; @@ -200,6 +218,8 @@ template <> inline bool GDALIsValueExactAs(double) /* GDALCopyWord() */ /************************************************************************/ +// Integer input and output: clamp the input + template struct sGDALCopyWord { static inline void f(const Tin tValueIn, Tout &tValueOut) @@ -211,6 +231,16 @@ template struct sGDALCopyWord } }; +// Integer input and floating point output: simply convert + +template struct sGDALCopyWord +{ + static inline void f(const Tin tValueIn, GFloat16 &hfValueOut) + { + hfValueOut = static_cast(tValueIn); + } +}; + template struct sGDALCopyWord { static inline void f(const Tin tValueIn, float &fValueOut) @@ -227,6 +257,50 @@ template struct sGDALCopyWord } }; +// Floating point input and output, converting between indentical types: simply copy + +template <> struct sGDALCopyWord +{ + static inline void f(const GFloat16 hfValueIn, GFloat16 &hfValueOut) + { + hfValueOut = hfValueIn; + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const float fValueIn, float &fValueOut) + { + fValueOut = fValueIn; + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const double dfValueIn, double &dfValueOut) + { + dfValueOut = dfValueIn; + } +}; + +// Floating point input and output, converting to a larger type: use implicit conversion + +template <> struct sGDALCopyWord +{ + static inline void f(const GFloat16 hfValueIn, float &dfValueOut) + { + dfValueOut = hfValueIn; + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const GFloat16 hfValueIn, double &dfValueOut) + { + dfValueOut = hfValueIn; + } +}; + template <> struct sGDALCopyWord { static inline void f(const float fValueIn, double &dfValueOut) @@ -235,6 +309,46 @@ template <> struct sGDALCopyWord } }; +// Floating point input and out, converting to a smaller type: ensure overflow results in infinity + +template <> struct sGDALCopyWord +{ + static inline void f(const float fValueIn, GFloat16 &hfValueOut) + { + if (fValueIn > cpl::NumericLimits::max()) + { + hfValueOut = cpl::NumericLimits::infinity(); + return; + } + if (fValueIn < -cpl::NumericLimits::max()) + { + hfValueOut = -cpl::NumericLimits::infinity(); + return; + } + + hfValueOut = static_cast(fValueIn); + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const double dfValueIn, GFloat16 &hfValueOut) + { + if (dfValueIn > cpl::NumericLimits::max()) + { + hfValueOut = cpl::NumericLimits::infinity(); + return; + } + if (dfValueIn < -cpl::NumericLimits::max()) + { + hfValueOut = -cpl::NumericLimits::infinity(); + return; + } + + hfValueOut = static_cast(dfValueIn); + } +}; + template <> struct sGDALCopyWord { static inline void f(const double dfValueIn, float &fValueOut) @@ -254,53 +368,37 @@ template <> struct sGDALCopyWord } }; -template struct sGDALCopyWord +// Floating point input to a small unsigned integer type: nan becomes zero, otherwise round and clamp + +template struct sGDALCopyWord { - static inline void f(const float fValueIn, Tout &tValueOut) + static inline void f(const GFloat16 hfValueIn, Tout &tValueOut) { - if (std::isnan(fValueIn)) + if (CPLIsNan(hfValueIn)) { tValueOut = 0; return; } - float fMaxVal, fMinVal; - GDALGetDataLimits(fMaxVal, fMinVal); + GFloat16 hfMaxVal, hfMinVal; + GDALGetDataLimits(hfMaxVal, hfMinVal); tValueOut = static_cast( - GDALClampValue(fValueIn + 0.5f, fMaxVal, fMinVal)); - } -}; - -template <> struct sGDALCopyWord -{ - static inline void f(const float fValueIn, short &nValueOut) - { - if (std::isnan(fValueIn)) - { - nValueOut = 0; - return; - } - float fMaxVal, fMinVal; - GDALGetDataLimits(fMaxVal, fMinVal); - float fValue = fValueIn >= 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f; - nValueOut = - static_cast(GDALClampValue(fValue, fMaxVal, fMinVal)); + GDALClampValue(hfValueIn + GFloat16(0.5f), hfMaxVal, hfMinVal)); } }; -template <> struct sGDALCopyWord +template struct sGDALCopyWord { - static inline void f(const float fValueIn, signed char &nValueOut) + static inline void f(const float fValueIn, Tout &tValueOut) { - if (std::isnan(fValueIn)) + if (CPLIsNan(fValueIn)) { - nValueOut = 0; + tValueOut = 0; return; } float fMaxVal, fMinVal; - GDALGetDataLimits(fMaxVal, fMinVal); - float fValue = fValueIn >= 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f; - nValueOut = - static_cast(GDALClampValue(fValue, fMaxVal, fMinVal)); + GDALGetDataLimits(fMaxVal, fMinVal); + tValueOut = static_cast( + GDALClampValue(fValueIn + 0.5f, fMaxVal, fMinVal)); } }; @@ -308,7 +406,7 @@ template struct sGDALCopyWord { static inline void f(const double dfValueIn, Tout &tValueOut) { - if (std::isnan(dfValueIn)) + if (CPLIsNan(dfValueIn)) { tValueOut = 0; return; @@ -320,45 +418,64 @@ template struct sGDALCopyWord } }; -template <> struct sGDALCopyWord +// Floating point input to a large unsigned integer type: nan becomes zero, otherwise round and clamp. +// Avoid roundoff while clamping. + +template <> struct sGDALCopyWord { - static inline void f(const double dfValueIn, int &nValueOut) + static inline void f(const GFloat16 hfValueIn, std::uint64_t &nValueOut) { - if (std::isnan(dfValueIn)) + if (!(hfValueIn > 0)) { nValueOut = 0; - return; } - double dfMaxVal, dfMinVal; - GDALGetDataLimits(dfMaxVal, dfMinVal); - double dfValue = dfValueIn >= 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5; - nValueOut = - static_cast(GDALClampValue(dfValue, dfMaxVal, dfMinVal)); + else if (CPLIsInf(hfValueIn)) + { + nValueOut = cpl::NumericLimits::max(); + } + else + { + nValueOut = static_cast(hfValueIn + GFloat16(0.5f)); + } } }; -template <> struct sGDALCopyWord +template <> struct sGDALCopyWord { - static inline void f(const double dfValueIn, std::int64_t &nValueOut) + static inline void f(const float fValueIn, unsigned int &nValueOut) { - if (std::isnan(dfValueIn)) + if (!(fValueIn > 0)) { nValueOut = 0; } - else if (dfValueIn >= - static_cast(std::numeric_limits::max())) + else if (fValueIn >= + static_cast(cpl::NumericLimits::max())) { - nValueOut = std::numeric_limits::max(); + nValueOut = cpl::NumericLimits::max(); } - else if (dfValueIn <= - static_cast(std::numeric_limits::min())) + else { - nValueOut = std::numeric_limits::min(); + nValueOut = static_cast(fValueIn + 0.5f); + } + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const float fValueIn, std::uint64_t &nValueOut) + { + if (!(fValueIn > 0)) + { + nValueOut = 0; + } + else if (fValueIn >= + static_cast(cpl::NumericLimits::max())) + { + nValueOut = cpl::NumericLimits::max(); } else { - nValueOut = static_cast( - dfValueIn > 0.0f ? dfValueIn + 0.5f : dfValueIn - 0.5f); + nValueOut = static_cast(fValueIn + 0.5f); } } }; @@ -372,9 +489,9 @@ template <> struct sGDALCopyWord nValueOut = 0; } else if (dfValueIn > - static_cast(std::numeric_limits::max())) + static_cast(cpl::NumericLimits::max())) { - nValueOut = std::numeric_limits::max(); + nValueOut = cpl::NumericLimits::max(); } else { @@ -383,20 +500,101 @@ template <> struct sGDALCopyWord } }; -template <> struct sGDALCopyWord +// Floating point input to a very large unsigned integer type: nan becomes zero, otherwise round and clamp. +// Avoid infinity while clamping when the maximum integer is too large for the floating-point type. +// Avoid roundoff while clamping. + +template <> struct sGDALCopyWord { - static inline void f(const double dfValueIn, short &nValueOut) + static inline void f(const GFloat16 hfValueIn, unsigned short &nValueOut) + { + if (!(hfValueIn > 0)) + { + nValueOut = 0; + } + else if (CPLIsInf(hfValueIn)) + { + nValueOut = cpl::NumericLimits::max(); + } + else + { + nValueOut = static_cast(hfValueIn + GFloat16(0.5f)); + } + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const GFloat16 hfValueIn, unsigned int &nValueOut) { - if (std::isnan(dfValueIn)) + if (!(hfValueIn > 0)) + { + nValueOut = 0; + } + else if (CPLIsInf(hfValueIn)) + { + nValueOut = cpl::NumericLimits::max(); + } + else + { + nValueOut = static_cast(hfValueIn + GFloat16(0.5f)); + } + } +}; + +// Floating point input to a small signed integer type: nan becomes zero, otherwise round and clamp. +// Rounding for signed integers is different than for the unsigned integers above. + +template <> struct sGDALCopyWord +{ + static inline void f(const GFloat16 hfValueIn, signed char &nValueOut) + { + if (CPLIsNan(hfValueIn)) { nValueOut = 0; return; } - double dfMaxVal, dfMinVal; - GDALGetDataLimits(dfMaxVal, dfMinVal); - double dfValue = dfValueIn > 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5; + GFloat16 hfMaxVal, hfMinVal; + GDALGetDataLimits(hfMaxVal, hfMinVal); + GFloat16 hfValue = hfValueIn >= GFloat16(0.0f) + ? hfValueIn + GFloat16(0.5f) + : hfValueIn - GFloat16(0.5f); + nValueOut = static_cast( + GDALClampValue(hfValue, hfMaxVal, hfMinVal)); + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const float fValueIn, signed char &nValueOut) + { + if (CPLIsNan(fValueIn)) + { + nValueOut = 0; + return; + } + float fMaxVal, fMinVal; + GDALGetDataLimits(fMaxVal, fMinVal); + float fValue = fValueIn >= 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f; nValueOut = - static_cast(GDALClampValue(dfValue, dfMaxVal, dfMinVal)); + static_cast(GDALClampValue(fValue, fMaxVal, fMinVal)); + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const float fValueIn, short &nValueOut) + { + if (CPLIsNan(fValueIn)) + { + nValueOut = 0; + return; + } + float fMaxVal, fMinVal; + GDALGetDataLimits(fMaxVal, fMinVal); + float fValue = fValueIn >= 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f; + nValueOut = + static_cast(GDALClampValue(fValue, fMaxVal, fMinVal)); } }; @@ -404,7 +602,7 @@ template <> struct sGDALCopyWord { static inline void f(const double dfValueIn, signed char &nValueOut) { - if (std::isnan(dfValueIn)) + if (CPLIsNan(dfValueIn)) { nValueOut = 0; return; @@ -417,25 +615,87 @@ template <> struct sGDALCopyWord } }; -// Roundoff occurs for Float32 -> int32 for max/min. Overload GDALCopyWord -// specifically for this case. +template <> struct sGDALCopyWord +{ + static inline void f(const double dfValueIn, short &nValueOut) + { + if (CPLIsNan(dfValueIn)) + { + nValueOut = 0; + return; + } + double dfMaxVal, dfMinVal; + GDALGetDataLimits(dfMaxVal, dfMinVal); + double dfValue = dfValueIn > 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5; + nValueOut = + static_cast(GDALClampValue(dfValue, dfMaxVal, dfMinVal)); + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const double dfValueIn, int &nValueOut) + { + if (CPLIsNan(dfValueIn)) + { + nValueOut = 0; + return; + } + double dfMaxVal, dfMinVal; + GDALGetDataLimits(dfMaxVal, dfMinVal); + double dfValue = dfValueIn >= 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5; + nValueOut = + static_cast(GDALClampValue(dfValue, dfMaxVal, dfMinVal)); + } +}; + +// Floating point input to a large signed integer type: nan becomes zero, otherwise round and clamp. +// Rounding for signed integers is different than for the unsigned integers above. +// Avoid roundoff while clamping. + +template <> struct sGDALCopyWord +{ + static inline void f(const GFloat16 hfValueIn, short &nValueOut) + { + if (CPLIsNan(hfValueIn)) + { + nValueOut = 0; + } + else if (hfValueIn >= + static_cast(cpl::NumericLimits::max())) + { + nValueOut = cpl::NumericLimits::max(); + } + else if (hfValueIn <= + static_cast(cpl::NumericLimits::lowest())) + { + nValueOut = cpl::NumericLimits::lowest(); + } + else + { + nValueOut = static_cast(hfValueIn > GFloat16(0.0f) + ? hfValueIn + GFloat16(0.5f) + : hfValueIn - GFloat16(0.5f)); + } + } +}; + template <> struct sGDALCopyWord { static inline void f(const float fValueIn, int &nValueOut) { - if (std::isnan(fValueIn)) + if (CPLIsNan(fValueIn)) { nValueOut = 0; } - else if (fValueIn >= - static_cast(std::numeric_limits::max())) + else if (fValueIn >= static_cast(cpl::NumericLimits::max())) { - nValueOut = std::numeric_limits::max(); + nValueOut = cpl::NumericLimits::max(); } else if (fValueIn <= - static_cast(std::numeric_limits::min())) + static_cast(cpl::NumericLimits::lowest())) { - nValueOut = std::numeric_limits::min(); + nValueOut = cpl::NumericLimits::lowest(); } else { @@ -445,74 +705,105 @@ template <> struct sGDALCopyWord } }; -// Roundoff occurs for Float32 -> uint32 for max. Overload GDALCopyWord -// specifically for this case. -template <> struct sGDALCopyWord +template <> struct sGDALCopyWord { - static inline void f(const float fValueIn, unsigned int &nValueOut) + static inline void f(const float fValueIn, std::int64_t &nValueOut) { - if (!(fValueIn > 0)) + if (CPLIsNan(fValueIn)) { nValueOut = 0; } else if (fValueIn >= - static_cast(std::numeric_limits::max())) + static_cast(cpl::NumericLimits::max())) + { + nValueOut = cpl::NumericLimits::max(); + } + else if (fValueIn <= + static_cast(cpl::NumericLimits::lowest())) { - nValueOut = std::numeric_limits::max(); + nValueOut = cpl::NumericLimits::lowest(); } else { - nValueOut = static_cast(fValueIn + 0.5f); + nValueOut = static_cast( + fValueIn > 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f); } } }; -// Roundoff occurs for Float32 -> std::int64_t for max/min. Overload -// GDALCopyWord specifically for this case. -template <> struct sGDALCopyWord +template <> struct sGDALCopyWord { - static inline void f(const float fValueIn, std::int64_t &nValueOut) + static inline void f(const double dfValueIn, std::int64_t &nValueOut) { - if (std::isnan(fValueIn)) + if (CPLIsNan(dfValueIn)) { nValueOut = 0; } - else if (fValueIn >= - static_cast(std::numeric_limits::max())) + else if (dfValueIn >= + static_cast(cpl::NumericLimits::max())) { - nValueOut = std::numeric_limits::max(); + nValueOut = cpl::NumericLimits::max(); } - else if (fValueIn <= - static_cast(std::numeric_limits::min())) + else if (dfValueIn <= + static_cast(cpl::NumericLimits::min())) { - nValueOut = std::numeric_limits::min(); + nValueOut = cpl::NumericLimits::min(); } else { nValueOut = static_cast( - fValueIn > 0.0f ? fValueIn + 0.5f : fValueIn - 0.5f); + dfValueIn > 0.0 ? dfValueIn + 0.5 : dfValueIn - 0.5); } } }; -// Roundoff occurs for Float32 -> std::uint64_t for max. Overload GDALCopyWord -// specifically for this case. -template <> struct sGDALCopyWord +// Floating point input to a very large signed integer type: nan becomes zero, otherwise round and clamp. +// Rounding for signed integers is different than for the unsigned integers above. +// Avoid infinity while clamping when the maximum integer is too large for the floating-point type. +// Avoid roundoff while clamping. + +template <> struct sGDALCopyWord { - static inline void f(const float fValueIn, std::uint64_t &nValueOut) + static inline void f(const GFloat16 hfValueIn, int &nValueOut) { - if (!(fValueIn > 0)) + if (CPLIsNan(hfValueIn)) { nValueOut = 0; } - else if (fValueIn >= - static_cast(std::numeric_limits::max())) + else if (CPLIsInf(hfValueIn)) { - nValueOut = std::numeric_limits::max(); + nValueOut = hfValueIn > GFloat16(0.0f) + ? cpl::NumericLimits::max() + : cpl::NumericLimits::lowest(); } else { - nValueOut = static_cast(fValueIn + 0.5f); + nValueOut = static_cast(hfValueIn > GFloat16(0.0f) + ? hfValueIn + GFloat16(0.5f) + : hfValueIn - GFloat16(0.5f)); + } + } +}; + +template <> struct sGDALCopyWord +{ + static inline void f(const GFloat16 hfValueIn, std::int64_t &nValueOut) + { + if (CPLIsNan(hfValueIn)) + { + nValueOut = 0; + } + else if (CPLIsInf(hfValueIn)) + { + nValueOut = hfValueIn > GFloat16(0.0f) + ? cpl::NumericLimits::max() + : cpl::NumericLimits::lowest(); + } + else + { + nValueOut = static_cast( + hfValueIn > GFloat16(0.0f) ? hfValueIn + GFloat16(0.5f) + : hfValueIn - GFloat16(0.5f)); } } }; diff --git a/gcore/gdal_typetraits.h b/gcore/gdal_typetraits.h index f7e341f49e7c..de0553f02357 100644 --- a/gcore/gdal_typetraits.h +++ b/gcore/gdal_typetraits.h @@ -34,7 +34,7 @@ #ifdef GDAL_ENABLE_FLOAT16 #if defined(__GNUC__) || defined(__clang__) #define __STDC_WANT_IEC_60559_TYPES_EXT__ -#include // Also brings in _Float16 +#include // Also brings in GFloat16 #endif #endif @@ -165,6 +165,20 @@ template <> struct CXXTypeTraits } }; +template <> struct CXXTypeTraits +{ + static constexpr GDALDataType gdal_type = GDT_Float16; + static constexpr size_t size = sizeof(GFloat16); + static constexpr OGRFieldType ogr_type = OFTReal; + // We could introduce OFSTFloat16 + static constexpr OGRFieldSubType ogr_subtype = OFSTNone; + + static inline GDALExtendedDataType GetExtendedDataType() + { + return GDALExtendedDataType::Create(GDT_Float16); + } +}; + template <> struct CXXTypeTraits { static constexpr GDALDataType gdal_type = GDT_Float32; @@ -191,46 +205,44 @@ template <> struct CXXTypeTraits } }; -template <> struct CXXTypeTraits> +template <> struct CXXTypeTraits> { - static constexpr GDALDataType gdal_type = GDT_CFloat32; - static constexpr size_t size = sizeof(float) * 2; + static constexpr GDALDataType gdal_type = GDT_CFloat16; + static constexpr size_t size = sizeof(GFloat16) * 2; static constexpr OGRFieldType ogr_type = OFTMaxType; static constexpr OGRFieldSubType ogr_subtype = OFSTNone; static inline GDALExtendedDataType GetExtendedDataType() { - return GDALExtendedDataType::Create(GDT_CFloat32); + return GDALExtendedDataType::Create(GDT_CFloat16); } }; -template <> struct CXXTypeTraits> +template <> struct CXXTypeTraits> { - static constexpr GDALDataType gdal_type = GDT_CFloat64; - static constexpr size_t size = sizeof(double) * 2; + static constexpr GDALDataType gdal_type = GDT_CFloat32; + static constexpr size_t size = sizeof(float) * 2; static constexpr OGRFieldType ogr_type = OFTMaxType; static constexpr OGRFieldSubType ogr_subtype = OFSTNone; static inline GDALExtendedDataType GetExtendedDataType() { - return GDALExtendedDataType::Create(GDT_CFloat64); + return GDALExtendedDataType::Create(GDT_CFloat32); } }; -#if defined(GDAL_ENABLE_FLOAT16) && defined(FLT16_MAX) && defined(FLT16_MIN) -template <> struct CXXTypeTraits<_Float16> +template <> struct CXXTypeTraits> { - static constexpr GDALDataType gdal_type = GDT_Float16; - static constexpr size_t size = sizeof(_Float16); - static constexpr OGRFieldType ogr_type = OFTReal; + static constexpr GDALDataType gdal_type = GDT_CFloat64; + static constexpr size_t size = sizeof(double) * 2; + static constexpr OGRFieldType ogr_type = OFTMaxType; static constexpr OGRFieldSubType ogr_subtype = OFSTNone; static inline GDALExtendedDataType GetExtendedDataType() { - return GDALExtendedDataType::Create(GDT_Float16); + return GDALExtendedDataType::Create(GDT_CFloat64); } }; -#endif template <> struct CXXTypeTraits { @@ -473,11 +485,13 @@ inline OGRFieldType GetOGRFieldType(const GDALDataType gdal_type) case GDT_Int64: return OFTInteger64; case GDT_UInt64: // Questionable + case GDT_Float16: case GDT_Float32: case GDT_Float64: return OFTReal; case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_Unknown: diff --git a/gcore/gdalcachedpixelaccessor.h b/gcore/gdalcachedpixelaccessor.h index a629a02a8336..ef6c1fcf8b69 100644 --- a/gcore/gdalcachedpixelaccessor.h +++ b/gcore/gdalcachedpixelaccessor.h @@ -351,6 +351,11 @@ template <> struct GDALCachedPixelAccessorGetDataType static constexpr GDALDataType DataType = GDT_Int64; }; +template <> struct GDALCachedPixelAccessorGetDataType +{ + static constexpr GDALDataType DataType = GDT_Float16; +}; + template <> struct GDALCachedPixelAccessorGetDataType { static constexpr GDALDataType DataType = GDT_Float32; diff --git a/gcore/gdalmultidim.cpp b/gcore/gdalmultidim.cpp index 2832b5609e21..18437311716b 100644 --- a/gcore/gdalmultidim.cpp +++ b/gcore/gdalmultidim.cpp @@ -23,6 +23,7 @@ #include // isalnum #include "cpl_error_internal.h" +#include "cpl_float.h" #include "gdal_priv.h" #include "gdal_pam.h" #include "gdal_utils.h" @@ -980,7 +981,8 @@ bool GDALGroup::CopyFrom(const std::shared_ptr &poDstRootGroup, bool bHasOffset = false; bool bHasScale = false; if (bAutoScale && srcArrayType.GetClass() == GEDTC_NUMERIC && - (srcArrayType.GetNumericDataType() == GDT_Float32 || + (srcArrayType.GetNumericDataType() == GDT_Float16 || + srcArrayType.GetNumericDataType() == GDT_Float32 || srcArrayType.GetNumericDataType() == GDT_Float64) && srcArray->GetOffset(&bHasOffset) == 0.0 && !bHasOffset && srcArray->GetScale(&bHasScale) == 1.0 && !bHasScale && @@ -1004,8 +1006,8 @@ bool GDALGroup::CopyFrom(const std::shared_ptr &poDstRootGroup, #define setDTMinMax(ctype) \ do \ { \ - dfDTMin = static_cast(std::numeric_limits::min()); \ - dfDTMax = static_cast(std::numeric_limits::max()); \ + dfDTMin = static_cast(cpl::NumericLimits::lowest()); \ + dfDTMax = static_cast(cpl::NumericLimits::max()); \ } while (0) switch (eAutoScaleType) @@ -1034,11 +1036,13 @@ bool GDALGroup::CopyFrom(const std::shared_ptr &poDstRootGroup, case GDT_Int64: setDTMinMax(std::int64_t); break; + case GDT_Float16: case GDT_Float32: case GDT_Float64: case GDT_Unknown: case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_TypeCount: @@ -1673,6 +1677,10 @@ bool GDALExtendedDataType::CopyValue(const void *pSrc, static_cast( *static_cast(pSrc))); break; + case GDT_Float16: + str = CPLSPrintf("%.5g", + double(*static_cast(pSrc))); + break; case GDT_Float32: str = CPLSPrintf("%.9g", *static_cast(pSrc)); break; @@ -1691,6 +1699,12 @@ bool GDALExtendedDataType::CopyValue(const void *pSrc, str = CPLSPrintf("%d+%dj", src[0], src[1]); break; } + case GDT_CFloat16: + { + const GFloat16 *src = static_cast(pSrc); + str = CPLSPrintf("%.5g+%.5gj", double(src[0]), double(src[1])); + break; + } case GDT_CFloat32: { const float *src = static_cast(pSrc); @@ -6559,9 +6573,13 @@ GDALMDArray::GetUnscaled(double dfOverriddenScale, double dfOverriddenOffset, GDALDataType eDT = GDALDataTypeIsComplex(GetDataType().GetNumericDataType()) ? GDT_CFloat64 : GDT_Float64; - if (dfOverriddenScale == -1 && dfOverriddenOffset == 0 && - GetDataType().GetNumericDataType() == GDT_Float32) - eDT = GDT_Float32; + if (dfOverriddenScale == -1 && dfOverriddenOffset == 0) + { + if (GetDataType().GetNumericDataType() == GDT_Float16) + eDT = GDT_Float16; + if (GetDataType().GetNumericDataType() == GDT_Float32) + eDT = GDT_Float32; + } return GDALMDArrayUnscaled::Create(self, dfScale, dfOffset, dfOverriddenDstNodata, eDT); @@ -7061,6 +7079,12 @@ bool GDALMDArrayMask::IRead(const GUInt64 *arrayStartIdx, const size_t *count, tmpBufferStrideVector); break; + case GDT_Float16: + ReadInternal(count, bufferStride, bufferDataType, + pDstBuffer, pTempBuffer, oTmpBufferDT, + tmpBufferStrideVector); + break; + case GDT_Float32: ReadInternal(count, bufferStride, bufferDataType, pDstBuffer, pTempBuffer, oTmpBufferDT, @@ -7075,6 +7099,7 @@ bool GDALMDArrayMask::IRead(const GUInt64 *arrayStartIdx, const size_t *count, case GDT_Unknown: case GDT_CInt16: case GDT_CInt32: + case GDT_CFloat16: case GDT_CFloat32: case GDT_CFloat64: case GDT_TypeCount: @@ -7095,9 +7120,9 @@ template static bool IsValidForDT(double dfVal) { if (std::isnan(dfVal)) return false; - if (dfVal < static_cast(std::numeric_limits::lowest())) + if (dfVal < static_cast(cpl::NumericLimits::lowest())) return false; - if (dfVal > static_cast(std::numeric_limits::max())) + if (dfVal > static_cast(cpl::NumericLimits::max())) return false; return static_cast(static_cast(dfVal)) == dfVal; } @@ -9911,8 +9936,8 @@ bool GDALMDArray::ComputeStatistics(bool bApproxOK, double *pdfMin, { const GDALMDArray *array = nullptr; std::shared_ptr poMask{}; - double dfMin = std::numeric_limits::max(); - double dfMax = -std::numeric_limits::max(); + double dfMin = cpl::NumericLimits::max(); + double dfMax = -cpl::NumericLimits::max(); double dfMean = 0.0; double dfM2 = 0.0; GUInt64 nValidCount = 0; diff --git a/gcore/gdalmultidim_rat.cpp b/gcore/gdalmultidim_rat.cpp index b05cab584eba..1234e097932b 100644 --- a/gcore/gdalmultidim_rat.cpp +++ b/gcore/gdalmultidim_rat.cpp @@ -77,6 +77,7 @@ class GDALRasterAttributeTableFromMDArrays final case GDT_UInt32: case GDT_Int64: case GDT_UInt64: + case GDT_Float16: case GDT_Float32: case GDT_Float64: return GFT_Real; diff --git a/gcore/gdalnodatamaskband.cpp b/gcore/gdalnodatamaskband.cpp index 953f100bc96a..3805ff010390 100644 --- a/gcore/gdalnodatamaskband.cpp +++ b/gcore/gdalnodatamaskband.cpp @@ -116,6 +116,8 @@ static GDALDataType GetWorkDataType(GDALDataType eDataType) eWrkDT = GDT_Int32; break; + case GDT_Float16: + case GDT_CFloat16: case GDT_Float32: case GDT_CFloat32: eWrkDT = GDT_Float32; @@ -184,6 +186,12 @@ bool GDALNoDataMaskBand::IsNoDataInRange(double dfNoDataValue, return GDALIsValueInRange(dfNoDataValue); } + case GDT_Float16: + { + return isnan(dfNoDataValue) || isinf(dfNoDataValue) || + GDALIsValueInRange(dfNoDataValue); + } + case GDT_Float32: { return std::isnan(dfNoDataValue) || std::isinf(dfNoDataValue) || diff --git a/gcore/gdalnodatavaluesmaskband.cpp b/gcore/gdalnodatavaluesmaskband.cpp index bf3eb0421c12..1843aa1b39f1 100644 --- a/gcore/gdalnodatavaluesmaskband.cpp +++ b/gcore/gdalnodatavaluesmaskband.cpp @@ -127,6 +127,8 @@ CPLErr GDALNoDataValuesMaskBand::IReadBlock(int nXBlockOff, int nYBlockOff, eWrkDT = GDT_Int32; break; + case GDT_Float16: + case GDT_CFloat16: case GDT_Float32: case GDT_CFloat32: eWrkDT = GDT_Float32; diff --git a/gcore/gdalrasterband.cpp b/gcore/gdalrasterband.cpp index 4185b65cab29..3e26a211bfe0 100644 --- a/gcore/gdalrasterband.cpp +++ b/gcore/gdalrasterband.cpp @@ -30,6 +30,7 @@ #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_progress.h" #include "cpl_string.h" #include "cpl_virtualmem.h" @@ -2810,6 +2811,10 @@ double GDALRasterBand::GetMaximum(int *pbSuccess) case GDT_UInt64: return static_cast(std::numeric_limits::max()); + case GDT_Float16: + case GDT_CFloat16: + return 65504.0; + case GDT_Float32: case GDT_CFloat32: return 4294967295.0; // Not actually accurate. @@ -2911,11 +2916,15 @@ double GDALRasterBand::GetMinimum(int *pbSuccess) return 0; case GDT_Int64: - return static_cast(std::numeric_limits::min()); + return static_cast(std::numeric_limits::lowest()); case GDT_UInt64: return 0; + case GDT_Float16: + case GDT_CFloat16: + return -65504.0; + case GDT_Float32: case GDT_CFloat32: return -4294967295.0; // Not actually accurate. @@ -3870,6 +3879,28 @@ GDALDatasetH CPL_STDCALL GDALGetBandDataset(GDALRasterBandH hBand) return GDALDataset::ToHandle(poBand->GetDataset()); } +/************************************************************************/ +/* ComputeFloat16NoDataValue() */ +/************************************************************************/ + +static inline void ComputeFloat16NoDataValue(GDALDataType eDataType, + double dfNoDataValue, + int &bGotNoDataValue, + GFloat16 &fNoDataValue, + bool &bGotFloat16NoDataValue) +{ + if (eDataType == GDT_Float16 && bGotNoDataValue) + { + dfNoDataValue = GDALAdjustNoDataCloseToFloatMax(dfNoDataValue); + if (GDALIsValueInRange(dfNoDataValue)) + { + fNoDataValue = static_cast(dfNoDataValue); + bGotFloat16NoDataValue = true; + bGotNoDataValue = false; + } + } +} + /************************************************************************/ /* ComputeFloatNoDataValue() */ /************************************************************************/ @@ -3892,6 +3923,45 @@ static inline void ComputeFloatNoDataValue(GDALDataType eDataType, } } +/************************************************************************/ +/* struct GDALNoDataValues */ +/************************************************************************/ + +/** + * \brief No-data-values for all types + * + * The functions below pass various no-data-values around. To avoid + * long argument lists, this struct collects the no-data-values for + * all types into a single, convenient place. + **/ + +struct GDALNoDataValues +{ + int bGotNoDataValue; + double dfNoDataValue; + + bool bGotFloatNoDataValue; + float fNoDataValue; + + bool bGotFloat16NoDataValue; + GFloat16 hfNoDataValue; + + GDALNoDataValues(GDALRasterBand *poRasterBand, GDALDataType eDataType) + : bGotNoDataValue(FALSE), dfNoDataValue(0.0), + bGotFloatNoDataValue(false), fNoDataValue(0.0f), + bGotFloat16NoDataValue(false), hfNoDataValue(GFloat16(0.0f)) + { + dfNoDataValue = poRasterBand->GetNoDataValue(&bGotNoDataValue); + bGotNoDataValue = bGotNoDataValue && !CPLIsNan(dfNoDataValue); + + ComputeFloatNoDataValue(eDataType, dfNoDataValue, bGotNoDataValue, + fNoDataValue, bGotFloatNoDataValue); + + ComputeFloat16NoDataValue(eDataType, dfNoDataValue, bGotNoDataValue, + hfNoDataValue, bGotFloat16NoDataValue); + } +}; + /************************************************************************/ /* GetHistogram() */ /************************************************************************/ @@ -3995,15 +4065,9 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, } memset(panHistogram, 0, sizeof(GUIntBig) * nBuckets); - int bGotNoDataValue = FALSE; - const double dfNoDataValue = GetNoDataValue(&bGotNoDataValue); - bGotNoDataValue = bGotNoDataValue && !std::isnan(dfNoDataValue); - bool bGotFloatNoDataValue = false; - float fNoDataValue = 0.0f; - ComputeFloatNoDataValue(eDataType, dfNoDataValue, bGotNoDataValue, - fNoDataValue, bGotFloatNoDataValue); + GDALNoDataValues sNoDataValues(this, eDataType); GDALRasterBand *poMaskBand = nullptr; - if (!bGotNoDataValue) + if (!sNoDataValues.bGotNoDataValue) { const int l_nMaskFlags = GetMaskFlags(); if (l_nMaskFlags != GMF_ALL_VALID && l_nMaskFlags != GMF_NODATA && @@ -4130,13 +4194,26 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, dfValue = static_cast( static_cast(pData)[iOffset]); break; + case GDT_Float16: + { + const GFloat16 hfValue = + static_cast(pData)[iOffset]; + if (CPLIsNan(hfValue) || + (sNoDataValues.bGotFloat16NoDataValue && + ARE_REAL_EQUAL(hfValue, + sNoDataValues.hfNoDataValue))) + continue; + dfValue = hfValue; + break; + } case GDT_Float32: { const float fValue = static_cast(pData)[iOffset]; - if (std::isnan(fValue) || - (bGotFloatNoDataValue && - ARE_REAL_EQUAL(fValue, fNoDataValue))) + if (CPLIsNan(fValue) || + (sNoDataValues.bGotFloatNoDataValue && + ARE_REAL_EQUAL(fValue, + sNoDataValues.fNoDataValue))) continue; dfValue = fValue; break; @@ -4168,6 +4245,17 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); } break; + case GDT_CFloat16: + { + const double dfReal = + static_cast(pData)[iOffset * 2]; + const double dfImag = + static_cast(pData)[iOffset * 2 + 1]; + if (CPLIsNan(dfReal) || CPLIsNan(dfImag)) + continue; + dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; + } case GDT_CFloat32: { const double dfReal = @@ -4177,8 +4265,8 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, if (std::isnan(dfReal) || std::isnan(dfImag)) continue; dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; } - break; case GDT_CFloat64: { const double dfReal = @@ -4188,15 +4276,16 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, if (std::isnan(dfReal) || std::isnan(dfImag)) continue; dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; } - break; case GDT_Unknown: case GDT_TypeCount: CPLAssert(false); } - if (eDataType != GDT_Float32 && bGotNoDataValue && - ARE_REAL_EQUAL(dfValue, dfNoDataValue)) + if (eDataType != GDT_Float16 && eDataType != GDT_Float32 && + sNoDataValues.bGotNoDataValue && + ARE_REAL_EQUAL(dfValue, sNoDataValues.dfNoDataValue)) continue; // Given that dfValue and dfMin are not NaN, and dfScale > 0 and @@ -4317,8 +4406,9 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, { if (pabyMaskData && pabyMaskData[i] == 0) continue; - if (!(bGotNoDataValue && - (pabyData[i] == static_cast(dfNoDataValue)))) + if (!(sNoDataValues.bGotNoDataValue && + (pabyData[i] == + static_cast(sNoDataValues.dfNoDataValue)))) { panHistogram[pabyData[i]]++; } @@ -4375,13 +4465,26 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, dfValue = static_cast( static_cast(pData)[iOffset]); break; + case GDT_Float16: + { + const GFloat16 hfValue = + static_cast(pData)[iOffset]; + if (CPLIsNan(hfValue) || + (sNoDataValues.bGotFloat16NoDataValue && + ARE_REAL_EQUAL(hfValue, + sNoDataValues.hfNoDataValue))) + continue; + dfValue = hfValue; + break; + } case GDT_Float32: { const float fValue = static_cast(pData)[iOffset]; - if (std::isnan(fValue) || - (bGotFloatNoDataValue && - ARE_REAL_EQUAL(fValue, fNoDataValue))) + if (CPLIsNan(fValue) || + (sNoDataValues.bGotFloatNoDataValue && + ARE_REAL_EQUAL(fValue, + sNoDataValues.fNoDataValue))) continue; dfValue = fValue; break; @@ -4398,8 +4501,8 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, double dfImag = static_cast(pData)[iOffset * 2 + 1]; dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; } - break; case GDT_CInt32: { double dfReal = @@ -4407,8 +4510,19 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, double dfImag = static_cast(pData)[iOffset * 2 + 1]; dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; + } + case GDT_CFloat16: + { + double dfReal = + static_cast(pData)[iOffset * 2]; + double dfImag = + static_cast(pData)[iOffset * 2 + 1]; + if (CPLIsNan(dfReal) || CPLIsNan(dfImag)) + continue; + dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; } - break; case GDT_CFloat32: { double dfReal = @@ -4418,8 +4532,8 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, if (std::isnan(dfReal) || std::isnan(dfImag)) continue; dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; } - break; case GDT_CFloat64: { double dfReal = @@ -4429,8 +4543,8 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, if (std::isnan(dfReal) || std::isnan(dfImag)) continue; dfValue = sqrt(dfReal * dfReal + dfImag * dfImag); + break; } - break; case GDT_Unknown: case GDT_TypeCount: CPLAssert(false); @@ -4438,8 +4552,9 @@ CPLErr GDALRasterBand::GetHistogram(double dfMin, double dfMax, int nBuckets, return CE_Failure; } - if (eDataType != GDT_Float32 && bGotNoDataValue && - ARE_REAL_EQUAL(dfValue, dfNoDataValue)) + if (eDataType != GDT_Float16 && eDataType != GDT_Float32 && + sNoDataValues.bGotNoDataValue && + ARE_REAL_EQUAL(dfValue, sNoDataValues.dfNoDataValue)) continue; // Given that dfValue and dfMin are not NaN, and dfScale > 0 @@ -5132,7 +5247,7 @@ struct ComputeStatisticsInternalGeneric nSampleCount += static_cast(nXCheck) * nYCheck; } } - else if (nMin == std::numeric_limits::min() && + else if (nMin == std::numeric_limits::lowest() && nMax == std::numeric_limits::max()) { if constexpr (COMPUTE_OTHER_STATS) @@ -6025,9 +6140,8 @@ struct ComputeStatisticsInternal static inline double GetPixelValue(GDALDataType eDataType, bool bSignedByte, const void *pData, GPtrDiff_t iOffset, - bool bGotNoDataValue, double dfNoDataValue, - bool bGotFloatNoDataValue, - float fNoDataValue, bool &bValid) + const GDALNoDataValues &sNoDataValues, + bool &bValid) { bValid = true; double dfValue = 0; @@ -6064,11 +6178,26 @@ static inline double GetPixelValue(GDALDataType eDataType, bool bSignedByte, dfValue = static_cast( static_cast(pData)[iOffset]); break; + case GDT_Float16: + { + const GFloat16 hfValue = + static_cast(pData)[iOffset]; + if (CPLIsNan(hfValue) || + (sNoDataValues.bGotFloat16NoDataValue && + ARE_REAL_EQUAL(hfValue, sNoDataValues.hfNoDataValue))) + { + bValid = false; + return 0.0; + } + dfValue = hfValue; + return dfValue; + } case GDT_Float32: { const float fValue = static_cast(pData)[iOffset]; - if (std::isnan(fValue) || - (bGotFloatNoDataValue && ARE_REAL_EQUAL(fValue, fNoDataValue))) + if (CPLIsNan(fValue) || + (sNoDataValues.bGotFloatNoDataValue && + ARE_REAL_EQUAL(fValue, sNoDataValues.fNoDataValue))) { bValid = false; return 0.0; @@ -6090,6 +6219,14 @@ static inline double GetPixelValue(GDALDataType eDataType, bool bSignedByte, case GDT_CInt32: dfValue = static_cast(pData)[iOffset * 2]; break; + case GDT_CFloat16: + dfValue = static_cast(pData)[iOffset * 2]; + if (isnan(dfValue)) + { + bValid = false; + return 0.0; + } + break; case GDT_CFloat32: dfValue = static_cast(pData)[iOffset * 2]; if (std::isnan(dfValue)) @@ -6112,7 +6249,8 @@ static inline double GetPixelValue(GDALDataType eDataType, bool bSignedByte, break; } - if (bGotNoDataValue && ARE_REAL_EQUAL(dfValue, dfNoDataValue)) + if (sNoDataValues.bGotNoDataValue && + ARE_REAL_EQUAL(dfValue, sNoDataValues.dfNoDataValue)) { bValid = false; return 0.0; @@ -6277,16 +6415,9 @@ CPLErr GDALRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin, GDALRasterIOExtraArg sExtraArg; INIT_RASTERIO_EXTRA_ARG(sExtraArg); - int bGotNoDataValue = FALSE; - const double dfNoDataValue = GetNoDataValue(&bGotNoDataValue); - bGotNoDataValue = bGotNoDataValue && !std::isnan(dfNoDataValue); - bool bGotFloatNoDataValue = false; - float fNoDataValue = 0.0f; - ComputeFloatNoDataValue(eDataType, dfNoDataValue, bGotNoDataValue, - fNoDataValue, bGotFloatNoDataValue); - + GDALNoDataValues sNoDataValues(this, eDataType); GDALRasterBand *poMaskBand = nullptr; - if (!bGotNoDataValue) + if (!sNoDataValues.bGotNoDataValue) { const int l_nMaskFlags = GetMaskFlags(); if (l_nMaskFlags != GMF_ALL_VALID && l_nMaskFlags != GMF_NODATA && @@ -6378,10 +6509,8 @@ CPLErr GDALRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin, continue; bool bValid = true; - double dfValue = - GetPixelValue(eDataType, bSignedByte, pData, iOffset, - CPL_TO_BOOL(bGotNoDataValue), dfNoDataValue, - bGotFloatNoDataValue, fNoDataValue, bValid); + double dfValue = GetPixelValue(eDataType, bSignedByte, pData, + iOffset, sNoDataValues, bValid); if (!bValid) continue; @@ -6452,11 +6581,13 @@ CPLErr GDALRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin, GUIntBig nSumSquare = 0; // If no valid nodata, map to invalid value (256 for Byte) const GUInt32 nNoDataValue = - (bGotNoDataValue && dfNoDataValue >= 0 && - dfNoDataValue <= nMaxValueType && - fabs(dfNoDataValue - - static_cast(dfNoDataValue + 1e-10)) < 1e-10) - ? static_cast(dfNoDataValue + 1e-10) + (sNoDataValues.bGotNoDataValue && + sNoDataValues.dfNoDataValue >= 0 && + sNoDataValues.dfNoDataValue <= nMaxValueType && + fabs(sNoDataValues.dfNoDataValue - + static_cast(sNoDataValues.dfNoDataValue + + 1e-10)) < 1e-10) + ? static_cast(sNoDataValues.dfNoDataValue + 1e-10) : nMaxValueType + 1; for (GIntBig iSampleBlock = 0; @@ -6630,10 +6761,9 @@ CPLErr GDALRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin, continue; bool bValid = true; - double dfValue = GetPixelValue( - eDataType, bSignedByte, pData, iOffset, - CPL_TO_BOOL(bGotNoDataValue), dfNoDataValue, - bGotFloatNoDataValue, fNoDataValue, bValid); + double dfValue = + GetPixelValue(eDataType, bSignedByte, pData, iOffset, + sNoDataValues, bValid); if (!bValid) continue; @@ -6857,12 +6987,10 @@ static void ComputeMinMax(const T *buffer, size_t nElts, T nodataValue, T *pMin, } template -static void ComputeMinMaxGeneric(const void *pData, int nXCheck, int nYCheck, - int nBlockXSize, bool bGotNoDataValue, - double dfNoDataValue, - bool bGotFloatNoDataValue, float fNoDataValue, - const GByte *pabyMaskData, double &dfMin, - double &dfMax) +static void +ComputeMinMaxGeneric(const void *pData, int nXCheck, int nYCheck, + int nBlockXSize, const GDALNoDataValues &sNoDataValues, + const GByte *pabyMaskData, double &dfMin, double &dfMax) { double dfLocalMin = dfMin; double dfLocalMax = dfMax; @@ -6876,9 +7004,8 @@ static void ComputeMinMaxGeneric(const void *pData, int nXCheck, int nYCheck, if (pabyMaskData && pabyMaskData[iOffset] == 0) continue; bool bValid = true; - double dfValue = GetPixelValue( - eDataType, bSignedByte, pData, iOffset, bGotNoDataValue, - dfNoDataValue, bGotFloatNoDataValue, fNoDataValue, bValid); + double dfValue = GetPixelValue(eDataType, bSignedByte, pData, + iOffset, sNoDataValues, bValid); if (!bValid) continue; @@ -6893,9 +7020,8 @@ static void ComputeMinMaxGeneric(const void *pData, int nXCheck, int nYCheck, static void ComputeMinMaxGeneric(const void *pData, GDALDataType eDataType, bool bSignedByte, int nXCheck, int nYCheck, - int nBlockXSize, bool bGotNoDataValue, - double dfNoDataValue, - bool bGotFloatNoDataValue, float fNoDataValue, + int nBlockXSize, + const GDALNoDataValues &sNoDataValues, const GByte *pabyMaskData, double &dfMin, double &dfMax) { @@ -6908,80 +7034,90 @@ static void ComputeMinMaxGeneric(const void *pData, GDALDataType eDataType, if (bSignedByte) { ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + pData, nXCheck, nYCheck, nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); } else { ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + pData, nXCheck, nYCheck, nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); } break; case GDT_Int8: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_UInt16: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_Int16: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_UInt32: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_Int32: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_UInt64: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_Int64: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); + break; + case GDT_Float16: + ComputeMinMaxGeneric( + pData, nXCheck, nYCheck, nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_Float32: ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, false, 0, - bGotFloatNoDataValue, fNoDataValue, pabyMaskData, dfMin, dfMax); + pData, nXCheck, nYCheck, nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_Float64: ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + pData, nXCheck, nYCheck, nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_CInt16: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_CInt32: - ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, nXCheck, nYCheck, + nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); + break; + case GDT_CFloat16: + ComputeMinMaxGeneric( + pData, nXCheck, nYCheck, nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_CFloat32: ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + pData, nXCheck, nYCheck, nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_CFloat64: ComputeMinMaxGeneric( - pData, nXCheck, nYCheck, nBlockXSize, bGotNoDataValue, - dfNoDataValue, false, 0, pabyMaskData, dfMin, dfMax); + pData, nXCheck, nYCheck, nBlockXSize, sNoDataValues, + pabyMaskData, dfMin, dfMax); break; case GDT_TypeCount: CPLAssert(false); @@ -6992,9 +7128,8 @@ static void ComputeMinMaxGeneric(const void *pData, GDALDataType eDataType, static bool ComputeMinMaxGenericIterBlocks( GDALRasterBand *poBand, GDALDataType eDataType, bool bSignedByte, GIntBig nTotalBlocks, int nSampleRate, int nBlocksPerRow, - bool bGotNoDataValue, double dfNoDataValue, bool bGotFloatNoDataValue, - float fNoDataValue, GDALRasterBand *poMaskBand, double &dfMin, - double &dfMax) + const GDALNoDataValues &sNoDataValues, GDALRasterBand *poMaskBand, + double &dfMin, double &dfMax) { GByte *pabyMaskData = nullptr; @@ -7041,9 +7176,8 @@ static bool ComputeMinMaxGenericIterBlocks( } ComputeMinMaxGeneric(pData, eDataType, bSignedByte, nXCheck, nYCheck, - nBlockXSize, CPL_TO_BOOL(bGotNoDataValue), - dfNoDataValue, bGotFloatNoDataValue, fNoDataValue, - pabyMaskData, dfMin, dfMax); + nBlockXSize, sNoDataValues, pabyMaskData, dfMin, + dfMax); poBlock->DropLock(); } @@ -7110,16 +7244,9 @@ CPLErr GDALRasterBand::ComputeRasterMinMax(int bApproxOK, double *adfMinMax) /* -------------------------------------------------------------------- */ /* Read actual data and compute minimum and maximum. */ /* -------------------------------------------------------------------- */ - int bGotNoDataValue = FALSE; - const double dfNoDataValue = GetNoDataValue(&bGotNoDataValue); - bGotNoDataValue = bGotNoDataValue && !std::isnan(dfNoDataValue); - bool bGotFloatNoDataValue = false; - float fNoDataValue = 0.0f; - ComputeFloatNoDataValue(eDataType, dfNoDataValue, bGotNoDataValue, - fNoDataValue, bGotFloatNoDataValue); - + GDALNoDataValues sNoDataValues(this, eDataType); GDALRasterBand *poMaskBand = nullptr; - if (!bGotNoDataValue) + if (!sNoDataValues.bGotNoDataValue) { const int l_nMaskFlags = GetMaskFlags(); if (l_nMaskFlags != GMF_ALL_VALID && l_nMaskFlags != GMF_NODATA && @@ -7154,23 +7281,26 @@ CPLErr GDALRasterBand::ComputeRasterMinMax(int bApproxOK, double *adfMinMax) double dfMin = std::numeric_limits::max(); // used for generic code path double dfMax = - -std::numeric_limits::max(); // used for generic code path + std::numeric_limits::lowest(); // used for generic code path const bool bUseOptimizedPath = !poMaskBand && ((eDataType == GDT_Byte && !bSignedByte) || eDataType == GDT_Int16 || eDataType == GDT_UInt16); const auto ComputeMinMaxForBlock = - [this, bSignedByte, bGotNoDataValue, dfNoDataValue, &nMin, &nMax, - &nMinInt16, &nMaxInt16](const void *pData, int nXCheck, - int nBufferWidth, int nYCheck) + [this, bSignedByte, &sNoDataValues, &nMin, &nMax, &nMinInt16, + &nMaxInt16](const void *pData, int nXCheck, int nBufferWidth, + int nYCheck) { if (eDataType == GDT_Byte && !bSignedByte) { const bool bHasNoData = - bGotNoDataValue && GDALIsValueInRange(dfNoDataValue) && - static_cast(dfNoDataValue) == dfNoDataValue; + sNoDataValues.bGotNoDataValue && + GDALIsValueInRange(sNoDataValues.dfNoDataValue) && + static_cast(sNoDataValues.dfNoDataValue) == + sNoDataValues.dfNoDataValue; const GUInt32 nNoDataValue = - bHasNoData ? static_cast(dfNoDataValue) : 0; + bHasNoData ? static_cast(sNoDataValues.dfNoDataValue) + : 0; GUIntBig nSum, nSumSquare, nSampleCount, nValidCount; // unused ComputeStatisticsInternal:: @@ -7181,10 +7311,13 @@ CPLErr GDALRasterBand::ComputeRasterMinMax(int bApproxOK, double *adfMinMax) else if (eDataType == GDT_UInt16) { const bool bHasNoData = - bGotNoDataValue && GDALIsValueInRange(dfNoDataValue) && - static_cast(dfNoDataValue) == dfNoDataValue; + sNoDataValues.bGotNoDataValue && + GDALIsValueInRange(sNoDataValues.dfNoDataValue) && + static_cast(sNoDataValues.dfNoDataValue) == + sNoDataValues.dfNoDataValue; const GUInt32 nNoDataValue = - bHasNoData ? static_cast(dfNoDataValue) : 0; + bHasNoData ? static_cast(sNoDataValues.dfNoDataValue) + : 0; GUIntBig nSum, nSumSquare, nSampleCount, nValidCount; // unused ComputeStatisticsInternal:: @@ -7195,12 +7328,14 @@ CPLErr GDALRasterBand::ComputeRasterMinMax(int bApproxOK, double *adfMinMax) else if (eDataType == GDT_Int16) { const bool bHasNoData = - bGotNoDataValue && GDALIsValueInRange(dfNoDataValue) && - static_cast(dfNoDataValue) == dfNoDataValue; + sNoDataValues.bGotNoDataValue && + GDALIsValueInRange(sNoDataValues.dfNoDataValue) && + static_cast(sNoDataValues.dfNoDataValue) == + sNoDataValues.dfNoDataValue; if (bHasNoData) { const int16_t nNoDataValue = - static_cast(dfNoDataValue); + static_cast(sNoDataValues.dfNoDataValue); for (int iY = 0; iY < nYCheck; iY++) { ComputeMinMax( @@ -7286,10 +7421,9 @@ CPLErr GDALRasterBand::ComputeRasterMinMax(int bApproxOK, double *adfMinMax) } else { - ComputeMinMaxGeneric( - pData, eDataType, bSignedByte, nXReduced, nYReduced, nXReduced, - CPL_TO_BOOL(bGotNoDataValue), dfNoDataValue, - bGotFloatNoDataValue, fNoDataValue, pabyMaskData, dfMin, dfMax); + ComputeMinMaxGeneric(pData, eDataType, bSignedByte, nXReduced, + nYReduced, nXReduced, sNoDataValues, + pabyMaskData, dfMin, dfMax); } CPLFree(pData); @@ -7357,9 +7491,7 @@ CPLErr GDALRasterBand::ComputeRasterMinMax(int bApproxOK, double *adfMinMax) static_cast(nBlocksPerRow) * nBlocksPerColumn; if (!ComputeMinMaxGenericIterBlocks( this, eDataType, bSignedByte, nTotalBlocks, nSampleRate, - nBlocksPerRow, CPL_TO_BOOL(bGotNoDataValue), dfNoDataValue, - bGotFloatNoDataValue, fNoDataValue, poMaskBand, dfMin, - dfMax)) + nBlocksPerRow, sNoDataValues, poMaskBand, dfMin, dfMax)) { return CE_Failure; } @@ -7478,16 +7610,9 @@ CPLErr GDALRasterBand::ComputeRasterMinMaxLocation(double *pdfMin, if (!InitBlockInfo()) return CE_Failure; - int bGotNoDataValue = FALSE; - const double dfNoDataValue = GetNoDataValue(&bGotNoDataValue); - bGotNoDataValue = bGotNoDataValue && !std::isnan(dfNoDataValue); - bool bGotFloatNoDataValue = false; - float fNoDataValue = 0.0f; - ComputeFloatNoDataValue(eDataType, dfNoDataValue, bGotNoDataValue, - fNoDataValue, bGotFloatNoDataValue); - + GDALNoDataValues sNoDataValues(this, eDataType); GDALRasterBand *poMaskBand = nullptr; - if (!bGotNoDataValue) + if (!sNoDataValues.bGotNoDataValue) { const int l_nMaskFlags = GetMaskFlags(); if (l_nMaskFlags != GMF_ALL_VALID && l_nMaskFlags != GMF_NODATA && @@ -7562,10 +7687,9 @@ CPLErr GDALRasterBand::ComputeRasterMinMaxLocation(double *pdfMin, if (pabyMaskData && pabyMaskData[iOffset] == 0) continue; bool bValid = true; - double dfValue = GetPixelValue( - eDataType, bSignedByte, pData, iOffset, bGotNoDataValue, - dfNoDataValue, bGotFloatNoDataValue, fNoDataValue, - bValid); + double dfValue = + GetPixelValue(eDataType, bSignedByte, pData, iOffset, + sNoDataValues, bValid); if (!bValid) continue; if (dfValue < dfMin) @@ -7592,19 +7716,22 @@ CPLErr GDALRasterBand::ComputeRasterMinMaxLocation(double *pdfMin, { std::tie(pos_min, pos_max) = gdal::minmax_element( pData, static_cast(nBlockXSize) * nBlockYSize, - eEffectiveDT, bGotNoDataValue, dfNoDataValue); + eEffectiveDT, sNoDataValues.bGotNoDataValue, + sNoDataValues.dfNoDataValue); } else if (bNeedsMin) { pos_min = gdal::min_element( pData, static_cast(nBlockXSize) * nBlockYSize, - eEffectiveDT, bGotNoDataValue, dfNoDataValue); + eEffectiveDT, sNoDataValues.bGotNoDataValue, + sNoDataValues.dfNoDataValue); } else if (bNeedsMax) { pos_max = gdal::max_element( pData, static_cast(nBlockXSize) * nBlockYSize, - eEffectiveDT, bGotNoDataValue, dfNoDataValue); + eEffectiveDT, sNoDataValues.bGotNoDataValue, + sNoDataValues.dfNoDataValue); } if (bNeedsMin) @@ -7612,9 +7739,9 @@ CPLErr GDALRasterBand::ComputeRasterMinMaxLocation(double *pdfMin, const int nMinXBlock = static_cast(pos_min % nBlockXSize); const int nMinYBlock = static_cast(pos_min / nBlockXSize); bool bValid = true; - const double dfMinValueBlock = GetPixelValue( - eDataType, bSignedByte, pData, pos_min, bGotNoDataValue, - dfNoDataValue, bGotFloatNoDataValue, fNoDataValue, bValid); + const double dfMinValueBlock = + GetPixelValue(eDataType, bSignedByte, pData, pos_min, + sNoDataValues, bValid); if (bValid && dfMinValueBlock < dfMin) { dfMin = dfMinValueBlock; @@ -7628,9 +7755,9 @@ CPLErr GDALRasterBand::ComputeRasterMinMaxLocation(double *pdfMin, const int nMaxXBlock = static_cast(pos_max % nBlockXSize); const int nMaxYBlock = static_cast(pos_max / nBlockXSize); bool bValid = true; - const double dfMaxValueBlock = GetPixelValue( - eDataType, bSignedByte, pData, pos_max, bGotNoDataValue, - dfNoDataValue, bGotFloatNoDataValue, fNoDataValue, bValid); + const double dfMaxValueBlock = + GetPixelValue(eDataType, bSignedByte, pData, pos_max, + sNoDataValues, bValid); if (bValid && dfMaxValueBlock > dfMax) { dfMax = dfMaxValueBlock; @@ -8039,14 +8166,15 @@ GDALRasterBand *GDALRasterBand::GetMaskBand() /* -------------------------------------------------------------------- */ if (poDS != nullptr) { - const char *pszNoDataValues = poDS->GetMetadataItem("NODATA_VALUES"); - if (pszNoDataValues != nullptr) + const char *pszGDALNoDataValues = + poDS->GetMetadataItem("NODATA_VALUES"); + if (pszGDALNoDataValues != nullptr) { - char **papszNoDataValues = - CSLTokenizeStringComplex(pszNoDataValues, " ", FALSE, FALSE); + char **papszGDALNoDataValues = CSLTokenizeStringComplex( + pszGDALNoDataValues, " ", FALSE, FALSE); // Make sure we have as many values as bands. - if (CSLCount(papszNoDataValues) == poDS->GetRasterCount() && + if (CSLCount(papszGDALNoDataValues) == poDS->GetRasterCount() && poDS->GetRasterCount() != 0) { // Make sure that all bands have the same data type @@ -8076,7 +8204,7 @@ GDALRasterBand *GDALRasterBand::GetMaskBand() CPLError(CE_Failure, CPLE_OutOfMemory, "Out of memory"); poMask.reset(); } - CSLDestroy(papszNoDataValues); + CSLDestroy(papszGDALNoDataValues); return poMask.get(); } else @@ -8096,7 +8224,7 @@ GDALRasterBand *GDALRasterBand::GetMaskBand() "Ignoring it for mask."); } - CSLDestroy(papszNoDataValues); + CSLDestroy(papszGDALNoDataValues); } } @@ -8757,7 +8885,7 @@ void GDALRasterBand::ReportError(CPLErr eErrClass, CPLErrorNum err_no, * uncompressed, scanline oriented (i.e. not tiled). Strips must be organized in * the file in sequential order, and be equally spaced (which is generally the * case). Only power-of-two bit depths are supported (8 for GDT_Bye, 16 for - * GDT_Int16/GDT_UInt16, 32 for GDT_Float32 and 64 for GDT_Float64) + * GDT_Int16/GDT_UInt16/GDT_Float16, 32 for GDT_Float32 and 64 for GDT_Float64) * * The pointer returned remains valid until CPLVirtualMemFree() is called. * CPLVirtualMemFree() must be called before the raster band object is diff --git a/gcore/overview.cpp b/gcore/overview.cpp index 4aa10da83ff8..8d1c4f3b2444 100644 --- a/gcore/overview.cpp +++ b/gcore/overview.cpp @@ -30,6 +30,7 @@ #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_progress.h" #include "cpl_vsi.h" #include "gdal.h" @@ -173,6 +174,7 @@ static CPLErr GDALResampleChunk_Near(const GDALOverviewResampleArgs &args, case GDT_Int16: case GDT_UInt16: + case GDT_Float16: { CPLAssert(GDALGetDataTypeSizeBytes(args.eWrkDataType) == 2); return GDALResampleChunk_NearT( @@ -181,6 +183,7 @@ static CPLErr GDALResampleChunk_Near(const GDALOverviewResampleArgs &args, } case GDT_CInt16: + case GDT_CFloat16: case GDT_Int32: case GDT_UInt32: case GDT_Float32: @@ -2422,6 +2425,7 @@ static CPLErr GDALResampleChunk_Mode(const GDALOverviewResampleArgs &args, case GDT_Int16: case GDT_UInt16: + case GDT_Float16: { CPLAssert(GDALGetDataTypeSizeBytes(args.eWrkDataType) == 2); return GDALResampleChunk_ModeT( @@ -2430,6 +2434,7 @@ static CPLErr GDALResampleChunk_Mode(const GDALOverviewResampleArgs &args, } case GDT_CInt16: + case GDT_CFloat16: case GDT_Int32: case GDT_UInt32: { @@ -3141,8 +3146,8 @@ static CPLErr GDALResampleChunk_ConvolutionT( constexpr int nWrkDataTypeSize = static_cast(sizeof(Twork)); // TODO: we should have some generic function to do this. - Twork fDstMin = -std::numeric_limits::max(); - Twork fDstMax = std::numeric_limits::max(); + Twork fDstMin = cpl::NumericLimits::lowest(); + Twork fDstMax = cpl::NumericLimits::max(); if (dstDataType == GDT_Byte) { fDstMin = std::numeric_limits::min(); diff --git a/gcore/rasterio.cpp b/gcore/rasterio.cpp index bc06c11a8d7c..957f9890f936 100644 --- a/gcore/rasterio.cpp +++ b/gcore/rasterio.cpp @@ -32,6 +32,7 @@ #include "cpl_conv.h" #include "cpl_cpu_features.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_progress.h" #include "cpl_string.h" #include "cpl_vsi.h" @@ -331,7 +332,6 @@ CPLErr GDALRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, else { // Type to type conversion. - if (eRWFlag == GF_Read) GDALCopyWords64( pabySrcBlock + nSrcByteOffset, eDataType, nBandDataSize, @@ -715,7 +715,6 @@ CPLErr GDALRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, { /* type to type conversion ... ouch, this is expensive way of handling single words */ - GDALCopyWords64(static_cast(pData) + iBufOffset, eBufType, 0, pabyDstBlock + iDstOffset, eDataType, 0, 1); @@ -2824,6 +2823,11 @@ inline void GDALCopyWordsFromT(const T *const CPL_RESTRICT pSrcData, static_cast(pDstData), nDstPixelStride, nWordCount); break; + case GDT_Float16: + GDALCopyWordsT(pSrcData, nSrcPixelStride, + static_cast(pDstData), nDstPixelStride, + nWordCount); + break; case GDT_Float32: GDALCopyWordsT(pSrcData, nSrcPixelStride, static_cast(pDstData), nDstPixelStride, @@ -2864,6 +2868,21 @@ inline void GDALCopyWordsFromT(const T *const CPL_RESTRICT pSrcData, nDstPixelStride, nWordCount); } break; + case GDT_CFloat16: + if (bInComplex) + { + GDALCopyWordsComplexT(pSrcData, nSrcPixelStride, + static_cast(pDstData), + nDstPixelStride, nWordCount); + } + else // input is not complex, so we need to promote to a complex + // buffer + { + GDALCopyWordsComplexOutT(pSrcData, nSrcPixelStride, + static_cast(pDstData), + nDstPixelStride, nWordCount); + } + break; case GDT_CFloat32: if (bInComplex) { @@ -3003,6 +3022,7 @@ static void GDALReplicateWord(const void *CPL_RESTRICT pSrcData, CASE_DUPLICATE_SIMPLE(GDT_Int32, GInt32) CASE_DUPLICATE_SIMPLE(GDT_UInt64, std::uint64_t) CASE_DUPLICATE_SIMPLE(GDT_Int64, std::int64_t) + CASE_DUPLICATE_SIMPLE(GDT_Float16, GFloat16) CASE_DUPLICATE_SIMPLE(GDT_Float32, float) CASE_DUPLICATE_SIMPLE(GDT_Float64, double) @@ -3023,6 +3043,7 @@ static void GDALReplicateWord(const void *CPL_RESTRICT pSrcData, CASE_DUPLICATE_COMPLEX(GDT_CInt16, GInt16) CASE_DUPLICATE_COMPLEX(GDT_CInt32, GInt32) + CASE_DUPLICATE_COMPLEX(GDT_CFloat16, GFloat16) CASE_DUPLICATE_COMPLEX(GDT_CFloat32, float) CASE_DUPLICATE_COMPLEX(GDT_CFloat64, double) @@ -3518,6 +3539,11 @@ void CPL_STDCALL GDALCopyWords64(const void *CPL_RESTRICT pSrcData, static_cast(pSrcData), nSrcPixelStride, false, pDstData, eDstType, nDstPixelStride, nWordCount); break; + case GDT_Float16: + GDALCopyWordsFromT( + static_cast(pSrcData), nSrcPixelStride, false, + pDstData, eDstType, nDstPixelStride, nWordCount); + break; case GDT_Float32: GDALCopyWordsFromT(static_cast(pSrcData), nSrcPixelStride, false, pDstData, @@ -3538,6 +3564,11 @@ void CPL_STDCALL GDALCopyWords64(const void *CPL_RESTRICT pSrcData, nSrcPixelStride, true, pDstData, eDstType, nDstPixelStride, nWordCount); break; + case GDT_CFloat16: + GDALCopyWordsFromT( + static_cast(pSrcData), nSrcPixelStride, true, + pDstData, eDstType, nDstPixelStride, nWordCount); + break; case GDT_CFloat32: GDALCopyWordsFromT(static_cast(pSrcData), nSrcPixelStride, true, pDstData, eDstType, @@ -5880,10 +5911,12 @@ static void GDALTranspose2D(const void *pSrc, GDALDataType eSrcType, DST *pDst, case GDT_Int32: CALL_GDALTranspose2D_internal(int32_t); break; case GDT_UInt64: CALL_GDALTranspose2D_internal(uint64_t); break; case GDT_Int64: CALL_GDALTranspose2D_internal(int64_t); break; + case GDT_Float16: CALL_GDALTranspose2D_internal(GFloat16); break; case GDT_Float32: CALL_GDALTranspose2D_internal(float); break; case GDT_Float64: CALL_GDALTranspose2D_internal(double); break; case GDT_CInt16: CALL_GDALTranspose2DComplex_internal(int16_t); break; case GDT_CInt32: CALL_GDALTranspose2DComplex_internal(int32_t); break; + case GDT_CFloat16: CALL_GDALTranspose2DComplex_internal(GFloat16); break; case GDT_CFloat32: CALL_GDALTranspose2DComplex_internal(float); break; case GDT_CFloat64: CALL_GDALTranspose2DComplex_internal(double); break; case GDT_Unknown: @@ -6140,10 +6173,12 @@ void GDALTranspose2D(const void *pSrc, GDALDataType eSrcType, void *pDst, case GDT_Int32: CALL_GDALTranspose2D_internal(int32_t, false); break; case GDT_UInt64: CALL_GDALTranspose2D_internal(uint64_t, false); break; case GDT_Int64: CALL_GDALTranspose2D_internal(int64_t, false); break; + case GDT_Float16: CALL_GDALTranspose2D_internal(GFloat16, false); break; case GDT_Float32: CALL_GDALTranspose2D_internal(float, false); break; case GDT_Float64: CALL_GDALTranspose2D_internal(double, false); break; case GDT_CInt16: CALL_GDALTranspose2D_internal(int16_t, true); break; case GDT_CInt32: CALL_GDALTranspose2D_internal(int32_t, true); break; + case GDT_CFloat16: CALL_GDALTranspose2D_internal(GFloat16, true); break; case GDT_CFloat32: CALL_GDALTranspose2D_internal(float, true); break; case GDT_CFloat64: CALL_GDALTranspose2D_internal(double, true); break; case GDT_Unknown: diff --git a/gcore/rawdataset.cpp b/gcore/rawdataset.cpp index 22bc2693fbf5..b5e62592bdff 100644 --- a/gcore/rawdataset.cpp +++ b/gcore/rawdataset.cpp @@ -458,6 +458,11 @@ void RawRasterBand::DoByteSwap(void *pBuffer, size_t nValues, int nByteSkip, nValues, nByteSkip); } } + else if (eDataType == GDT_Float16 || eDataType == GDT_CFloat16) + { + // No VAX support for GFloat16 + std::abort(); + } else if (eDataType == GDT_Float32 || eDataType == GDT_CFloat32) { GByte *pPtr = static_cast(pBuffer); diff --git a/ogr/ogrsf_frmts/gpkg/gdalgeopackagerasterband.cpp b/ogr/ogrsf_frmts/gpkg/gdalgeopackagerasterband.cpp index c29fb8b9639d..e4886a9457aa 100644 --- a/ogr/ogrsf_frmts/gpkg/gdalgeopackagerasterband.cpp +++ b/ogr/ogrsf_frmts/gpkg/gdalgeopackagerasterband.cpp @@ -15,6 +15,7 @@ #include "gdal_alg_priv.h" #include "ogrsqlitevfs.h" #include "cpl_error.h" +#include "cpl_float.h" #include #include diff --git a/port/CMakeLists.txt b/port/CMakeLists.txt index 3a66abb533f1..406a6cdd6d69 100644 --- a/port/CMakeLists.txt +++ b/port/CMakeLists.txt @@ -14,6 +14,7 @@ target_public_header( cpl_conv.h cpl_csv.h cpl_error.h + cpl_float.h cpl_hash_set.h cpl_http.h cpl_json.h diff --git a/port/cpl_float.cpp b/port/cpl_float.cpp index d92dde8eb0bb..7fec7d49ae3e 100644 --- a/port/cpl_float.cpp +++ b/port/cpl_float.cpp @@ -46,6 +46,8 @@ #include "cpl_float.h" #include "cpl_error.h" +#include + /************************************************************************/ /* HalfToFloat() */ /* */ @@ -288,3 +290,19 @@ GUInt16 CPLFloatToHalf(GUInt32 iFloat32, bool &bHasWarned) // coverity[overflow_sink] return static_cast((iSign << 15) | (iExponent << 10) | iMantissa); } + +GUInt16 CPLConvertFloatToHalf(float fFloat32) +{ + GUInt32 nFloat32; + std::memcpy(&nFloat32, &fFloat32, sizeof nFloat32); + bool bHasWarned = true; + return CPLFloatToHalf(nFloat32, bHasWarned); +} + +float CPLConvertHalfToFloat(GUInt16 nHalf) +{ + GUInt32 nFloat32 = CPLHalfToFloat(nHalf); + float fFloat32; + std::memcpy(&fFloat32, &nFloat32, sizeof fFloat32); + return fFloat32; +} diff --git a/port/cpl_float.h b/port/cpl_float.h index a1277a7c06c7..b202bd4f576c 100644 --- a/port/cpl_float.h +++ b/port/cpl_float.h @@ -49,13 +49,648 @@ #include "cpl_port.h" +#ifdef __cplusplus +#include +#include +#include +#include +#include +#ifdef HAVE_STD_FLOAT16_T +#include +#endif +#endif + CPL_C_START GUInt32 CPL_DLL CPLHalfToFloat(GUInt16 iHalf); GUInt32 CPL_DLL CPLTripleToFloat(GUInt32 iTriple); CPL_C_END #ifdef __cplusplus + GUInt16 CPL_DLL CPLFloatToHalf(GUInt32 iFloat32, bool &bHasWarned); + +GUInt16 CPL_DLL CPLConvertFloatToHalf(float fFloat32); +float CPL_DLL CPLConvertHalfToFloat(GUInt16 nHalf); + +namespace cpl +{ + +// We define our own version of `std::numeric_limits` so that we can +// specialized it for `cpl::Float16` if necessary. Specializing +// `std::numeric_limits` doesn't always work because some libraries +// use `std::numeric_limits`, and one cannot specialize a type +// template after it has been used. +template struct NumericLimits : std::numeric_limits +{ +}; + +#ifndef HAVE_STD_FLOAT16_T + +// Define a type `cpl::Float16`. If the compiler supports it natively +// (as `_Float16`), then this class is a simple wrapper. Otherwise we +// store the values in a `GUInt16` as bit pattern. + +//! @cond Doxygen_Suppress +struct Float16 +{ + struct make_from_bits_and_value + { + }; + +#ifdef HAVE__FLOAT16 + + // How we represent a `Float16` internally + using repr = _Float16; + + // How we compute on `Float16` values + using compute = _Float16; + + // Create a Float16 in a constexpr manner. Since we can't convert + // bits in a constexpr function, we need to take both the bit + // pattern and a float value as input, and can then choose which + // of the two to use. + constexpr Float16(make_from_bits_and_value, CPL_UNUSED std::uint16_t bits, + float fValue) + : rValue(repr(fValue)) + { + } + + static constexpr repr computeToRepr(compute fValue) + { + return fValue; + } + + static constexpr compute reprToCompute(repr rValue) + { + return rValue; + } + + template static constexpr repr toRepr(T fValue) + { + return static_cast(fValue); + } + + template static constexpr T fromRepr(repr rValue) + { + return static_cast(rValue); + } + +#else // #ifndef HAVE__FLOAT16 + + // How we represent a `Float16` internally + using repr = std::uint16_t; + + // How we compute on `Float16` values + using compute = float; + + // Create a Float16 in a constexpr manner. Since we can't convert + // bits in a constexpr function, we need to take both the bit + // pattern and a float value as input, and can then choose which + // of the two to use. + constexpr Float16(make_from_bits_and_value, std::uint16_t bits, + CPL_UNUSED float fValue) + : rValue(bits) + { + } + + static unsigned float2unsigned(float f) + { + unsigned u; + std::memcpy(&u, &f, 4); + return u; + } + + static float unsigned2float(unsigned u) + { + float f; + std::memcpy(&f, &u, 4); + return f; + } + + // Copied from cpl_float.cpp so that we can inline for performance + static std::uint16_t computeToRepr(float fFloat32) + { + std::uint32_t iFloat32 = float2unsigned(fFloat32); + + std::uint32_t iSign = (iFloat32 >> 31) & 0x00000001; + std::uint32_t iExponent = (iFloat32 >> 23) & 0x000000ff; + std::uint32_t iMantissa = iFloat32 & 0x007fffff; + + if (iExponent == 255) + { + if (iMantissa == 0) + { + // Positive or negative infinity. + return static_cast((iSign << 15) | 0x7C00); + } + + // NaN -- preserve sign and significand bits. + if (iMantissa >> 13) + return static_cast((iSign << 15) | 0x7C00 | + (iMantissa >> 13)); + return static_cast((iSign << 15) | 0x7E00); + } + + if (iExponent <= 127 - 15) + { + // Zero, float32 denormalized number or float32 too small normalized + // number + if (13 + 1 + 127 - 15 - iExponent >= 32) + return static_cast(iSign << 15); + + // Return a denormalized number + return static_cast( + (iSign << 15) | + ((iMantissa | 0x00800000) >> (13 + 1 + 127 - 15 - iExponent))); + } + + if (iExponent - (127 - 15) >= 31) + { + return static_cast((iSign << 15) | + 0x7C00); // Infinity + } + + // Normalized number. + iExponent = iExponent - (127 - 15); + iMantissa = iMantissa >> 13; + + // Assemble sign, exponent and mantissa. + // coverity[overflow_sink] + return static_cast((iSign << 15) | (iExponent << 10) | + iMantissa); + } + + // Copied from cpl_float.cpp so that we can inline for performance + static float reprToCompute(std::uint16_t iHalf) + { + std::uint32_t iSign = (iHalf >> 15) & 0x00000001; + int iExponent = (iHalf >> 10) & 0x0000001f; + std::uint32_t iMantissa = iHalf & 0x000003ff; + + if (iExponent == 31) + { + if (iMantissa == 0) + { + // Positive or negative infinity. + return unsigned2float((iSign << 31) | 0x7f800000); + } + + // NaN -- preserve sign and significand bits. + return unsigned2float((iSign << 31) | 0x7f800000 | + (iMantissa << 13)); + } + + if (iExponent == 0) + { + if (iMantissa == 0) + { + // Plus or minus zero. + return unsigned2float(iSign << 31); + } + + // Denormalized number -- renormalize it. + while (!(iMantissa & 0x00000400)) + { + iMantissa <<= 1; + iExponent -= 1; + } + + iExponent += 1; + iMantissa &= ~0x00000400U; + } + + // Normalized number. + iExponent = iExponent + (127 - 15); + iMantissa = iMantissa << 13; + + // Assemble sign, exponent and mantissa. + /* coverity[overflow_sink] */ + return unsigned2float((iSign << 31) | + (static_cast(iExponent) << 23) | + iMantissa); + } + + template static repr toRepr(T fValue) + { + return computeToRepr(static_cast(fValue)); + } + + template static T fromRepr(repr rValue) + { + return static_cast(reprToCompute(rValue)); + } + +#endif // #ifndef HAVE__FLOAT16 + + private: + repr rValue; + + public: + compute get() const + { + return reprToCompute(rValue); + } + + Float16() = default; + Float16(const Float16 &) = default; + Float16(Float16 &&) = default; + Float16 &operator=(const Float16 &) = default; + Float16 &operator=(Float16 &&) = default; + + // Constructors and conversion operators + +#ifdef HAVE__FLOAT16 + // cppcheck-suppress noExplicitConstructor + constexpr Float16(_Float16 hfValue) : rValue(hfValue) + { + } + + constexpr operator _Float16() const + { + return rValue; + } +#endif + + // cppcheck-suppress-macro noExplicitConstructor +#define GDAL_DEFINE_CONVERSION(TYPE) \ + \ + Float16(TYPE fValue) : rValue(toRepr(fValue)) \ + { \ + } \ + \ + operator TYPE() const \ + { \ + return fromRepr(rValue); \ + } + + GDAL_DEFINE_CONVERSION(float) + GDAL_DEFINE_CONVERSION(double) + GDAL_DEFINE_CONVERSION(char) + GDAL_DEFINE_CONVERSION(signed char) + GDAL_DEFINE_CONVERSION(short) + GDAL_DEFINE_CONVERSION(int) + GDAL_DEFINE_CONVERSION(long) + GDAL_DEFINE_CONVERSION(long long) + GDAL_DEFINE_CONVERSION(unsigned char) + GDAL_DEFINE_CONVERSION(unsigned short) + GDAL_DEFINE_CONVERSION(unsigned int) + GDAL_DEFINE_CONVERSION(unsigned long) + GDAL_DEFINE_CONVERSION(unsigned long long) + +#undef GDAL_DEFINE_CONVERSION + + // Arithmetic operators + + friend Float16 operator+(Float16 x) + { + return +x.get(); + } + + friend Float16 operator-(Float16 x) + { + return -x.get(); + } + +#define GDAL_DEFINE_ARITHOP(OP) \ + \ + friend Float16 operator OP(Float16 x, Float16 y) \ + { \ + return x.get() OP y.get(); \ + } \ + \ + friend double operator OP(double x, Float16 y) \ + { \ + return x OP y.get(); \ + } \ + \ + friend float operator OP(float x, Float16 y) \ + { \ + return x OP y.get(); \ + } \ + \ + friend Float16 operator OP(int x, Float16 y) \ + { \ + return x OP y.get(); \ + } \ + \ + friend double operator OP(Float16 x, double y) \ + { \ + return x.get() OP y; \ + } \ + \ + friend float operator OP(Float16 x, float y) \ + { \ + return x.get() OP y; \ + } \ + \ + friend Float16 operator OP(Float16 x, int y) \ + { \ + return x.get() OP y; \ + } + + GDAL_DEFINE_ARITHOP(+) + GDAL_DEFINE_ARITHOP(-) + GDAL_DEFINE_ARITHOP(*) + GDAL_DEFINE_ARITHOP(/) + +#undef GDAL_DEFINE_ARITHOP + + // Comparison operators + +#define GDAL_DEFINE_COMPARISON(OP) \ + \ + friend bool operator OP(Float16 x, Float16 y) \ + { \ + return x.get() OP y.get(); \ + } \ + \ + friend bool operator OP(float x, Float16 y) \ + { \ + return x OP y.get(); \ + } \ + \ + friend bool operator OP(double x, Float16 y) \ + { \ + return x OP y.get(); \ + } \ + \ + friend bool operator OP(int x, Float16 y) \ + { \ + return x OP y.get(); \ + } \ + \ + friend bool operator OP(Float16 x, float y) \ + { \ + return x.get() OP y; \ + } \ + \ + friend bool operator OP(Float16 x, double y) \ + { \ + return x.get() OP y; \ + } \ + \ + friend bool operator OP(Float16 x, int y) \ + { \ + return x.get() OP y; \ + } + + GDAL_DEFINE_COMPARISON(==) + GDAL_DEFINE_COMPARISON(!=) + GDAL_DEFINE_COMPARISON(<) + GDAL_DEFINE_COMPARISON(>) + GDAL_DEFINE_COMPARISON(<=) + GDAL_DEFINE_COMPARISON(>=) + +#undef GDAL_DEFINE_COMPARISON + + // Standard math functions + + friend bool isfinite(Float16 x) + { + using std::isfinite; + return isfinite(float(x)); + } + + friend bool isinf(Float16 x) + { + using std::isinf; + return isinf(float(x)); + } + + friend bool isnan(Float16 x) + { + using std::isnan; + return isnan(float(x)); + } + + friend bool isnormal(Float16 x) + { + using std::isnormal; + return isnormal(float(x)); + } + + friend bool signbit(Float16 x) + { + using std::signbit; + return signbit(float(x)); + } + + friend Float16 abs(Float16 x) + { + using std::abs; + return Float16(abs(float(x))); + } + + friend Float16 cbrt(Float16 x) + { + using std::cbrt; + return Float16(cbrt(float(x))); + } + + friend Float16 ceil(Float16 x) + { + using std::ceil; + return Float16(ceil(float(x))); + } + + friend Float16 copysign(Float16 x, Float16 y) + { + using std::copysign; + return Float16(copysign(float(x), float(y))); + } + + friend Float16 fabs(Float16 x) + { + using std::fabs; + return Float16(fabs(float(x))); + } + + friend Float16 floor(Float16 x) + { + using std::floor; + return Float16(floor(float(x))); + } + + friend Float16 fmax(Float16 x, Float16 y) + { + using std::fmax; + return Float16(fmax(float(x), float(y))); + } + + friend Float16 fmin(Float16 x, Float16 y) + { + using std::fmin; + return Float16(fmin(float(x), float(y))); + } + + friend Float16 hypot(Float16 x, Float16 y) + { + using std::hypot; + return Float16(hypot(float(x), float(y))); + } + + friend Float16 max(Float16 x, Float16 y) + { + using std::max; + return Float16(max(float(x), float(y))); + } + + friend Float16 min(Float16 x, Float16 y) + { + using std::min; + return Float16(min(float(x), float(y))); + } + + // Adapted from the LLVM Project, under the Apache License v2.0 + friend Float16 nextafter(Float16 x, Float16 y) + { + if (isnan(x)) + return x; + if (isnan(y)) + return y; + if (x == y) + return y; + + std::uint16_t bits; + if (x != Float16(0)) + { + std::memcpy(&bits, &x.rValue, 2); + if ((x < y) == (x > Float16(0))) + ++bits; + else + --bits; + } + else + { + bits = (signbit(y) << 15) | 0x0001; + } + + Float16 r; + std::memcpy(&r.rValue, &bits, 2); + + return r; + } + + friend Float16 pow(Float16 x, Float16 y) + { + using std::pow; + return Float16(pow(float(x), float(y))); + } + + friend Float16 pow(Float16 x, int n) + { + using std::pow; + return Float16(pow(float(x), n)); + } + + friend Float16 round(Float16 x) + { + using std::round; + return Float16(round(float(x))); + } + + friend Float16 sqrt(Float16 x) + { + using std::sqrt; + return Float16(sqrt(float(x))); + } +}; + +template <> struct NumericLimits +{ + static constexpr bool is_specialized = true; + static constexpr bool is_signed = true; + static constexpr bool is_integer = false; + static constexpr bool is_exact = false; + static constexpr bool has_infinity = true; + static constexpr bool has_quiet_NaN = true; + static constexpr bool has_signaling_NaN = true; + static constexpr bool has_denorm = true; + static constexpr bool is_iec559 = true; + + static constexpr int digits = 11; + static constexpr int digits10 = 3; + static constexpr int max_digits10 = 5; + static constexpr int radix = 2; + + static constexpr Float16 epsilon() + { + return Float16(Float16::make_from_bits_and_value{}, 0x1400, 0.000977f); + } + + static constexpr Float16 min() + { + return Float16(Float16::make_from_bits_and_value{}, 0x0001, 6.0e-8f); + } + + static constexpr Float16 lowest() + { + return Float16(Float16::make_from_bits_and_value{}, 0xfbff, -65504.0f); + } + + static constexpr Float16 max() + { + return Float16(Float16::make_from_bits_and_value{}, 0x7bff, +65504.0f); + } + + static constexpr Float16 infinity() + { + return Float16(Float16::make_from_bits_and_value{}, 0x7c00, + std::numeric_limits::infinity()); + } + + static constexpr Float16 quiet_NaN() + { + return Float16(Float16::make_from_bits_and_value{}, 0x7e00, + std::numeric_limits::quiet_NaN()); + } + + static constexpr Float16 signaling_NaN() + { + return Float16(Float16::make_from_bits_and_value{}, 0xfe00, + std::numeric_limits::signaling_NaN()); + } +}; + +//! @endcond + +#endif // #ifndef HAVE_STD_FLOAT16_T + +} // namespace cpl + +#ifdef HAVE_STD_FLOAT16_T +using GFloat16 = std::float16_t; +#else +using GFloat16 = cpl::Float16; #endif +// Define some GDAL wrappers. Their C equivalents are defined in `cpl_port.h`. +// (These wrappers are not necessary any mroe in C++, one can always +// call `isnan` etc directly.) + +template constexpr int CPLIsNan(T x) +{ + // We need to write `using std::isnan` instead of directly using + // `std::isnan` because `std::isnan` only supports the types + // `float` and `double`. The `isnan` for `cpl::Float16` is found in the + // `cpl` namespace via argument-dependent lookup + // . + using std::isnan; + return isnan(x); +} + +template constexpr int CPLIsInf(T x) +{ + using std::isinf; + return isinf(x); +} + +template constexpr int CPLIsFinite(T x) +{ + using std::isfinite; + return isfinite(x); +} + +#endif // #ifdef __cplusplus + #endif // CPL_FLOAT_H_INCLUDED diff --git a/port/cpl_json_streaming_writer.cpp b/port/cpl_json_streaming_writer.cpp index 9888ccd0540c..897640d8752c 100644 --- a/port/cpl_json_streaming_writer.cpp +++ b/port/cpl_json_streaming_writer.cpp @@ -17,6 +17,7 @@ #include #include "cpl_conv.h" +#include "cpl_float.h" #include "cpl_string.h" #include "cpl_json_streaming_writer.h" @@ -222,6 +223,25 @@ void CPLJSonStreamingWriter::Add(std::uint64_t nVal) Print(CPLSPrintf(CPL_FRMT_GUIB, static_cast(nVal))); } +void CPLJSonStreamingWriter::Add(GFloat16 hfVal, int nPrecision) +{ + EmitCommaIfNeeded(); + if (CPLIsNan(hfVal)) + { + Print("\"NaN\""); + } + else if (CPLIsInf(hfVal)) + { + Print(hfVal > 0 ? "\"Infinity\"" : "\"-Infinity\""); + } + else + { + char szFormatting[10]; + snprintf(szFormatting, sizeof(szFormatting), "%%.%dg", nPrecision); + Print(CPLSPrintf(szFormatting, float(hfVal))); + } +} + void CPLJSonStreamingWriter::Add(float fVal, int nPrecision) { EmitCommaIfNeeded(); diff --git a/port/cpl_json_streaming_writer.h b/port/cpl_json_streaming_writer.h index 03e9307eb4c0..56f4e6b5ac65 100644 --- a/port/cpl_json_streaming_writer.h +++ b/port/cpl_json_streaming_writer.h @@ -20,6 +20,8 @@ #include #include #include + +#include "cpl_float.h" #include "cpl_port.h" class CPL_DLL CPLJSonStreamingWriter @@ -93,8 +95,9 @@ class CPL_DLL CPLJSonStreamingWriter void Add(std::int64_t nVal); void Add(std::uint64_t nVal); + void Add(GFloat16 hfVal, int nPrecision = 5); void Add(float fVal, int nPrecision = 9); - void Add(double dfVal, int nPrecision = 18); + void Add(double dfVal, int nPrecision = 17); void AddNull(); void StartObj(); diff --git a/port/cpl_port.h b/port/cpl_port.h index 960d40befec8..0ed92aef1626 100644 --- a/port/cpl_port.h +++ b/port/cpl_port.h @@ -123,6 +123,13 @@ #include #endif +#ifdef __cplusplus +extern "C++" +{ +#include +} +#endif + /* ==================================================================== */ /* Base portability stuff ... this stuff may need to be */ /* modified for new platforms. */ @@ -551,8 +558,8 @@ static inline char *CPL_afl_friendly_strstr(const char *haystack, #endif /*! @endcond */ -#ifndef GDAL_COMPILATION /*! @cond Doxygen_Suppress */ +#ifndef __cplusplus /* -------------------------------------------------------------------- */ /* Handle isnan() and isinf(). Note that isinf() and isnan() */ /* are supposed to be macros according to C99, defined in math.h */ @@ -574,99 +581,25 @@ static inline char *CPL_afl_friendly_strstr(const char *haystack, #define CPLIsNan(x) __builtin_isnan(x) #define CPLIsInf(x) __builtin_isinf(x) #define CPLIsFinite(x) __builtin_isfinite(x) -#elif defined(__cplusplus) && defined(HAVE_STD_IS_NAN) && HAVE_STD_IS_NAN -extern "C++" -{ -#ifndef DOXYGEN_SKIP -#include -#endif - static inline int CPLIsNan(float f) - { - return std::isnan(f); - } - - static inline int CPLIsNan(double f) - { - return std::isnan(f); - } - - static inline int CPLIsInf(float f) - { - return std::isinf(f); - } - - static inline int CPLIsInf(double f) - { - return std::isinf(f); - } - - static inline int CPLIsFinite(float f) - { - return std::isfinite(f); - } - - static inline int CPLIsFinite(double f) - { - return std::isfinite(f); - } -} -#else -/** Return whether a floating-pointer number is NaN */ -#if defined(__cplusplus) && defined(__GNUC__) && defined(__linux) && \ - !defined(__ANDROID__) && !defined(CPL_SUPRESS_CPLUSPLUS) -/* so to not get warning about conversion from double to float with */ -/* gcc -Wfloat-conversion when using isnan()/isinf() macros */ -extern "C++" -{ - static inline int CPLIsNan(float f) - { - return __isnanf(f); - } - - static inline int CPLIsNan(double f) - { - return __isnan(f); - } - - static inline int CPLIsInf(float f) - { - return __isinff(f); - } - - static inline int CPLIsInf(double f) - { - return __isinf(f); - } - - static inline int CPLIsFinite(float f) - { - return !__isnanf(f) && !__isinff(f); - } - - static inline int CPLIsFinite(double f) - { - return !__isnan(f) && !__isinf(f); - } -} -#else +#elif defined(isinf) || defined(__FreeBSD__) +/** Return whether a floating-pointer number is nan */ #define CPLIsNan(x) isnan(x) -#if defined(isinf) || defined(__FreeBSD__) /** Return whether a floating-pointer number is +/- infinity */ #define CPLIsInf(x) isinf(x) /** Return whether a floating-pointer number is finite */ #define CPLIsFinite(x) (!isnan(x) && !isinf(x)) #elif defined(__sun__) #include +#define CPLIsNan(x) isnan(x) #define CPLIsInf(x) (!finite(x) && !isnan(x)) #define CPLIsFinite(x) finite(x) #else +#define CPLIsNan(x) ((x) != (x)) #define CPLIsInf(x) (0) #define CPLIsFinite(x) (!isnan(x)) #endif #endif -#endif /*! @endcond */ -#endif // GDAL_COMPILATION /*! @cond Doxygen_Suppress */ /*--------------------------------------------------------------------- diff --git a/scripts/cppcheck.sh b/scripts/cppcheck.sh index f3b9932c2077..d73f03318fd4 100755 --- a/scripts/cppcheck.sh +++ b/scripts/cppcheck.sh @@ -376,7 +376,7 @@ if grep "noConstructor" ${LOG_FILE} ; then ret_code=1 fi -if grep "noExplicitConstructor" ${LOG_FILE} ; then +if grep "noExplicitConstructor" ${LOG_FILE} | grep -v -e port/cpl_float.h ; then echo "noExplicitConstructor check failed" ret_code=1 fi @@ -413,7 +413,7 @@ if grep "functionStatic" ${LOG_FILE} | grep -v -e OGRSQLiteDataSource::OpenRaste ret_code=1 fi -if grep "knownConditionTrueFalse" ${LOG_FILE} ; then +if grep "knownConditionTrueFalse" ${LOG_FILE} | grep -v -e gcore/gdal_priv_templates.hpp ; then echo "knownConditionTrueFalse check failed" ret_code=1 fi diff --git a/swig/include/gdal.i b/swig/include/gdal.i index b9d0b9c25f23..273807be61ac 100644 --- a/swig/include/gdal.i +++ b/swig/include/gdal.i @@ -139,13 +139,15 @@ typedef enum { /*! Thirty two bit signed integer */ GDT_Int32 = 5, /*! 64 bit unsigned integer */ GDT_UInt64 = 12, /*! 64 bit signed integer */ GDT_Int64 = 13, + /*! Sixteen bit floating point */ GDT_Float16 = 15, /*! Thirty two bit floating point */ GDT_Float32 = 6, /*! Sixty four bit floating point */ GDT_Float64 = 7, /*! Complex Int16 */ GDT_CInt16 = 8, /*! Complex Int32 */ GDT_CInt32 = 9, + /*! Complex Float16 */ GDT_CFloat16 = 16, /*! Complex Float32 */ GDT_CFloat32 = 10, /*! Complex Float64 */ GDT_CFloat64 = 11, - GDT_TypeCount = 15 /* maximum type # + 1 */ + GDT_TypeCount = 17 /* maximum type # + 1 */ } GDALDataType; /*! Types of color interpretation for raster bands. */ diff --git a/swig/include/gdal_array.i b/swig/include/gdal_array.i index 5cefdda97231..ad836a44d564 100644 --- a/swig/include/gdal_array.i +++ b/swig/include/gdal_array.i @@ -343,12 +343,18 @@ static GDALDataType NumpyTypeToGDALType(PyArrayObject *psArray) case NPY_CFLOAT: return GDT_CFloat32; + // case NPY_CHALF + // return GDT_CFloat16; + case NPY_DOUBLE: return GDT_Float64; case NPY_FLOAT: return GDT_Float32; + case NPY_HALF: + return GDT_Float16; + case NPY_INT32: return GDT_Int32; @@ -997,6 +1003,9 @@ static bool AddNumpyArrayToDict(PyObject *dict, { 'e', NPY_FLOAT16, 2 }, { 'f', NPY_FLOAT32, 4 }, { 'g', NPY_FLOAT64, 8 }, + // { 'E', NPY_COMPLEX32, 4 }, + // { 'F', NPY_COMPLEX64, 8 }, + // { 'G', NPY_COMPLEX128, 16 }, }; const size_t nEltsInMapArrowTypeToNumpyType = sizeof(MapArrowTypeToNumpyType) / sizeof(MapArrowTypeToNumpyType[0]); @@ -2003,9 +2012,9 @@ PyObject* _RecordBatchAsNumpy(VoidPtrAsLong recordBatchPtr, GIntBig nLineSpace = virtualmem->nLineSpace; /* if bAuto == TRUE */ int numpytype; - if( datatype == GDT_CInt16 || datatype == GDT_CInt32 ) + if( datatype == GDT_CInt16 || datatype == GDT_CInt32 || datatype == GDT_CFloat16 ) { - PyErr_SetString( PyExc_RuntimeError, "GDT_CInt16 and GDT_CInt32 not supported for now" ); + PyErr_SetString( PyExc_RuntimeError, "GDT_CInt16, GDT_CInt32, and GDT_CFloat16 not supported for now" ); SWIG_fail; } @@ -2019,10 +2028,12 @@ PyObject* _RecordBatchAsNumpy(VoidPtrAsLong recordBatchPtr, case GDT_UInt32: numpytype = NPY_UINT32; break; case GDT_Int64: numpytype = NPY_INT64; break; case GDT_UInt64: numpytype = NPY_UINT64; break; + case GDT_Float16: numpytype = NPY_FLOAT16; break; case GDT_Float32: numpytype = NPY_FLOAT32; break; case GDT_Float64: numpytype = NPY_FLOAT64; break; //case GDT_CInt16: numpytype = NPY_INT16; break; //case GDT_CInt32: numpytype = NPY_INT32; break; + //case GDT_CFloat16: numpytype = NPY_CHALF; break; case GDT_CFloat32: numpytype = NPY_CFLOAT; break; case GDT_CFloat64: numpytype = NPY_CDOUBLE; break; default: numpytype = NPY_UBYTE; break; @@ -2365,11 +2376,13 @@ codes = {gdalconst.GDT_Byte: numpy.uint8, gdalconst.GDT_Int32: numpy.int32, gdalconst.GDT_UInt64: numpy.uint64, gdalconst.GDT_Int64: numpy.int64, + gdalconst.GDT_Float16: numpy.float16, gdalconst.GDT_Float32: numpy.float32, gdalconst.GDT_Float64: numpy.float64, gdalconst.GDT_CInt16: numpy.complex64, gdalconst.GDT_CInt32: numpy.complex64, - gdalconst.GDT_CFloat32: numpy.complex64, + gdalconst.GDT_CFloat16: numpy.complex64, + gdalconst.GDT_CFloat32: numpy.complex64, gdalconst.GDT_CFloat64: numpy.complex128} np_class_to_gdal_code = { v : k for k, v in codes.items() } diff --git a/swig/include/gdalconst.i b/swig/include/gdalconst.i index 3b0db9f0a946..dd5ca775e20f 100644 --- a/swig/include/gdalconst.i +++ b/swig/include/gdalconst.i @@ -43,10 +43,12 @@ %constant GDT_Int32 = GDT_Int32; %constant GDT_UInt64 = GDT_UInt64; %constant GDT_Int64 = GDT_Int64; +%constant GDT_Float16 = GDT_Float16; %constant GDT_Float32 = GDT_Float32; %constant GDT_Float64 = GDT_Float64; %constant GDT_CInt16 = GDT_CInt16; %constant GDT_CInt32 = GDT_CInt32; +%constant GDT_CFloat16 = GDT_CFloat16; %constant GDT_CFloat32 = GDT_CFloat32; %constant GDT_CFloat64 = GDT_CFloat64; %constant GDT_TypeCount = GDT_TypeCount; diff --git a/swig/include/java/gdal_java.i b/swig/include/java/gdal_java.i index 70637de4a2ba..8b61d37e0e39 100644 --- a/swig/include/java/gdal_java.i +++ b/swig/include/java/gdal_java.i @@ -307,12 +307,12 @@ import java.lang.Integer; } %enddef - DEFINE_READ_MDA_DATA(char, GDT_Byte) - DEFINE_READ_MDA_DATA(short, GDT_Int16) - DEFINE_READ_MDA_DATA(int, GDT_Int32) - DEFINE_READ_MDA_DATA(int64_t, GDT_Int64) - DEFINE_READ_MDA_DATA(float, GDT_Float32) - DEFINE_READ_MDA_DATA(double, GDT_Float64) + DEFINE_READ_MDA_DATA(char, GDT_Byte) + DEFINE_READ_MDA_DATA(short, GDT_Int16) + DEFINE_READ_MDA_DATA(int, GDT_Int32) + DEFINE_READ_MDA_DATA(int64_t, GDT_Int64) + DEFINE_READ_MDA_DATA(float, GDT_Float32) + DEFINE_READ_MDA_DATA(double, GDT_Float64) %define DEFINE_WRITE_MDA_DATA(ctype, buffer_type_code) %apply(int nList, GInt64 *pList) { (int starts, GInt64 *startsValues) }; @@ -358,12 +358,12 @@ import java.lang.Integer; } %enddef - DEFINE_WRITE_MDA_DATA(char , GDT_Byte) - DEFINE_WRITE_MDA_DATA(short , GDT_Int16) - DEFINE_WRITE_MDA_DATA(int , GDT_Int32) - DEFINE_WRITE_MDA_DATA(int64_t, GDT_Int64) - DEFINE_WRITE_MDA_DATA(float , GDT_Float32) - DEFINE_WRITE_MDA_DATA(double , GDT_Float64) + DEFINE_WRITE_MDA_DATA(char , GDT_Byte) + DEFINE_WRITE_MDA_DATA(short , GDT_Int16) + DEFINE_WRITE_MDA_DATA(int , GDT_Int32) + DEFINE_WRITE_MDA_DATA(int64_t , GDT_Int64) + DEFINE_WRITE_MDA_DATA(float , GDT_Float32) + DEFINE_WRITE_MDA_DATA(double , GDT_Float64) } /* extend */ diff --git a/swig/include/python/gdal_python.i b/swig/include/python/gdal_python.i index bfd4a429dd90..3421377200fb 100644 --- a/swig/include/python/gdal_python.i +++ b/swig/include/python/gdal_python.i @@ -34,6 +34,7 @@ static int getAlignment(GDALDataType ntype) return 1; case GDT_Int16: case GDT_UInt16: + case GDT_Float16: return 2; case GDT_Int32: case GDT_UInt32: @@ -44,6 +45,7 @@ static int getAlignment(GDALDataType ntype) case GDT_UInt64: return 8; case GDT_CInt16: + case GDT_CFloat16: return 2; case GDT_CInt32: case GDT_CFloat32: @@ -184,8 +186,12 @@ static void readraster_releasebuffer(CPLErr eErr, gdalconst.GDT_UInt16: ("%su2" % byteorders[sys.byteorder]), gdalconst.GDT_Int32: ("%si4" % byteorders[sys.byteorder]), gdalconst.GDT_UInt32: ("%su4" % byteorders[sys.byteorder]), + gdalconst.GDT_Int64: ("%si8" % byteorders[sys.byteorder]), + gdalconst.GDT_UInt64: ("%su8" % byteorders[sys.byteorder]), + gdalconst.GDT_Float16: ("%sf2" % byteorders[sys.byteorder]), gdalconst.GDT_Float32: ("%sf4" % byteorders[sys.byteorder]), gdalconst.GDT_Float64: ("%sf8" % byteorders[sys.byteorder]), + gdalconst.GDT_CFloat16: ("%sf2" % byteorders[sys.byteorder]), gdalconst.GDT_CFloat32: ("%sf4" % byteorders[sys.byteorder]), gdalconst.GDT_CFloat64: ("%sf8" % byteorders[sys.byteorder]), gdalconst.GDT_Byte: ("%st8" % byteorders[sys.byteorder]), @@ -1845,6 +1851,10 @@ def GetMDArrayNames(self, options = []) -> "list[str]": buffer_stride.reverse() if not buffer_datatype: buffer_datatype = self.GetDataType() + if buffer_datatype.GetClass() == GEDTC_NUMERIC and buffer_datatype.GetNumericDataType() == gdalconst.GDT_Float16: + buffer_datatype = ExtendedDataType.Create(GDT_Float32) + elif buffer_datatype.GetClass() == GEDTC_NUMERIC and buffer_datatype.GetNumericDataType() == gdalconst.GDT_CFloat16: + buffer_datatype = ExtendedDataType.Create(GDT_CFloat32) return _gdal.MDArray_Read(self, array_start_idx, count, array_step, buffer_stride, buffer_datatype) def ReadAsArray(self, @@ -1922,8 +1932,12 @@ def GetMDArrayNames(self, options = []) -> "list[str]": ('l', 4): GDT_Int32, ('q', 8): GDT_Int64, ('Q', 8): GDT_UInt64, + ('e', 2): GDT_Float16, ('f', 4): GDT_Float32, - ('d', 8): GDT_Float64 + ('d', 8): GDT_Float64, + # ('E', 2): GDT_CFloat16, + ('F', 4): GDT_CFloat32, + ('D', 8): GDT_CFloat64 } key = (buffer.typecode, buffer.itemsize) if key not in map_typecode_itemsize_to_gdal: diff --git a/swig/include/python/typemaps_python.i b/swig/include/python/typemaps_python.i index e8a6cbbe9f0b..724b04399b49 100644 --- a/swig/include/python/typemaps_python.i +++ b/swig/include/python/typemaps_python.i @@ -2579,6 +2579,11 @@ DecomposeSequenceOf4DCoordinates( PyObject *seq, int nCount, double *x, double * buf->format = (char*) "I"; buf->itemsize = 4; } + else if( *($3) == GDT_Float16 ) + { + buf->format = (char*) "f"; + buf->itemsize = 2; + } else if( *($3) == GDT_Float32 ) { buf->format = (char*) "f"; diff --git a/swig/python/gdal-utils/osgeo_utils/gdal_calc.py b/swig/python/gdal-utils/osgeo_utils/gdal_calc.py index 51bfd2cd7e4f..c89256b95f0e 100644 --- a/swig/python/gdal-utils/osgeo_utils/gdal_calc.py +++ b/swig/python/gdal-utils/osgeo_utils/gdal_calc.py @@ -55,10 +55,12 @@ gdal.GDT_Int32: -2147483647, gdal.GDT_UInt64: None, gdal.GDT_Int64: None, + gdal.GDT_Float16: 65504.0, gdal.GDT_Float32: 3.402823466e38, gdal.GDT_Float64: 1.7976931348623158e308, gdal.GDT_CInt16: None, gdal.GDT_CInt32: None, + gdal.GDT_CFloat16: None, gdal.GDT_CFloat32: None, gdal.GDT_CFloat64: None, } @@ -634,6 +636,11 @@ def Calc( elif not isinstance(myResult, numpy.ndarray): myResult = numpy.ones((nYValid, nXValid)) * myResult + # Convert float16 to float32 if necessary + # (While numpy probably supports float16, GDAL may not) + if myResult.dtype == "float16": + myResult = numpy.float32(myResult) + # write data block to the output file myOutB = myOut.GetRasterBand(bandNo) if gdal_array.BandWriteArray(myOutB, myResult, xoff=myX, yoff=myY) != 0: diff --git a/swig/python/gdal-utils/osgeo_utils/samples/fft.py b/swig/python/gdal-utils/osgeo_utils/samples/fft.py index 81db9fd6e84a..bee0c09eecc2 100644 --- a/swig/python/gdal-utils/osgeo_utils/samples/fft.py +++ b/swig/python/gdal-utils/osgeo_utils/samples/fft.py @@ -42,6 +42,8 @@ def ParseType(typ): return gdal.GDT_Int32 elif typ == "UInt32": return gdal.GDT_UInt32 + elif typ == "Float16": + return gdal.GDT_Float16 elif typ == "Float32": return gdal.GDT_Float32 elif typ == "Float64": @@ -50,6 +52,8 @@ def ParseType(typ): return gdal.GDT_CInt16 elif typ == "CInt32": return gdal.GDT_CInt32 + elif typ == "CFloat16": + return gdal.GDT_CFloat16 elif typ == "CFloat32": return gdal.GDT_CFloat32 elif typ == "CFloat64": diff --git a/swig/python/gdal-utils/osgeo_utils/samples/rel.py b/swig/python/gdal-utils/osgeo_utils/samples/rel.py index 75d01aacb197..ad91569f2a6f 100644 --- a/swig/python/gdal-utils/osgeo_utils/samples/rel.py +++ b/swig/python/gdal-utils/osgeo_utils/samples/rel.py @@ -66,6 +66,8 @@ def ParseType(typ): return gdal.GDT_Int32 if typ == "UInt32": return gdal.GDT_UInt32 + if typ == "Float16": + return gdal.GDT_Float16 if typ == "Float32": return gdal.GDT_Float32 if typ == "Float64": @@ -74,6 +76,8 @@ def ParseType(typ): return gdal.GDT_CInt16 if typ == "CInt32": return gdal.GDT_CInt32 + if typ == "CFloat16": + return gdal.GDT_CFloat16 if typ == "CFloat32": return gdal.GDT_CFloat32 if typ == "CFloat64":