Initial implementation of float16 support

.github/workflows/icc/build.sh

@@ -2,10 +2,17 @@
 
 set -eu
 
+# 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 \
+    -DCMAKE_C_FLAGS=-DGDAL_DISABLE_FLOAT16 \
+    -DCMAKE_CXX_FLAGS=-DGDAL_DISABLE_FLOAT16 \
     -DUSE_CCACHE=ON
 make -j$(nproc)
 

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"

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"
@@ -1008,6 +1009,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,
@@ -1089,6 +1096,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,
@@ -1764,6 +1777,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"

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"
@@ -1669,6 +1670,11 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand,
                 dfDstImag = 0.0;
                 break;
 
+            case GDT_Float16:
+                dfDstReal = reinterpret_cast<GFloat16 *>(pabyDst)[iDstOffset];
+                dfDstImag = 0.0;
+                break;
+
             case GDT_Float32:
                 dfDstReal = reinterpret_cast<float *>(pabyDst)[iDstOffset];
                 dfDstImag = 0.0;
@@ -1691,6 +1697,13 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand,
                     reinterpret_cast<GInt32 *>(pabyDst)[iDstOffset * 2 + 1];
                 break;
 
+            case GDT_CFloat16:
+                dfDstReal =
+                    reinterpret_cast<GFloat16 *>(pabyDst)[iDstOffset * 2];
+                dfDstImag =
+                    reinterpret_cast<GFloat16 *>(pabyDst)[iDstOffset * 2 + 1];
+                break;
+
             case GDT_CFloat32:
                 dfDstReal = reinterpret_cast<float *>(pabyDst)[iDstOffset * 2];
                 dfDstImag =
@@ -1790,6 +1803,11 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand,
             CLAMP(std::int64_t);
             break;
 
+        case GDT_Float16:
+            reinterpret_cast<GFloat16 *>(pabyDst)[iDstOffset] =
+                static_cast<GFloat16>(dfReal);
+            break;
+
         case GDT_Float32:
             reinterpret_cast<float *>(pabyDst)[iDstOffset] =
                 static_cast<float>(dfReal);
@@ -1851,6 +1869,13 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand,
             break;
         }
 
+        case GDT_CFloat16:
+            reinterpret_cast<GFloat16 *>(pabyDst)[iDstOffset * 2] =
+                static_cast<GFloat16>(dfReal);
+            reinterpret_cast<GFloat16 *>(pabyDst)[iDstOffset * 2 + 1] =
+                static_cast<GFloat16>(dfImag);
+            break;
+
         case GDT_CFloat32:
             reinterpret_cast<float *>(pabyDst)[iDstOffset * 2] =
                 static_cast<float>(dfReal);
@@ -1943,6 +1968,10 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand,
                     reinterpret_cast<std::uint64_t *>(pabyDst)[iDstOffset]);
                 break;
 
+            case GDT_Float16:
+                dfDstReal = reinterpret_cast<GFloat16 *>(pabyDst)[iDstOffset];
+                break;
+
             case GDT_Float32:
                 dfDstReal = reinterpret_cast<float *>(pabyDst)[iDstOffset];
                 break;
@@ -1953,6 +1982,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:
@@ -2010,6 +2040,11 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand,
             CLAMP(std::int64_t);
             break;
 
+        case GDT_Float16:
+            reinterpret_cast<GFloat16 *>(pabyDst)[iDstOffset] =
+                static_cast<GFloat16>(dfReal);
+            break;
+
         case GDT_Float32:
             reinterpret_cast<float *>(pabyDst)[iDstOffset] =
                 static_cast<float>(dfReal);
@@ -2021,6 +2056,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;
@@ -2103,6 +2139,11 @@ static bool GWKGetPixelValue(const GDALWarpKernel *poWK, int iBand,
             *pdfImag = 0.0;
             break;
 
+        case GDT_Float16:
+            *pdfReal = reinterpret_cast<GFloat16 *>(pabySrc)[iSrcOffset];
+            *pdfImag = 0.0;
+            break;
+
         case GDT_Float32:
             *pdfReal = reinterpret_cast<float *>(pabySrc)[iSrcOffset];
             *pdfImag = 0.0;
@@ -2123,6 +2164,12 @@ static bool GWKGetPixelValue(const GDALWarpKernel *poWK, int iBand,
             *pdfImag = reinterpret_cast<GInt32 *>(pabySrc)[iSrcOffset * 2 + 1];
             break;
 
+        case GDT_CFloat16:
+            *pdfReal = reinterpret_cast<GFloat16 *>(pabySrc)[iSrcOffset * 2];
+            *pdfImag =
+                reinterpret_cast<GFloat16 *>(pabySrc)[iSrcOffset * 2 + 1];
+            break;
+
         case GDT_CFloat32:
             *pdfReal = reinterpret_cast<float *>(pabySrc)[iSrcOffset * 2];
             *pdfImag = reinterpret_cast<float *>(pabySrc)[iSrcOffset * 2 + 1];
@@ -2203,6 +2250,10 @@ static bool GWKGetPixelValueReal(const GDALWarpKernel *poWK, int iBand,
                 reinterpret_cast<std::uint64_t *>(pabySrc)[iSrcOffset]);
             break;
 
+        case GDT_Float16:
+            *pdfReal = reinterpret_cast<GFloat16 *>(pabySrc)[iSrcOffset];
+            break;
+
         case GDT_Float32:
             *pdfReal = reinterpret_cast<float *>(pabySrc)[iSrcOffset];
             break;
@@ -2213,6 +2264,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:
@@ -2409,6 +2461,19 @@ static bool GWKGetPixelRow(const GDALWarpKernel *poWK, int iBand,
             break;
         }
 
+        case GDT_Float16:
+        {
+            GFloat16 *pSrc =
+                reinterpret_cast<GFloat16 *>(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 =
@@ -2467,6 +2532,22 @@ static bool GWKGetPixelRow(const GDALWarpKernel *poWK, int iBand,
             break;
         }
 
+        case GDT_CFloat16:
+        {
+            GFloat16 *pSrc =
+                reinterpret_cast<GFloat16 *>(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 =

apps/gdaldem_lib.cpp

@@ -99,6 +99,7 @@
 #include <limits>
 
 #include "cpl_error.h"
+#include "cpl_float.h"
 #include "cpl_progress.h"
 #include "cpl_string.h"
 #include "cpl_vsi.h"

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;

apps/gdalmdiminfo_lib.cpp

@@ -163,6 +163,9 @@ static void DumpValue(CPLJSonStreamingWriter &serializer, const GByte *bytes,
         case GDT_UInt64:
             DumpValue<std::uint64_t>(serializer, bytes);
             break;
+        case GDT_Float16:
+            DumpValue<GFloat16>(serializer, bytes);
+            break;
         case GDT_Float32:
             DumpValue<float>(serializer, bytes);
             break;
@@ -175,6 +178,9 @@ static void DumpValue(CPLJSonStreamingWriter &serializer, const GByte *bytes,
         case GDT_CInt32:
             DumpComplexValue<GInt32>(serializer, bytes);
             break;
+        case GDT_CFloat16:
+            DumpComplexValue<GFloat16>(serializer, bytes);
+            break;
         case GDT_CFloat32:
             DumpComplexValue<float>(serializer, bytes);
             break;

autotest/cpp/CMakeLists.txt

@@ -290,6 +290,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
@@ -303,6 +304,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)