Skip to content

Commit

Permalink
HDF5: add support for libhdf5 >= 1.14.4.2 when built with Float16
Browse files Browse the repository at this point in the history 
Fixes HDFGroup/hdf5#4527
  • Loading branch information
@rouault
rouault committed May 29, 2024
1 parent 3757cd2 commit 16ade82
Show file tree
Hide file tree
Showing 7 changed files with 293 additions and 16 deletions.
64 changes: 64 additions & 0 deletions autotest/gdrivers/hdf5multidim.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
)
11 changes: 11 additions & 0 deletions frmts/hdf5/gh5_convenience.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
* DEALINGS IN THE SOFTWARE.
****************************************************************************/

#include "cpl_float.h"
#include "gh5_convenience.h"

/************************************************************************/
Expand Down Expand Up @@ -209,6 +210,16 @@ bool GH5_FetchAttribute(hid_t loc_id, const char *pszAttrName, double &dfResult,
pszAttrName, static_cast<GUIntBig>(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))
Expand Down
4 changes: 4 additions & 0 deletions frmts/hdf5/hdf5_api.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
68 changes: 68 additions & 0 deletions frmts/hdf5/hdf5dataset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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))
Expand All @@ -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() */
/* */
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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<uint16_t *>(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++)
Expand Down
4 changes: 2 additions & 2 deletions frmts/hdf5/hdf5dataset.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -275,6 +273,8 @@ class HDF5Dataset CPL_NON_FINAL : public GDALPamDataset
static std::shared_ptr<GDALGroup> OpenGroup(
const std::shared_ptr<GDAL::HDF5SharedResources> &poSharedResources);

static bool IsNativeCFloat16(hid_t hDataType);
static const char *GetDataTypeName(hid_t);
static GDALDataType GetDataType(hid_t);
};

Expand Down
95 changes: 83 additions & 12 deletions frmts/hdf5/hdf5imagedataset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@

#include "hdf5_api.h"

#include "cpl_float.h"
#include "cpl_string.h"
#include "gdal_frmts.h"
#include "gdal_pam.h"
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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);
Expand All @@ -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);

Expand Down Expand Up @@ -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<size_t>(nBlockXSize) * nBlockYSize; i > 0;
/* do nothing */)
{
--i;
uint16_t nVal16;
memcpy(&nVal16, static_cast<uint16_t *>(pImage) + i,
sizeof(nVal16));
const uint32_t nVal32 = CPLHalfToFloat(nVal16);
float fVal;
memcpy(&fVal, &nVal32, sizeof(fVal));
*(static_cast<float *>(pImage) + i) = fVal;
}
}
else if (eDataType == GDT_CFloat32 && poGDS->m_bConvertFromFloat16)
{
for (size_t i = static_cast<size_t>(nBlockXSize) * nBlockYSize; i > 0;
/* do nothing */)
{
--i;
for (int j = 1; j >= 0; --j)
{
uint16_t nVal16;
memcpy(&nVal16, static_cast<uint16_t *>(pImage) + 2 * i + j,
sizeof(nVal16));
const uint32_t nVal32 = CPLHalfToFloat(nVal16);
float fVal;
memcpy(&fVal, &nVal32, sizeof(fVal));
*(static_cast<float *>(pImage) + 2 * i + j) = fVal;
}
}
}
#endif

return CE_None;
}

Expand All @@ -482,6 +518,15 @@ CPLErr HDF5ImageRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
{
HDF5ImageDataset *poGDS = static_cast<HDF5ImageDataset *>(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;
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -1021,9 +1076,25 @@ GDALDataset *HDF5ImageDataset::Open(GDALOpenInfo *poOpenInfo)
static_cast<hsize_t *>(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.
Expand Down Expand Up @@ -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);

Expand Down
Loading

0 comments on commit 16ade82

Please sign in to comment.