/******************************************************************************
 *
 * Project:  Hierarchical Data Format Release 4 (HDF4)
 * Purpose:  HDF4 Datasets. Open HDF4 file, fetch metadata and list of
 *           subdatasets.
 *           This driver initially based on code supplied by Markus Neteler
 * Author:   Andrey Kiselev, dron@ak4719.spb.edu
 *
 ******************************************************************************
 * Copyright (c) 2002, Andrey Kiselev <dron@ak4719.spb.edu>
 * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_multiproc.h"
#include "cpl_string.h"
#include "gdal_frmts.h"
#include "gdal_priv.h"

#include "hdf.h"
#include "mfhdf.h"

#include "HdfEosDef.h"

#include "hdf4compat.h"
#include "hdf4dataset.h"
#include <cctype>

#include "hdf4drivercore.h"

extern const char *const pszGDALSignature;

CPLMutex *hHDF4Mutex = nullptr;

/************************************************************************/
/* ==================================================================== */
/*                              HDF4Dataset                             */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                           HDF4Dataset()                              */
/************************************************************************/

HDF4Dataset::HDF4Dataset()
    : bIsHDFEOS(false), hGR(0), hSD(0), nImages(0),
      iSubdatasetType(H4ST_UNKNOWN), pszSubdatasetType(nullptr),
      papszGlobalMetadata(nullptr), papszSubDatasets(nullptr)
{
}

/************************************************************************/
/*                            ~HDF4Dataset()                            */
/************************************************************************/

HDF4Dataset::~HDF4Dataset()

{
    CPLMutexHolderD(&hHDF4Mutex);

    if (hSD)
        SDend(hSD);
    if (hGR)
        GRend(hGR);
    if (papszSubDatasets)
        CSLDestroy(papszSubDatasets);
    if (papszGlobalMetadata)
        CSLDestroy(papszGlobalMetadata);
}

/************************************************************************/
/*                      GetMetadataDomainList()                         */
/************************************************************************/

char **HDF4Dataset::GetMetadataDomainList()
{
    return BuildMetadataDomainList(GDALPamDataset::GetMetadataDomainList(),
                                   TRUE, "SUBDATASETS", nullptr);
}

/************************************************************************/
/*                            GetMetadata()                             */
/************************************************************************/

char **HDF4Dataset::GetMetadata(const char *pszDomain)

{
    if (pszDomain != nullptr && STARTS_WITH_CI(pszDomain, "SUBDATASETS"))
        return papszSubDatasets;

    return GDALDataset::GetMetadata(pszDomain);
}

/************************************************************************/
/*                           SPrintArray()                              */
/*      Prints numerical arrays in string buffer.                       */
/*      This function takes pfaDataArray as a pointer to printed array, */
/*      nValues as a number of values to print and pszDelimiter as a    */
/*      field delimiting strings.                                       */
/*      Pointer to filled buffer will be returned.                      */
/************************************************************************/

char *SPrintArray(GDALDataType eDataType, const void *paDataArray, int nValues,
                  const char *pszDelimiter)
{
    const int iFieldSize = 32 + static_cast<int>(strlen(pszDelimiter));
    char *pszField = static_cast<char *>(CPLMalloc(iFieldSize + 1));
    const int iStringSize = nValues * iFieldSize + 1;
    char *pszString = static_cast<char *>(CPLMalloc(iStringSize));
    memset(pszString, 0, iStringSize);
    for (int i = 0; i < nValues; i++)
    {
        switch (eDataType)
        {
            case GDT_Byte:
                snprintf(pszField, iFieldSize + 1, "%d%s",
                         reinterpret_cast<GByte *>(
                             const_cast<void *>(paDataArray))[i],
                         (i < nValues - 1) ? pszDelimiter : "");
                break;
            case GDT_Int8:
                snprintf(pszField, iFieldSize + 1, "%d%s",
                         reinterpret_cast<GInt8 *>(
                             const_cast<void *>(paDataArray))[i],
                         (i < nValues - 1) ? pszDelimiter : "");
                break;
            case GDT_UInt16:
                snprintf(pszField, iFieldSize + 1, "%u%s",
                         reinterpret_cast<GUInt16 *>(
                             const_cast<void *>(paDataArray))[i],
                         (i < nValues - 1) ? pszDelimiter : "");
                break;
            case GDT_Int16:
            default:
                snprintf(pszField, iFieldSize + 1, "%d%s",
                         reinterpret_cast<GInt16 *>(
                             const_cast<void *>(paDataArray))[i],
                         (i < nValues - 1) ? pszDelimiter : "");
                break;
            case GDT_UInt32:
                snprintf(pszField, iFieldSize + 1, "%u%s",
                         reinterpret_cast<GUInt32 *>(
                             const_cast<void *>(paDataArray))[i],
                         (i < nValues - 1) ? pszDelimiter : "");
                break;
            case GDT_Int32:
                snprintf(pszField, iFieldSize + 1, "%d%s",
                         reinterpret_cast<GInt32 *>(
                             const_cast<void *>(paDataArray))[i],
                         (i < nValues - 1) ? pszDelimiter : "");
                break;
            case GDT_Float32:
                CPLsnprintf(pszField, iFieldSize + 1, "%.10g%s",
                            reinterpret_cast<float *>(
                                const_cast<void *>(paDataArray))[i],
                            (i < nValues - 1) ? pszDelimiter : "");
                break;
            case GDT_Float64:
                CPLsnprintf(pszField, iFieldSize + 1, "%.15g%s",
                            reinterpret_cast<double *>(
                                const_cast<void *>(paDataArray))[i],
                            (i < nValues - 1) ? pszDelimiter : "");
                break;
        }
        strcat(pszString, pszField);
    }

    CPLFree(pszField);
    return pszString;
}

/************************************************************************/
/*              Translate HDF4 data type into GDAL data type            */
/************************************************************************/
GDALDataType HDF4Dataset::GetDataType(int32 iNumType)
{
    switch (iNumType)
    {
        case DFNT_CHAR8:   // The same as DFNT_CHAR
        case DFNT_UCHAR8:  // The same as DFNT_UCHAR
        case DFNT_UINT8:
            return GDT_Byte;
        case DFNT_INT8:
            return GDT_Int8;
        case DFNT_INT16:
            return GDT_Int16;
        case DFNT_UINT16:
            return GDT_UInt16;
        case DFNT_INT32:
            return GDT_Int32;
        case DFNT_UINT32:
            return GDT_UInt32;
        case DFNT_INT64:
            return GDT_Int64;
        case DFNT_UINT64:
            return GDT_UInt64;
        case DFNT_FLOAT32:
            return GDT_Float32;
        case DFNT_FLOAT64:
            return GDT_Float64;
        default:
            return GDT_Unknown;
    }
}

/************************************************************************/
/*              Return the human readable name of data type             */
/************************************************************************/

const char *HDF4Dataset::GetDataTypeName(int32 iNumType)
{
    switch (iNumType)
    {
        case DFNT_CHAR8:  // The same as DFNT_CHAR
            return "8-bit character";
        case DFNT_UCHAR8:  // The same as DFNT_UCHAR
            return "8-bit unsigned character";
        case DFNT_INT8:
            return "8-bit integer";
        case DFNT_UINT8:
            return "8-bit unsigned integer";
        case DFNT_INT16:
            return "16-bit integer";
        case DFNT_UINT16:
            return "16-bit unsigned integer";
        case DFNT_INT32:
            return "32-bit integer";
        case DFNT_UINT32:
            return "32-bit unsigned integer";
        case DFNT_INT64:
            return "64-bit integer";
        case DFNT_UINT64:
            return "64-bit unsigned integer";
        case DFNT_FLOAT32:
            return "32-bit floating-point";
        case DFNT_FLOAT64:
            return "64-bit floating-point";
        default:
        {
            CPLError(CE_Warning, CPLE_AppDefined, "Unknown type %d",
                     static_cast<int>(iNumType));

            return "unknown type";
        }
    }
}

/************************************************************************/
/*  Return the size of data type in bytes                               */
/************************************************************************/

int HDF4Dataset::GetDataTypeSize(int32 iNumType)
{
    switch (iNumType)
    {
        case DFNT_CHAR8:   // The same as DFNT_CHAR
        case DFNT_UCHAR8:  // The same as DFNT_UCHAR
        case DFNT_INT8:
        case DFNT_UINT8:
            return 1;
        case DFNT_INT16:
        case DFNT_UINT16:
            return 2;
        case DFNT_INT32:
        case DFNT_UINT32:
        case DFNT_FLOAT32:
            return 4;
        case DFNT_INT64:
        case DFNT_UINT64:
        case DFNT_FLOAT64:
            return 8;
        default:
        {
            CPLError(CE_Warning, CPLE_AppDefined, "Unknown type %d",
                     static_cast<int>(iNumType));
            return 0;
        }
    }
}

/************************************************************************/
/*  Convert value stored in the input buffer to double value.           */
/************************************************************************/

double HDF4Dataset::AnyTypeToDouble(int32 iNumType, void *pData)
{
    switch (iNumType)
    {
        case DFNT_INT8:
            return static_cast<double>(*reinterpret_cast<signed char *>(pData));
        case DFNT_UINT8:
            return static_cast<double>(*reinterpret_cast<GByte *>(pData));
        case DFNT_INT16:
            return static_cast<double>(*reinterpret_cast<GInt16 *>(pData));
        case DFNT_UINT16:
            return static_cast<double>(*reinterpret_cast<GUInt16 *>(pData));
        case DFNT_INT32:
            return static_cast<double>(*reinterpret_cast<GInt32 *>(pData));
        case DFNT_UINT32:
            return static_cast<double>(*reinterpret_cast<GUInt32 *>(pData));
        case DFNT_INT64:
            return static_cast<double>(*reinterpret_cast<GInt64 *>(pData));
        case DFNT_UINT64:
            return static_cast<double>(*reinterpret_cast<GUInt64 *>(pData));
        case DFNT_FLOAT32:
            return static_cast<double>(*reinterpret_cast<float *>(pData));
        case DFNT_FLOAT64:
            return static_cast<double>(*reinterpret_cast<double *>(pData));
        default:
        {
            CPLError(CE_Warning, CPLE_AppDefined, "Unknown type %d",
                     static_cast<int>(iNumType));
            return 0.0;
        }
    }
}

/************************************************************************/
/*         Tokenize HDF-EOS attributes.                                 */
/************************************************************************/

char **HDF4Dataset::HDF4EOSTokenizeAttrs(const char *pszString)

{
    const char *const pszDelimiters = " \t\n\r";
    char **papszRetList = nullptr;

    char *pszToken = static_cast<char *>(CPLCalloc(10, 1));
    int nTokenMax = 10;

    while (pszString != nullptr && *pszString != '\0')
    {
        bool bInString = false;
        bool bInBracket = false;

        int nTokenLen = 0;

        // Try to find the next delimiter, marking end of token.
        for (; *pszString != '\0'; pszString++)
        {

            // End if this is a delimiter skip it and break.
            if (!bInBracket && !bInString &&
                strchr(pszDelimiters, *pszString) != nullptr)
            {
                pszString++;
                break;
            }

            // Sometimes in bracketed tokens we may found a sort of
            // paragraph formatting. We will remove unneeded spaces and new
            // lines.
            if (bInBracket)
                if (strchr("\r\n", *pszString) != nullptr ||
                    (*pszString == ' ' &&
                     strchr(" \r\n", *(pszString - 1)) != nullptr))
                    continue;

            if (*pszString == '"')
            {
                bInString = !bInString;
                continue;
            }
            else if (*pszString == '(')
            {
                bInBracket = true;
                continue;
            }
            else if (*pszString == ')')
            {
                bInBracket = false;
                continue;
            }

            if (nTokenLen >= nTokenMax - 2)
            {
                nTokenMax = nTokenMax * 2 + 10;
                pszToken = static_cast<char *>(CPLRealloc(pszToken, nTokenMax));
            }

            pszToken[nTokenLen] = *pszString;
            nTokenLen++;
        }

        pszToken[nTokenLen] = '\0';

        if (pszToken[0] != '\0')
        {
            papszRetList = CSLAddString(papszRetList, pszToken);
        }

        // If the last token is an empty token, then we have to catch
        // it now, otherwise we won't reenter the loop and it will be lost.
        if (*pszString == '\0' && strchr(pszDelimiters, *(pszString - 1)))
        {
            papszRetList = CSLAddString(papszRetList, "");
        }
    }

    if (papszRetList == nullptr)
        papszRetList = static_cast<char **>(CPLCalloc(sizeof(char *), 1));

    CPLFree(pszToken);

    return papszRetList;
}

/************************************************************************/
/*     Find object name, class value in HDF-EOS attributes.             */
/*     Function returns pointer to the string in list next behind       */
/*     recognized object.                                               */
/************************************************************************/

char **HDF4Dataset::HDF4EOSGetObject(char **papszAttrList, char **ppszAttrName,
                                     char **ppszAttrClass, char **ppszAttrValue)
{
    *ppszAttrName = nullptr;
    *ppszAttrClass = nullptr;
    *ppszAttrValue = nullptr;

    const int iCount = CSLCount(papszAttrList);
    for (int i = 0; i < iCount - 2; i++)
    {
        if (EQUAL(papszAttrList[i], "OBJECT"))
        {
            i += 2;
            for (int j = 1; i + j < iCount - 2; j++)
            {
                if (EQUAL(papszAttrList[i + j], "END_OBJECT") ||
                    EQUAL(papszAttrList[i + j], "OBJECT"))
                    return &papszAttrList[i + j];
                else if (EQUAL(papszAttrList[i + j], "CLASS"))
                {
                    *ppszAttrClass = papszAttrList[i + j + 2];
                    continue;
                }
                else if (EQUAL(papszAttrList[i + j], "VALUE"))
                {
                    *ppszAttrName = papszAttrList[i];
                    *ppszAttrValue = papszAttrList[i + j + 2];
                    continue;
                }
            }
        }
    }

    return nullptr;
}

/************************************************************************/
/*         Translate HDF4-EOS attributes in GDAL metadata items         */
/************************************************************************/

char **HDF4Dataset::TranslateHDF4EOSAttributes(int32 iHandle, int32 iAttribute,
                                               int32 nValues,
                                               char **papszMetadata)
{
    char *pszData =
        static_cast<char *>(CPLMalloc((nValues + 1) * sizeof(char)));
    pszData[nValues] = '\0';
    SDreadattr(iHandle, iAttribute, pszData);
    // HDF4-EOS attributes has followed structure:
    //
    // GROUP = <name>
    //   GROUPTYPE = <name>
    //
    //   GROUP = <name>
    //
    //     OBJECT = <name>
    //       CLASS = <string>
    //       NUM_VAL = <number>
    //       VALUE = <string> or <number>
    //     END_OBJECT = <name>
    //
    //     .......
    //     .......
    //     .......
    //
    //   END_GROUP = <name>
    //
    // .......
    // .......
    // .......
    //
    // END_GROUP = <name>
    // END
    //
    // Used data types:
    // <name>   --- unquoted character strings
    // <string> --- quoted character strings
    // <number> --- numerical value
    // If NUM_VAL != 1 then values combined in lists:
    // (<string>,<string>,...)
    // or
    // (<number>,<number>,...)
    //
    // Records within objects could come in any order, objects could contain
    // other objects (and lack VALUE record), groups could contain other groups
    // and objects. Names of groups and objects are not unique and may repeat.
    // Objects may contains other types of records.
    //
    // We are interested in OBJECTS structures only. To avoid multiple items
    // with the same name, names will be suffixed with the class values, e.g.
    //
    //  OBJECT                 = PARAMETERNAME
    //    CLASS                = "9"
    //    NUM_VAL              = 1
    //    VALUE                = "Spectral IR Surf Bidirect Reflectivity"
    //  END_OBJECT             = PARAMETERNAME
    //
    //  will be translated into metadata record:
    //
    //  PARAMETERNAME.9 = "Spectral IR Surf Bidirect Reflectivity"

    char *pszAttrName = nullptr;
    char *pszAttrClass = nullptr;
    char *pszAttrValue = nullptr;
    char *pszAddAttrName = nullptr;

    char **const papszAttrList = HDF4EOSTokenizeAttrs(pszData);
    char **papszAttrs = papszAttrList;
    while (papszAttrs)
    {
        papszAttrs = HDF4EOSGetObject(papszAttrs, &pszAttrName, &pszAttrClass,
                                      &pszAttrValue);
        if (pszAttrName && pszAttrValue)
        {
            // Now we should recognize special type of HDF EOS metastructures:
            // ADDITIONALATTRIBUTENAME = <name>
            // PARAMETERVALUE = <value>
            if (EQUAL(pszAttrName, "ADDITIONALATTRIBUTENAME"))
            {
                pszAddAttrName = pszAttrValue;
            }
            else if (pszAddAttrName && EQUAL(pszAttrName, "PARAMETERVALUE"))
            {
                papszMetadata = CSLAddNameValue(papszMetadata, pszAddAttrName,
                                                pszAttrValue);
                pszAddAttrName = nullptr;
            }
            else
            {
                // Add class suffix to the key name if applicable.
                papszMetadata = CSLAddNameValue(
                    papszMetadata,
                    pszAttrClass
                        ? CPLSPrintf("%s.%s", pszAttrName, pszAttrClass)
                        : pszAttrName,
                    pszAttrValue);
            }
        }
    }

    CSLDestroy(papszAttrList);
    CPLFree(pszData);

    return papszMetadata;
}

/************************************************************************/
/*         Translate HDF4 attributes in GDAL metadata items             */
/************************************************************************/

char **HDF4Dataset::TranslateHDF4Attributes(int32 iHandle, int32 iAttribute,
                                            char *pszAttrName, int32 iNumType,
                                            int32 nValues, char **papszMetadata)
{

    /* -------------------------------------------------------------------- */
    /*     Allocate a buffer to hold the attribute data.                    */
    /* -------------------------------------------------------------------- */
    void *pData = nullptr;
    if (iNumType == DFNT_CHAR8 || iNumType == DFNT_UCHAR8)
        pData = CPLMalloc((nValues + 1) * GetDataTypeSize(iNumType));
    else
        pData = CPLMalloc(nValues * GetDataTypeSize(iNumType));

    /* -------------------------------------------------------------------- */
    /*     Read the attribute data.                                         */
    /* -------------------------------------------------------------------- */
    SDreadattr(iHandle, iAttribute, pData);
    if (iNumType == DFNT_CHAR8 || iNumType == DFNT_UCHAR8)
    {
        reinterpret_cast<char *>(pData)[nValues] = '\0';
        papszMetadata = CSLAddNameValue(
            papszMetadata, pszAttrName,
            const_cast<const char *>(reinterpret_cast<char *>(pData)));
    }
    else
    {
        char *pszTemp =
            SPrintArray(GetDataType(iNumType), pData, nValues, ", ");
        papszMetadata = CSLAddNameValue(papszMetadata, pszAttrName, pszTemp);
        CPLFree(pszTemp);
    }

    CPLFree(pData);

    return papszMetadata;
}

/************************************************************************/
/*                       ReadGlobalAttributes()                         */
/************************************************************************/

CPLErr HDF4Dataset::ReadGlobalAttributes(int32 iHandler)
{
    /* -------------------------------------------------------------------- */
    /*     Obtain number of SDSs and global attributes in input file.       */
    /* -------------------------------------------------------------------- */
    int32 nDatasets = 0;
    int32 nAttributes = 0;
    if (SDfileinfo(iHandler, &nDatasets, &nAttributes) != 0)
        return CE_Failure;

    char szAttrName[H4_MAX_NC_NAME] = {};  // TODO: Get this off the stack.

    // Loop through the all attributes
    for (int32 iAttribute = 0; iAttribute < nAttributes; iAttribute++)
    {
        int32 iNumType = 0;
        int32 nValues = 0;

        // Get information about the attribute. Note that the first
        // parameter is an SD interface identifier.
        SDattrinfo(iHandler, iAttribute, szAttrName, &iNumType, &nValues);

        if (STARTS_WITH_CI(szAttrName, "coremetadata") ||
            STARTS_WITH_CI(szAttrName, "archivemetadata.") ||
            STARTS_WITH_CI(szAttrName, "productmetadata.") ||
            STARTS_WITH_CI(szAttrName, "badpixelinformation") ||
            STARTS_WITH_CI(szAttrName, "product_summary") ||
            STARTS_WITH_CI(szAttrName, "dem_specific") ||
            STARTS_WITH_CI(szAttrName, "bts_specific") ||
            STARTS_WITH_CI(szAttrName, "etse_specific") ||
            STARTS_WITH_CI(szAttrName, "dst_specific") ||
            STARTS_WITH_CI(szAttrName, "acv_specific") ||
            STARTS_WITH_CI(szAttrName, "act_specific") ||
            STARTS_WITH_CI(szAttrName, "etst_specific") ||
            STARTS_WITH_CI(szAttrName, "level_1_carryover"))
        {
            bIsHDFEOS = true;
            papszGlobalMetadata = TranslateHDF4EOSAttributes(
                iHandler, iAttribute, nValues, papszGlobalMetadata);
        }

        // Skip "StructMetadata.N" records. We will fetch information
        // from them using HDF-EOS API
        else if (STARTS_WITH_CI(szAttrName, "structmetadata."))
        {
            bIsHDFEOS = true;
            continue;
        }

        else
        {
            papszGlobalMetadata =
                TranslateHDF4Attributes(iHandler, iAttribute, szAttrName,
                                        iNumType, nValues, papszGlobalMetadata);
        }
    }

    return CE_None;
}

/************************************************************************/
/*                            QuoteIfNeeded()                           */
/************************************************************************/

static CPLString QuoteIfNeeded(const CPLString &osStr)
{
    if (osStr.find(' ') != std::string::npos ||
        osStr.find(':') != std::string::npos ||
        osStr.find('"') != std::string::npos ||
        osStr.find('\\') != std::string::npos)
    {
        CPLString osRet;
        for (size_t i = 0; i < osStr.size(); i++)
        {
            if (osStr[i] == '"')
                osRet += "\\\"";
            else if (osStr[i] == '\\')
                osRet += "\\\\";
            else
                osRet += osStr[i];
        }
        return '"' + osRet + '"';
    }
    return osStr;
}

/************************************************************************/
/*                                Open()                                */
/************************************************************************/

GDALDataset *HDF4Dataset::Open(GDALOpenInfo *poOpenInfo)

{
#ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
    // During fuzzing, do not use Identify to reject crazy content.
    if (!HDF4DatasetIdentify(poOpenInfo))
        return nullptr;
#endif

    CPLMutexHolderD(&hHDF4Mutex);

    /* -------------------------------------------------------------------- */
    /*      Try opening the dataset.                                        */
    /* -------------------------------------------------------------------- */

    // Attempt to increase maximum number of opened HDF files.
#ifdef HDF4_HAS_MAXOPENFILES
    intn nCurrMax = 0;
    intn nSysLimit = 0;

    if (SDget_maxopenfiles(&nCurrMax, &nSysLimit) >= 0 && nCurrMax < nSysLimit)
    {
        /*intn res = */ SDreset_maxopenfiles(nSysLimit);
    }
#endif /* HDF4_HAS_MAXOPENFILES */

    int32 hHDF4 = Hopen(poOpenInfo->pszFilename, DFACC_READ, 0);

    if (hHDF4 <= 0)
        return nullptr;

    Hclose(hHDF4);

    /* -------------------------------------------------------------------- */
    /*      Create a corresponding GDALDataset.                             */
    /* -------------------------------------------------------------------- */
    // Release mutex otherwise we will deadlock with GDALDataset own mutex.
    CPLReleaseMutex(hHDF4Mutex);
    HDF4Dataset *poDS = new HDF4Dataset();
    CPLAcquireMutex(hHDF4Mutex, 1000.0);

    if (poOpenInfo->fpL != nullptr)
    {
        VSIFCloseL(poOpenInfo->fpL);
        poOpenInfo->fpL = nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*          Open HDF SDS Interface.                                     */
    /* -------------------------------------------------------------------- */
    poDS->hSD = SDstart(poOpenInfo->pszFilename, DFACC_READ);

    if (poDS->hSD == -1)
    {
        // Release mutex otherwise we will deadlock with GDALDataset own mutex.
        CPLReleaseMutex(hHDF4Mutex);
        delete poDS;
        CPLAcquireMutex(hHDF4Mutex, 1000.0);
        CPLError(CE_Failure, CPLE_OpenFailed,
                 "Failed to open HDF4 file \"%s\" for SDS reading.",
                 poOpenInfo->pszFilename);
        return nullptr;
    }

    if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
    {
        poDS->OpenMultiDim(poOpenInfo->pszFilename,
                           poOpenInfo->papszOpenOptions);
        return poDS;
    }

    /* -------------------------------------------------------------------- */
    /*              Now read Global Attributes.                             */
    /* -------------------------------------------------------------------- */
    if (poDS->ReadGlobalAttributes(poDS->hSD) != CE_None)
    {
        // Release mutex otherwise we will deadlock with GDALDataset own mutex.
        CPLReleaseMutex(hHDF4Mutex);
        delete poDS;
        CPLAcquireMutex(hHDF4Mutex, 1000.0);
        CPLError(CE_Failure, CPLE_OpenFailed,
                 "Failed to read global attributes from HDF4 file \"%s\".",
                 poOpenInfo->pszFilename);
        return nullptr;
    }

    poDS->SetMetadata(poDS->papszGlobalMetadata, "");

    /* -------------------------------------------------------------------- */
    /*              Determine type of file we read.                         */
    /* -------------------------------------------------------------------- */
    const char *pszValue =
        CSLFetchNameValue(poDS->papszGlobalMetadata, "Signature");

    if (pszValue != nullptr && EQUAL(pszValue, pszGDALSignature))
    {
        poDS->iSubdatasetType = H4ST_GDAL;
        poDS->pszSubdatasetType = "GDAL_HDF4";
    }

    else if ((pszValue = CSLFetchNameValue(poDS->papszGlobalMetadata,
                                           "Title")) != nullptr &&
             EQUAL(pszValue, "SeaWiFS Level-1A Data"))
    {
        poDS->iSubdatasetType = H4ST_SEAWIFS_L1A;
        poDS->pszSubdatasetType = "SEAWIFS_L1A";
    }

    else if ((pszValue = CSLFetchNameValue(poDS->papszGlobalMetadata,
                                           "Title")) != nullptr &&
             EQUAL(pszValue, "SeaWiFS Level-2 Data"))
    {
        poDS->iSubdatasetType = H4ST_SEAWIFS_L2;
        poDS->pszSubdatasetType = "SEAWIFS_L2";
    }

    else if ((pszValue = CSLFetchNameValue(poDS->papszGlobalMetadata,
                                           "Title")) != nullptr &&
             EQUAL(pszValue, "SeaWiFS Level-3 Standard Mapped Image"))
    {
        poDS->iSubdatasetType = H4ST_SEAWIFS_L3;
        poDS->pszSubdatasetType = "SEAWIFS_L3";
    }

    else if ((pszValue = CSLFetchNameValue(poDS->papszGlobalMetadata,
                                           "L1 File Generated By")) !=
                 nullptr &&
             STARTS_WITH_CI(pszValue, "HYP version "))
    {
        poDS->iSubdatasetType = H4ST_HYPERION_L1;
        poDS->pszSubdatasetType = "HYPERION_L1";
    }

    else
    {
        poDS->iSubdatasetType = H4ST_UNKNOWN;
        poDS->pszSubdatasetType = "UNKNOWN";
    }

    /* -------------------------------------------------------------------- */
    /*  If we have HDF-EOS dataset, process it here.                        */
    /* -------------------------------------------------------------------- */
    int32 aiDimSizes[H4_MAX_VAR_DIMS] = {};  // TODO: Get this off of the stack.
    int32 iRank = 0;
    int32 iNumType = 0;
    int32 nAttrs = 0;

    // Sometimes "HDFEOSVersion" attribute is not defined and we will
    // determine HDF-EOS datasets using other records
    // (see ReadGlobalAttributes() method).
    if (poDS->bIsHDFEOS ||
        CSLFetchNameValue(poDS->papszGlobalMetadata, "HDFEOSVersion"))
    {
        /* --------------------------------------------------------------------
         */
        /*  Process swath layers. */
        /* --------------------------------------------------------------------
         */
        hHDF4 = SWopen(poOpenInfo->pszFilename, DFACC_READ);
        if (hHDF4 < 0)
        {
            // Release mutex otherwise we will deadlock with GDALDataset own
            // mutex.
            CPLReleaseMutex(hHDF4Mutex);
            delete poDS;
            CPLAcquireMutex(hHDF4Mutex, 1000.0);
            CPLError(CE_Failure, CPLE_OpenFailed,
                     "Failed to open HDF-EOS file \"%s\" for swath reading.",
                     poOpenInfo->pszFilename);
            return nullptr;
        }
        int32 nStrBufSize = 0;
        int32 nSubDatasets =
            SWinqswath(poOpenInfo->pszFilename, nullptr, &nStrBufSize);

#ifdef DEBUG
        CPLDebug("HDF4", "Number of HDF-EOS swaths: %d",
                 static_cast<int>(nSubDatasets));
#endif

        if (nSubDatasets > 0 && nStrBufSize > 0)
        {
            char *pszSwathList =
                static_cast<char *>(CPLMalloc(nStrBufSize + 1));
            SWinqswath(poOpenInfo->pszFilename, pszSwathList, &nStrBufSize);
            pszSwathList[nStrBufSize] = '\0';

#ifdef DEBUG
            CPLDebug("HDF4", "List of HDF-EOS swaths: %s", pszSwathList);
#endif

            char **papszSwaths =
                CSLTokenizeString2(pszSwathList, ",", CSLT_HONOURSTRINGS);
            CPLFree(pszSwathList);

            if (nSubDatasets != CSLCount(papszSwaths))
            {
                CSLDestroy(papszSwaths);
                // Release mutex otherwise we will deadlock with GDALDataset own
                // mutex.
                CPLReleaseMutex(hHDF4Mutex);
                delete poDS;
                CPLAcquireMutex(hHDF4Mutex, 1000.0);
                CPLDebug("HDF4", "Cannot parse list of HDF-EOS grids.");
                return nullptr;
            }

            for (int32 i = 0; i < nSubDatasets; i++)
            {
                const int32 hSW = SWattach(hHDF4, papszSwaths[i]);

                const int32 nFields =
                    SWnentries(hSW, HDFE_NENTDFLD, &nStrBufSize);
                char *pszFieldList =
                    static_cast<char *>(CPLMalloc(nStrBufSize + 1));
                int32 *paiRank =
                    static_cast<int32 *>(CPLMalloc(nFields * sizeof(int32)));
                int32 *paiNumType =
                    static_cast<int32 *>(CPLMalloc(nFields * sizeof(int32)));

                SWinqdatafields(hSW, pszFieldList, paiRank, paiNumType);

#ifdef DEBUG
                {
                    char *const pszTmp =
                        SPrintArray(GDT_UInt32, paiRank, nFields, ",");

                    CPLDebug("HDF4", "Number of data fields in swath %d: %d",
                             static_cast<int>(i), static_cast<int>(nFields));
                    CPLDebug("HDF4", "List of data fields in swath %d: %s",
                             static_cast<int>(i), pszFieldList);
                    CPLDebug("HDF4", "Data fields ranks: %s", pszTmp);

                    CPLFree(pszTmp);
                }
#endif

                char **papszFields =
                    CSLTokenizeString2(pszFieldList, ",", CSLT_HONOURSTRINGS);

                char szTemp[256] = {'\0'};  // TODO: Get this off the stack.
                for (int32 j = 0; j < nFields; j++)
                {
                    SWfieldinfo(hSW, papszFields[j], &iRank, aiDimSizes,
                                &iNumType, nullptr);

                    if (iRank < 2)
                        continue;

                    // Add field to the list of GDAL subdatasets.
                    const int nCount = CSLCount(poDS->papszSubDatasets) / 2;
                    snprintf(szTemp, sizeof(szTemp), "SUBDATASET_%d_NAME",
                             nCount + 1);
                    // We will use the field index as an identificator.
                    poDS->papszSubDatasets = CSLSetNameValue(
                        poDS->papszSubDatasets, szTemp,
                        CPLSPrintf("HDF4_EOS:EOS_SWATH:\"%s\":%s:%s",
                                   poOpenInfo->pszFilename,
                                   QuoteIfNeeded(papszSwaths[i]).c_str(),
                                   QuoteIfNeeded(papszFields[j]).c_str()));

                    snprintf(szTemp, sizeof(szTemp), "SUBDATASET_%d_DESC",
                             nCount + 1);
                    char *pszString =
                        SPrintArray(GDT_UInt32, aiDimSizes, iRank, "x");
                    poDS->papszSubDatasets = CSLSetNameValue(
                        poDS->papszSubDatasets, szTemp,
                        CPLSPrintf("[%s] %s %s (%s)", pszString, papszFields[j],
                                   papszSwaths[i],
                                   poDS->GetDataTypeName(iNumType)));
                    CPLFree(pszString);
                    szTemp[0] = '\0';
                }

                CSLDestroy(papszFields);
                CPLFree(paiNumType);
                CPLFree(paiRank);
                CPLFree(pszFieldList);
                SWdetach(hSW);
            }

            CSLDestroy(papszSwaths);
        }
        SWclose(hHDF4);

        /* --------------------------------------------------------------------
         */
        /*  Process grid layers. */
        /* --------------------------------------------------------------------
         */
        hHDF4 = GDopen(poOpenInfo->pszFilename, DFACC_READ);
        nSubDatasets =
            GDinqgrid(poOpenInfo->pszFilename, nullptr, &nStrBufSize);

#ifdef DEBUG
        CPLDebug("HDF4", "Number of HDF-EOS grids: %d",
                 static_cast<int>(nSubDatasets));
#endif

        if (nSubDatasets > 0 && nStrBufSize > 0)
        {
            char *pszGridList = static_cast<char *>(CPLMalloc(nStrBufSize + 1));
            GDinqgrid(poOpenInfo->pszFilename, pszGridList, &nStrBufSize);

#ifdef DEBUG
            CPLDebug("HDF4", "List of HDF-EOS grids: %s", pszGridList);
#endif

            char **papszGrids =
                CSLTokenizeString2(pszGridList, ",", CSLT_HONOURSTRINGS);
            CPLFree(pszGridList);

            if (nSubDatasets != CSLCount(papszGrids))
            {
                CSLDestroy(papszGrids);
                GDclose(hHDF4);
                // Release mutex otherwise we will deadlock with GDALDataset own
                // mutex.
                CPLReleaseMutex(hHDF4Mutex);
                delete poDS;
                CPLAcquireMutex(hHDF4Mutex, 1000.0);
                CPLDebug("HDF4", "Cannot parse list of HDF-EOS grids.");
                return nullptr;
            }

            for (int32 i = 0; i < nSubDatasets; i++)
            {
                const int32 hGD = GDattach(hHDF4, papszGrids[i]);

                const int32 nFields =
                    GDnentries(hGD, HDFE_NENTDFLD, &nStrBufSize);
                char *pszFieldList =
                    static_cast<char *>(CPLMalloc(nStrBufSize + 1));
                int32 *paiRank =
                    static_cast<int32 *>(CPLMalloc(nFields * sizeof(int32)));
                int32 *paiNumType =
                    static_cast<int32 *>(CPLMalloc(nFields * sizeof(int32)));

                GDinqfields(hGD, pszFieldList, paiRank, paiNumType);

#ifdef DEBUG
                {
                    char *pszTmp =
                        SPrintArray(GDT_UInt32, paiRank, nFields, ",");
                    CPLDebug("HDF4", "Number of fields in grid %d: %d",
                             static_cast<int>(i), static_cast<int>(nFields));
                    CPLDebug("HDF4", "List of fields in grid %d: %s",
                             static_cast<int>(i), pszFieldList);
                    CPLDebug("HDF4", "Fields ranks: %s", pszTmp);
                    CPLFree(pszTmp);
                }
#endif

                char **papszFields =
                    CSLTokenizeString2(pszFieldList, ",", CSLT_HONOURSTRINGS);

                char szTemp[256];
                for (int32 j = 0; j < nFields; j++)
                {
                    GDfieldinfo(hGD, papszFields[j], &iRank, aiDimSizes,
                                &iNumType, nullptr);

                    if (iRank < 2)
                        continue;

                    // Add field to the list of GDAL subdatasets
                    const int nCount = CSLCount(poDS->papszSubDatasets) / 2;
                    snprintf(szTemp, sizeof(szTemp), "SUBDATASET_%d_NAME",
                             nCount + 1);
                    // We will use the field index as an identificator.
                    poDS->papszSubDatasets = CSLSetNameValue(
                        poDS->papszSubDatasets, szTemp,
                        CPLSPrintf("HDF4_EOS:EOS_GRID:\"%s\":%s:%s",
                                   poOpenInfo->pszFilename,
                                   QuoteIfNeeded(papszGrids[i]).c_str(),
                                   QuoteIfNeeded(papszFields[j]).c_str()));

                    snprintf(szTemp, sizeof(szTemp), "SUBDATASET_%d_DESC",
                             nCount + 1);
                    char *pszString =
                        SPrintArray(GDT_UInt32, aiDimSizes, iRank, "x");
                    poDS->papszSubDatasets = CSLSetNameValue(
                        poDS->papszSubDatasets, szTemp,
                        CPLSPrintf("[%s] %s %s (%s)", pszString, papszFields[j],
                                   papszGrids[i],
                                   poDS->GetDataTypeName(iNumType)));
                    CPLFree(pszString);
                }

                CSLDestroy(papszFields);
                CPLFree(paiNumType);
                CPLFree(paiRank);
                CPLFree(pszFieldList);
                GDdetach(hGD);
            }

            CSLDestroy(papszGrids);
        }
        GDclose(hHDF4);
    }

    char szName[VSNAMELENMAX + 1];

    const char *pszListSDS =
        CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "LIST_SDS", "AUTO");
    if ((poDS->papszSubDatasets == nullptr && EQUAL(pszListSDS, "AUTO")) ||
        (!EQUAL(pszListSDS, "AUTO") && CPLTestBool(pszListSDS)))
    {

        /* --------------------------------------------------------------------
         */
        /*  Make a list of subdatasets from SDSs contained in input HDF file. */
        /* --------------------------------------------------------------------
         */
        int32 nDatasets = 0;

        if (SDfileinfo(poDS->hSD, &nDatasets, &nAttrs) != 0)
            return nullptr;

        char szTemp[256] = {'\0'};  // TODO: Get this off the stack.
        const char *pszName = nullptr;

        for (int32 i = 0; i < nDatasets; i++)
        {
            const int32 iSDS = SDselect(poDS->hSD, i);
            if (SDgetinfo(iSDS, szName, &iRank, aiDimSizes, &iNumType,
                          &nAttrs) != 0)
                return nullptr;

            if (iRank == 1)  // Skip 1D datasets
                continue;

            // Do sort of known datasets. We will display only image bands
            if ((poDS->iSubdatasetType == H4ST_SEAWIFS_L1A) &&
                !STARTS_WITH_CI(szName, "l1a_data"))
                continue;
            else
                pszName = szName;

            // Add datasets with multiple dimensions to the list of GDAL
            // subdatasets.
            const int nCount = CSLCount(poDS->papszSubDatasets) / 2;
            snprintf(szTemp, sizeof(szTemp), "SUBDATASET_%d_NAME", nCount + 1);
            // We will use SDS index as an identificator, because SDS names
            // are not unique. Filename also needed for further file opening
            poDS->papszSubDatasets = CSLSetNameValue(
                poDS->papszSubDatasets, szTemp,
                CPLSPrintf("HDF4_SDS:%s:\"%s\":%ld", poDS->pszSubdatasetType,
                           poOpenInfo->pszFilename, static_cast<long>(i)));
            snprintf(szTemp, sizeof(szTemp), "SUBDATASET_%d_DESC", nCount + 1);
            char *pszString = SPrintArray(GDT_UInt32, aiDimSizes, iRank, "x");
            poDS->papszSubDatasets =
                CSLSetNameValue(poDS->papszSubDatasets, szTemp,
                                CPLSPrintf("[%s] %s (%s)", pszString, pszName,
                                           poDS->GetDataTypeName(iNumType)));
            CPLFree(pszString);

            SDendaccess(iSDS);
            szTemp[0] = '\0';
        }

        SDend(poDS->hSD);
        poDS->hSD = 0;
    }

    /* -------------------------------------------------------------------- */
    /*      Build a list of raster images. Note, that HDF-EOS dataset may   */
    /*      contain a raster image as well.                                 */
    /* -------------------------------------------------------------------- */

    hHDF4 = Hopen(poOpenInfo->pszFilename, DFACC_READ, 0);
    poDS->hGR = GRstart(hHDF4);

    if (poDS->hGR != -1)
    {
        if (GRfileinfo(poDS->hGR, &poDS->nImages, &nAttrs) == -1)
        {
            // Release mutex otherwise we will deadlock with GDALDataset own
            // mutex.
            CPLReleaseMutex(hHDF4Mutex);
            GRend(poDS->hGR);
            poDS->hGR = 0;
            Hclose(hHDF4);
            delete poDS;
            CPLAcquireMutex(hHDF4Mutex, 1000.0);
            return nullptr;
        }

        char szTemp[256] = {'\0'};  // TODO: Get this off the stack.
        for (int32 i = 0; i < poDS->nImages; i++)
        {
            const int32 iGR = GRselect(poDS->hGR, i);

            // iRank in GR interface has another meaning. It represents number
            // of samples per pixel. aiDimSizes has only two dimensions.
            int32 iInterlaceMode = 0;
            if (GRgetiminfo(iGR, szName, &iRank, &iNumType, &iInterlaceMode,
                            aiDimSizes, &nAttrs) != 0)
            {
                // Release mutex otherwise we will deadlock with GDALDataset
                // own mutex.
                CPLReleaseMutex(hHDF4Mutex);
                GRend(poDS->hGR);
                poDS->hGR = 0;
                Hclose(hHDF4);
                delete poDS;
                CPLAcquireMutex(hHDF4Mutex, 1000.0);
                return nullptr;
            }
            const int nCount = CSLCount(poDS->papszSubDatasets) / 2;
            snprintf(szTemp, sizeof(szTemp), "SUBDATASET_%d_NAME", nCount + 1);
            poDS->papszSubDatasets = CSLSetNameValue(
                poDS->papszSubDatasets, szTemp,
                CPLSPrintf("HDF4_GR:UNKNOWN:\"%s\":%ld",
                           poOpenInfo->pszFilename, static_cast<long>(i)));
            snprintf(szTemp, sizeof(szTemp), "SUBDATASET_%d_DESC", nCount + 1);
            char *pszString = SPrintArray(GDT_UInt32, aiDimSizes, 2, "x");
            poDS->papszSubDatasets =
                CSLSetNameValue(poDS->papszSubDatasets, szTemp,
                                CPLSPrintf("[%sx%ld] %s (%s)", pszString,
                                           static_cast<long>(iRank), szName,
                                           poDS->GetDataTypeName(iNumType)));
            CPLFree(pszString);

            GRendaccess(iGR);
            szTemp[0] = '\0';
        }

        GRend(poDS->hGR);
        poDS->hGR = 0;
    }

    Hclose(hHDF4);

    poDS->nRasterXSize = 512;  // XXX: bogus value
    poDS->nRasterYSize = 512;  // XXX: bogus value

    // Make sure we don't try to do any pam stuff with this dataset.
    poDS->nPamFlags |= GPF_NOSAVE;

    /* -------------------------------------------------------------------- */
    /*      If we have single subdataset only, open it immediately          */
    /* -------------------------------------------------------------------- */
    if (CSLCount(poDS->papszSubDatasets) / 2 == 1)
    {
        char *pszSDSName = CPLStrdup(
            CSLFetchNameValue(poDS->papszSubDatasets, "SUBDATASET_1_NAME"));
        // Release mutex otherwise we will deadlock with GDALDataset own mutex.
        CPLReleaseMutex(hHDF4Mutex);
        delete poDS;
        poDS = nullptr;

        GDALDataset *poRetDS =
            GDALDataset::FromHandle(GDALOpen(pszSDSName, poOpenInfo->eAccess));
        CPLFree(pszSDSName);

        CPLAcquireMutex(hHDF4Mutex, 1000.0);

        if (poRetDS)
        {
            poRetDS->SetDescription(poOpenInfo->pszFilename);
        }

        return poRetDS;
    }
    else
    {
        /* --------------------------------------------------------------------
         */
        /*      Confirm the requested access is supported. */
        /* --------------------------------------------------------------------
         */
        if (poOpenInfo->eAccess == GA_Update)
        {
            // Release mutex otherwise we will deadlock with GDALDataset own
            // mutex.
            CPLReleaseMutex(hHDF4Mutex);
            delete poDS;
            CPLAcquireMutex(hHDF4Mutex, 1000.0);

            ReportUpdateNotSupportedByDriver(HDF4_DRIVER_NAME);
            return nullptr;
        }
    }

    return poDS;
}

/************************************************************************/
/*                           HDF4UnloadDriver()                         */
/************************************************************************/

static void HDF4UnloadDriver(GDALDriver * /* poDriver */)
{
    if (hHDF4Mutex != nullptr)
        CPLDestroyMutex(hHDF4Mutex);
    hHDF4Mutex = nullptr;
}

/************************************************************************/
/*                        GDALRegister_HDF4()                           */
/************************************************************************/

void GDALRegister_HDF4()

{
    if (!GDAL_CHECK_VERSION("HDF4 driver"))
        return;

    if (GDALGetDriverByName(HDF4_DRIVER_NAME) != nullptr)
        return;

    GDALDriver *poDriver = new GDALDriver();
    HDF4DriverSetCommonMetadata(poDriver);
    poDriver->pfnOpen = HDF4Dataset::Open;
    poDriver->pfnUnloadDriver = HDF4UnloadDriver;

    GetGDALDriverManager()->RegisterDriver(poDriver);

#ifdef HDF4_HAS_MAXOPENFILES
    poDriver->SetMetadataItem("HDF4_HAS_MAXOPENFILES", "YES");
#endif

#ifdef HDF4_PLUGIN
    GDALRegister_HDF4Image();
#endif
}