Initial implementation of float16 support (RFC100) (#11180)

.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)
-

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 <cmath>
+    #include <stdfloat>
+    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

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.

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)

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<GUInt32>(nXSize) >
-        std::numeric_limits<T>::max() / static_cast<GUInt32>(nYSize))
+        cpl::NumericLimits<T>::max() / static_cast<GUInt32>(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<GUInt32>(nXSize) <=
-            std::numeric_limits<T>::max() / static_cast<GUInt32>(nYSize))
+            cpl::NumericLimits<T>::max() / static_cast<GUInt32>(nYSize))
     {
         nPixels = static_cast<T>(nXSize) * static_cast<T>(nYSize);
     }

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<T>::max();
+        nMaxValue = cpl::NumericLimits<T>::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<GFloat16 *>(pDataBuf), nValues,
+                           nBandValues, nMaxValue);
+            break;
+
         case GDT_Float32:
             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
                            static_cast<float *>(pDataBuf), nValues, nBandValues,
@@ -1090,6 +1097,12 @@ CPLErr GDALPansharpenOperation::WeightedBrovey(
                 static_cast<std::int64_t *>(pDataBuf), nValues, nBandValues, 0);
             break;
 
+        case GDT_Float16:
+            WeightedBrovey3<WorkDataType, GFloat16, FALSE>(
+                pPanBuffer, pUpsampledSpectralBuffer,
+                static_cast<GFloat16 *>(pDataBuf), nValues, nBandValues, 0);
+            break;
+
         case GDT_Float32:
             WeightedBrovey3<WorkDataType, float, FALSE>(
                 pPanBuffer, pUpsampledSpectralBuffer,
@@ -1769,6 +1782,13 @@ CPLErr GDALPansharpenOperation::PansharpenChunk(
                 pDataBuf, eBufDataType, nValues, nBandValues);
             break;
 
+        case GDT_Float16:
+            eErr = WeightedBrovey(
+                static_cast<const GFloat16 *>(pPanBuffer),
+                static_cast<const GFloat16 *>(pUpsampledSpectralBuffer),
+                pDataBuf, eBufDataType, nValues, nBandValues);
+            break;
+
         case GDT_Float32:
             eErr = WeightedBrovey(
                 static_cast<const float *>(pPanBuffer),

alg/gdalrasterize.cpp

@@ -244,6 +244,10 @@ static void gvBurnScanline(void *pCBData, int nY, int nXStart, int nXEnd,
             gvBurnScanlineBasic<std::uint64_t>(psInfo, nY, nXStart, nXEnd,
                                                dfVariant);
             break;
+        case GDT_Float16:
+            gvBurnScanlineBasic<GFloat16>(psInfo, nY, nXStart, nXEnd,
+                                          dfVariant);
+            break;
         case GDT_Float32:
             gvBurnScanlineBasic<float>(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<std::uint64_t>(psInfo, nY, nX, dfVariant);
             break;
+        case GDT_Float16:
+            gvBurnPointBasic<GFloat16>(psInfo, nY, nX, dfVariant);
+            break;
         case GDT_Float32:
             gvBurnPointBasic<float>(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:

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<T>::min() ||
-        padfNoData[0] > std::numeric_limits<T>::max() + 0.000001 ||
+    if (padfNoData[0] < cpl::NumericLimits<T>::min() ||
+        padfNoData[0] > cpl::NumericLimits<T>::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<float>::max()) == pszValueIn)
             {
-                return -std::numeric_limits<float>::max();
+                return std::numeric_limits<float>::lowest();
             }
             else if (eDataType == GDT_Float32 &&
                      CPLString().Printf("%.16g",