diff --git a/autotest/gdrivers/hdf5multidim.py b/autotest/gdrivers/hdf5multidim.py index 6f1f003c326f..31260e9ab766 100755 --- a/autotest/gdrivers/hdf5multidim.py +++ b/autotest/gdrivers/hdf5multidim.py @@ -840,3 +840,67 @@ def test_hdf5_multidim_block_size_structural_info(): var = rg.OpenMDArray("Band1") assert var.GetBlockSize() == [1, 2] assert var.GetStructuralInfo() == {"COMPRESSION": "DEFLATE", "FILTER": "SHUFFLE"} + + +############################################################################### +# Test reading a compound data type made of 2 Float16 values + + +def test_hdf5_multidim_read_cfloat16(): + + ds = gdal.OpenEx("data/hdf5/complex.h5", gdal.OF_MULTIDIM_RASTER) + rg = ds.GetRootGroup() + var = rg.OpenMDArray("f16") + assert var.GetDataType().GetNumericDataType() == gdal.GDT_CFloat32 + assert struct.unpack("f" * (5 * 5 * 2), var.Read()) == ( + 0.0, + 0.0, + 1.0, + 1.0, + 2.0, + 2.0, + 3.0, + 3.0, + 4.0, + 4.0, + 5.0, + 5.0, + 6.0, + 6.0, + 7.0, + 7.0, + 8.0, + 8.0, + 9.0, + 9.0, + 10.0, + 10.0, + 11.0, + 11.0, + 12.0, + 12.0, + 13.0, + 13.0, + 14.0, + 14.0, + 15.0, + 15.0, + 16.0, + 16.0, + 17.0, + 17.0, + 18.0, + 18.0, + 19.0, + 19.0, + 20.0, + 20.0, + 21.0, + 21.0, + 22.0, + 22.0, + 23.0, + 23.0, + 24.0, + 24.0, + ) diff --git a/frmts/hdf5/gh5_convenience.cpp b/frmts/hdf5/gh5_convenience.cpp index 304032e08e3f..e48b5c630d21 100644 --- a/frmts/hdf5/gh5_convenience.cpp +++ b/frmts/hdf5/gh5_convenience.cpp @@ -26,6 +26,7 @@ * DEALINGS IN THE SOFTWARE. ****************************************************************************/ +#include "cpl_float.h" #include "gh5_convenience.h" /************************************************************************/ @@ -209,6 +210,16 @@ bool GH5_FetchAttribute(hid_t loc_id, const char *pszAttrName, double &dfResult, pszAttrName, static_cast(nVal), dfResult); } } +#ifdef HDF5_HAVE_FLOAT16 + else if (H5Tequal(H5T_NATIVE_FLOAT16, hAttrNativeType)) + { + const uint16_t nVal16 = *((uint16_t *)buf); + const uint32_t nVal32 = CPLHalfToFloat(nVal16); + float fVal; + memcpy(&fVal, &nVal32, sizeof(fVal)); + dfResult = fVal; + } +#endif else if (H5Tequal(H5T_NATIVE_FLOAT, hAttrNativeType)) dfResult = *((float *)buf); else if (H5Tequal(H5T_NATIVE_DOUBLE, hAttrNativeType)) diff --git a/frmts/hdf5/hdf5_api.h b/frmts/hdf5/hdf5_api.h index 013534310488..0b3110d2cc2e 100644 --- a/frmts/hdf5/hdf5_api.h +++ b/frmts/hdf5/hdf5_api.h @@ -46,6 +46,10 @@ #include "hdf5.h" +#if defined(H5T_NATIVE_FLOAT16) && defined(H5_HAVE__FLOAT16) +#define HDF5_HAVE_FLOAT16 +#endif + #ifdef _MSC_VER #pragma warning(pop) #endif diff --git a/frmts/hdf5/hdf5dataset.cpp b/frmts/hdf5/hdf5dataset.cpp index 499b8a86c877..64c1fea2e92b 100644 --- a/frmts/hdf5/hdf5dataset.cpp +++ b/frmts/hdf5/hdf5dataset.cpp @@ -44,6 +44,7 @@ #include "cpl_conv.h" #include "cpl_error.h" +#include "cpl_float.h" #include "cpl_string.h" #include "gdal.h" #include "gdal_frmts.h" @@ -201,6 +202,10 @@ GDALDataType HDF5Dataset::GetDataType(hid_t TypeID) return GDT_Unknown; #endif } +#ifdef HDF5_HAVE_FLOAT16 + else if (H5Tequal(H5T_NATIVE_FLOAT16, TypeID)) + return GDT_Float32; +#endif else if (H5Tequal(H5T_NATIVE_FLOAT, TypeID)) return GDT_Float32; else if (H5Tequal(H5T_NATIVE_DOUBLE, TypeID)) @@ -258,6 +263,10 @@ GDALDataType HDF5Dataset::GetDataType(hid_t TypeID) eDataType = GDT_Unknown; #endif } +#ifdef HDF5_HAVE_FLOAT16 + else if (H5Tequal(H5T_NATIVE_FLOAT16, ElemTypeID)) + eDataType = GDT_CFloat32; +#endif else if (H5Tequal(H5T_NATIVE_FLOAT, ElemTypeID)) eDataType = GDT_CFloat32; else if (H5Tequal(H5T_NATIVE_DOUBLE, ElemTypeID)) @@ -272,6 +281,32 @@ GDALDataType HDF5Dataset::GetDataType(hid_t TypeID) return GDT_Unknown; } +/************************************************************************/ +/* IsNativeCFloat16() */ +/************************************************************************/ + +/* static*/ bool HDF5Dataset::IsNativeCFloat16(hid_t hDataType) +{ +#ifdef HDF5_HAVE_FLOAT16 + // For complex the compound type must contain 2 elements + if (H5Tget_class(hDataType) != H5T_COMPOUND || + H5Tget_nmembers(hDataType) != 2) + return false; + + // For complex the native types of both elements should be the same + hid_t ElemTypeID = H5Tget_member_type(hDataType, 0); + hid_t Elem2TypeID = H5Tget_member_type(hDataType, 1); + const bool bRet = H5Tequal(ElemTypeID, H5T_NATIVE_FLOAT16) > 0 && + H5Tequal(Elem2TypeID, H5T_NATIVE_FLOAT16) > 0; + H5Tclose(ElemTypeID); + H5Tclose(Elem2TypeID); + return bRet; +#else + CPL_IGNORE_RET_VAL(hDataType); + return false; +#endif +} + /************************************************************************/ /* GetDataTypeName() */ /* */ @@ -304,6 +339,10 @@ const char *HDF5Dataset::GetDataTypeName(hid_t TypeID) return "32/64-bit integer"; else if (H5Tequal(H5T_NATIVE_ULONG, TypeID)) return "32/64-bit unsigned integer"; +#ifdef HDF5_HAVE_FLOAT16 + else if (H5Tequal(H5T_NATIVE_FLOAT16, TypeID)) + return "16-bit floating-point"; +#endif else if (H5Tequal(H5T_NATIVE_FLOAT, TypeID)) return "32-bit floating-point"; else if (H5Tequal(H5T_NATIVE_DOUBLE, TypeID)) @@ -348,6 +387,13 @@ const char *HDF5Dataset::GetDataTypeName(hid_t TypeID) H5Tclose(ElemTypeID); return "complex, 32/64-bit integer"; } +#ifdef HDF5_HAVE_FLOAT16 + else if (H5Tequal(H5T_NATIVE_FLOAT16, ElemTypeID)) + { + H5Tclose(ElemTypeID); + return "complex, 16-bit floating-point"; + } +#endif else if (H5Tequal(H5T_NATIVE_FLOAT, ElemTypeID)) { H5Tclose(ElemTypeID); @@ -1132,6 +1178,28 @@ static herr_t HDF5AttrIterate(hid_t hH5ObjID, const char *pszAttrName, psContext->m_osValue += szData; } } +#ifdef HDF5_HAVE_FLOAT16 + else if (H5Tequal(H5T_NATIVE_FLOAT16, hAttrNativeType) > 0) + { + for (hsize_t i = 0; i < nAttrElmts; i++) + { + const uint16_t nVal16 = static_cast(buf)[i]; + const uint32_t nVal32 = CPLHalfToFloat(nVal16); + float fVal; + memcpy(&fVal, &nVal32, sizeof(fVal)); + CPLsnprintf(szData, nDataLen, "%.8g", fVal); + if (psContext->m_osValue.size() > MAX_METADATA_LEN) + { + CPLError(CE_Warning, CPLE_OutOfMemory, + "Header data too long. Truncated"); + break; + } + if (i > 0) + psContext->m_osValue += ' '; + psContext->m_osValue += szData; + } + } +#endif else if (H5Tequal(H5T_NATIVE_FLOAT, hAttrNativeType) > 0) { for (hsize_t i = 0; i < nAttrElmts; i++) diff --git a/frmts/hdf5/hdf5dataset.h b/frmts/hdf5/hdf5dataset.h index dc9245b4ca31..646f7c6b6297 100644 --- a/frmts/hdf5/hdf5dataset.h +++ b/frmts/hdf5/hdf5dataset.h @@ -240,8 +240,6 @@ class HDF5Dataset CPL_NON_FINAL : public GDALPamDataset char *CreatePath(HDF5GroupObjects *); static void DestroyH5Objects(HDF5GroupObjects *); - static const char *GetDataTypeName(hid_t); - /** * Reads an array of double attributes from the HDF5 metadata. * It reads the attributes directly on its binary form directly, @@ -275,6 +273,8 @@ class HDF5Dataset CPL_NON_FINAL : public GDALPamDataset static std::shared_ptr OpenGroup( const std::shared_ptr &poSharedResources); + static bool IsNativeCFloat16(hid_t hDataType); + static const char *GetDataTypeName(hid_t); static GDALDataType GetDataType(hid_t); }; diff --git a/frmts/hdf5/hdf5imagedataset.cpp b/frmts/hdf5/hdf5imagedataset.cpp index 64a2c33917ea..4d140c523f28 100644 --- a/frmts/hdf5/hdf5imagedataset.cpp +++ b/frmts/hdf5/hdf5imagedataset.cpp @@ -29,6 +29,7 @@ #include "hdf5_api.h" +#include "cpl_float.h" #include "cpl_string.h" #include "gdal_frmts.h" #include "gdal_pam.h" @@ -72,9 +73,10 @@ class HDF5ImageDataset final : public HDF5Dataset int dimensions; hid_t dataset_id; hid_t dataspace_id; - hsize_t size; - hid_t datatype; hid_t native; +#ifdef HDF5_HAVE_FLOAT16 + bool m_bConvertFromFloat16 = false; +#endif Hdf5ProductType iSubdatasetType; HDF5CSKProductEnum iCSKProductType; double adfGeoTransform[6]; @@ -209,9 +211,9 @@ class HDF5ImageDataset final : public HDF5Dataset /************************************************************************/ HDF5ImageDataset::HDF5ImageDataset() : dims(nullptr), maxdims(nullptr), poH5Objects(nullptr), ndims(0), - dimensions(0), dataset_id(-1), dataspace_id(-1), size(0), datatype(-1), - native(-1), iSubdatasetType(UNKNOWN_PRODUCT), - iCSKProductType(PROD_UNKNOWN), bHasGeoTransform(false) + dimensions(0), dataset_id(-1), dataspace_id(-1), native(-1), + iSubdatasetType(UNKNOWN_PRODUCT), iCSKProductType(PROD_UNKNOWN), + bHasGeoTransform(false) { m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); @@ -236,8 +238,6 @@ HDF5ImageDataset::~HDF5ImageDataset() H5Dclose(dataset_id); if (dataspace_id > 0) H5Sclose(dataspace_id); - if (datatype > 0) - H5Tclose(datatype); if (native > 0) H5Tclose(native); @@ -465,6 +465,42 @@ CPLErr HDF5ImageRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, return CE_Failure; } +#ifdef HDF5_HAVE_FLOAT16 + if (eDataType == GDT_Float32 && poGDS->m_bConvertFromFloat16) + { + for (size_t i = static_cast(nBlockXSize) * nBlockYSize; i > 0; + /* do nothing */) + { + --i; + uint16_t nVal16; + memcpy(&nVal16, static_cast(pImage) + i, + sizeof(nVal16)); + const uint32_t nVal32 = CPLHalfToFloat(nVal16); + float fVal; + memcpy(&fVal, &nVal32, sizeof(fVal)); + *(static_cast(pImage) + i) = fVal; + } + } + else if (eDataType == GDT_CFloat32 && poGDS->m_bConvertFromFloat16) + { + for (size_t i = static_cast(nBlockXSize) * nBlockYSize; i > 0; + /* do nothing */) + { + --i; + for (int j = 1; j >= 0; --j) + { + uint16_t nVal16; + memcpy(&nVal16, static_cast(pImage) + 2 * i + j, + sizeof(nVal16)); + const uint32_t nVal32 = CPLHalfToFloat(nVal16); + float fVal; + memcpy(&fVal, &nVal32, sizeof(fVal)); + *(static_cast(pImage) + 2 * i + j) = fVal; + } + } + } +#endif + return CE_None; } @@ -482,6 +518,15 @@ CPLErr HDF5ImageRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, { HDF5ImageDataset *poGDS = static_cast(poDS); +#ifdef HDF5_HAVE_FLOAT16 + if (poGDS->m_bConvertFromFloat16) + { + return GDALPamRasterBand::IRasterIO( + eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize, + eBufType, nPixelSpace, nLineSpace, psExtraArg); + } +#endif + const bool bIsBandInterleavedData = poGDS->ndims == 3 && poGDS->m_nOtherDimIndex == 0 && poGDS->GetYIndex() == 1 && poGDS->GetXIndex() == 2; @@ -761,6 +806,16 @@ CPLErr HDF5ImageDataset::IRasterIO( GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg) { +#ifdef HDF5_HAVE_FLOAT16 + if (m_bConvertFromFloat16) + { + return HDF5Dataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, + pData, nBufXSize, nBufYSize, eBufType, + nBandCount, panBandMap, nPixelSpace, + nLineSpace, nBandSpace, psExtraArg); + } +#endif + const auto IsConsecutiveBands = [](const int *panVals, int nCount) { for (int i = 1; i < nCount; ++i) @@ -1021,9 +1076,25 @@ GDALDataset *HDF5ImageDataset::Open(GDALOpenInfo *poOpenInfo) static_cast(CPLCalloc(poDS->ndims, sizeof(hsize_t))); poDS->dimensions = H5Sget_simple_extent_dims(poDS->dataspace_id, poDS->dims, poDS->maxdims); - poDS->datatype = H5Dget_type(poDS->dataset_id); - poDS->size = H5Tget_size(poDS->datatype); - poDS->native = H5Tget_native_type(poDS->datatype, H5T_DIR_ASCEND); + auto datatype = H5Dget_type(poDS->dataset_id); + poDS->native = H5Tget_native_type(datatype, H5T_DIR_ASCEND); + H5Tclose(datatype); + + const auto eGDALDataType = poDS->GetDataType(poDS->native); + if (eGDALDataType == GDT_Unknown) + { + CPLError(CE_Failure, CPLE_AppDefined, "Unhandled HDF5 data type"); + delete poDS; + return nullptr; + } + +#ifdef HDF5_HAVE_FLOAT16 + if (H5Tequal(H5T_NATIVE_FLOAT16, poDS->native) || + IsNativeCFloat16(poDS->native)) + { + poDS->m_bConvertFromFloat16 = true; + } +#endif // CSK code in IdentifyProductType() and CreateProjections() // uses dataset metadata. @@ -1239,8 +1310,8 @@ GDALDataset *HDF5ImageDataset::Open(GDALOpenInfo *poOpenInfo) for (int i = 0; i < nBands; i++) { - HDF5ImageRasterBand *const poBand = new HDF5ImageRasterBand( - poDS, i + 1, poDS->GetDataType(poDS->native)); + HDF5ImageRasterBand *const poBand = + new HDF5ImageRasterBand(poDS, i + 1, eGDALDataType); poDS->SetBand(i + 1, poBand); diff --git a/frmts/hdf5/hdf5multidim.cpp b/frmts/hdf5/hdf5multidim.cpp index fa27ff79a051..f5a0704c458f 100644 --- a/frmts/hdf5/hdf5multidim.cpp +++ b/frmts/hdf5/hdf5multidim.cpp @@ -29,6 +29,8 @@ #include "hdf5eosparser.h" #include "s100.h" +#include "cpl_float.h" + #include #include #include @@ -159,7 +161,16 @@ BuildDataType(hid_t hDataType, bool &bHasString, bool &bNonNativeDataType, const auto klass = H5Tget_class(hDataType); GDALDataType eDT = ::HDF5Dataset::GetDataType(hDataType); if (eDT != GDT_Unknown) + { +#ifdef HDF5_HAVE_FLOAT16 + if (H5Tequal(hDataType, H5T_NATIVE_FLOAT16) || + HDF5Dataset::IsNativeCFloat16(hDataType)) + { + bNonNativeDataType = true; + } +#endif return GDALExtendedDataType::Create(eDT); + } else if (klass == H5T_STRING) { bHasString = true; @@ -2326,9 +2337,44 @@ static void CopyValue(const GByte *pabySrcBuffer, hid_t hSrcDataType, { if (dstDataType.GetClass() != GEDTC_COMPOUND) { + const auto eSrcDataType = ::HDF5Dataset::GetDataType(hSrcDataType); // Typically source is complex data type - auto srcDataType(GDALExtendedDataType::Create( - ::HDF5Dataset::GetDataType(hSrcDataType))); +#ifdef HDF5_HAVE_FLOAT16 + if (eSrcDataType == GDT_CFloat32 && + ::HDF5Dataset::IsNativeCFloat16(hSrcDataType)) + { + if (dstDataType.GetNumericDataType() == GDT_CFloat32) + { + for (int j = 0; j <= 1; ++j) + { + uint16_t nVal16; + memcpy(&nVal16, pabySrcBuffer + j * sizeof(nVal16), + sizeof(nVal16)); + const uint32_t nVal32 = CPLHalfToFloat(nVal16); + memcpy(pabyDstBuffer + j * sizeof(float), &nVal32, + sizeof(nVal32)); + } + } + else if (dstDataType.GetNumericDataType() == GDT_CFloat64) + { + for (int j = 0; j <= 1; ++j) + { + uint16_t nVal16; + memcpy(&nVal16, pabySrcBuffer + j * sizeof(nVal16), + sizeof(nVal16)); + const uint32_t nVal32 = CPLHalfToFloat(nVal16); + float fVal; + memcpy(&fVal, &nVal32, sizeof(fVal)); + double dfVal = fVal; + memcpy(pabyDstBuffer + j * sizeof(double), &dfVal, + sizeof(dfVal)); + } + } + return; + } + +#endif + auto srcDataType(GDALExtendedDataType::Create(eSrcDataType)); if (srcDataType.GetClass() == GEDTC_NUMERIC && srcDataType.GetNumericDataType() != GDT_Unknown) { @@ -2364,6 +2410,19 @@ static void CopyValue(const GByte *pabySrcBuffer, hid_t hSrcDataType, CopyValue(pabySrcBuffer, hParent, pabyDstBuffer, dstDataType, {}); H5Tclose(hParent); } +#ifdef HDF5_HAVE_FLOAT16 + else if (H5Tequal(hSrcDataType, H5T_NATIVE_FLOAT16)) + { + uint16_t nVal16; + memcpy(&nVal16, pabySrcBuffer, sizeof(nVal16)); + const uint32_t nVal32 = CPLHalfToFloat(nVal16); + float fVal; + memcpy(&fVal, &nVal32, sizeof(fVal)); + GDALExtendedDataType::CopyValue( + &fVal, GDALExtendedDataType::Create(GDT_Float32), pabyDstBuffer, + dstDataType); + } +#endif else { GDALDataType eDT = ::HDF5Dataset::GetDataType(hSrcDataType);