text: /******************************************************************************
   *
   * Project:  Erdas Imagine (.img) Translator
   * Purpose:  Supporting functions for HFA (.img) ... main (C callable) API
   *           that is not dependent on GDAL (just CPL).
   * Author:   Frank Warmerdam, warmerdam@pobox.com
   *
   ******************************************************************************
   * Copyright (c) 1999, Intergraph Corporation
   * Copyright (c) 2007-2011, Even Rouault <even dot rouault at spatialys.com>
   *
   * Permission is hereby granted, free of charge, to any person obtaining a
   * copy of this software and associated documentation files (the "Software"),
   * to deal in the Software without restriction, including without limitation
   * the rights to use, copy, modify, merge, publish, distribute, sublicense,
   * and/or sell copies of the Software, and to permit persons to whom the
   * Software is furnished to do so, subject to the following conditions:
   *
   * The above copyright notice and this permission notice shall be included
   * in all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
   * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
   * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
   * DEALINGS IN THE SOFTWARE.
   ******************************************************************************
   *
   * hfaopen.cpp
   *
   * Supporting routines for reading Erdas Imagine (.imf) Hierarchical
   * File Architecture files.  This is intended to be a library independent
   * of the GDAL core, but dependent on the Common Portability Library.
   *
   */
  
  #include "cpl_port.h"
  #include "hfa_p.h"
  
  #include <cerrno>
  #include <climits>
  #include <cmath>
  #include <cstddef>
  #include <cstdio>
  #include <cstdlib>
  #include <cstring>
  #if HAVE_FCNTL_H
  #  include <fcntl.h>
  #endif
  #include <algorithm>
  #include <memory>
  #include <string>
  #include <vector>
  
  #include "cpl_conv.h"
  #include "cpl_error.h"
  #include "cpl_string.h"
  #include "cpl_vsi.h"
  #include "gdal_priv.h"
  #include "hfa.h"
  
  CPL_CVSID("$Id$")
  
  static const char * const apszAuxMetadataItems[] = {
  // node/entry            field_name                  metadata_key       type
   "Statistics",           "dminimum",              "STATISTICS_MINIMUM",     "Esta_Statistics",
   "Statistics",           "dmaximum",              "STATISTICS_MAXIMUM",     "Esta_Statistics",
   "Statistics",           "dmean",                 "STATISTICS_MEAN",        "Esta_Statistics",
   "Statistics",           "dmedian",               "STATISTICS_MEDIAN",      "Esta_Statistics",
   "Statistics",           "dmode",                 "STATISTICS_MODE",        "Esta_Statistics",
   "Statistics",           "dstddev",               "STATISTICS_STDDEV",      "Esta_Statistics",
   "HistogramParameters",  "lBinFunction.numBins",  "STATISTICS_HISTONUMBINS","Eimg_StatisticsParameters830",
   "HistogramParameters",  "dBinFunction.minLimit", "STATISTICS_HISTOMIN",    "Eimg_StatisticsParameters830",
   "HistogramParameters",  "dBinFunction.maxLimit", "STATISTICS_HISTOMAX",    "Eimg_StatisticsParameters830",
   "StatisticsParameters", "lSkipFactorX",          "STATISTICS_SKIPFACTORX", "",
   "StatisticsParameters", "lSkipFactorY",          "STATISTICS_SKIPFACTORY", "",
   "StatisticsParameters", "dExcludedValues",       "STATISTICS_EXCLUDEDVALUES","",
   "",                     "elayerType",            "LAYER_TYPE",             "",
   "RRDInfoList",          "salgorithm.string",     "OVERVIEWS_ALGORITHM",    "Emif_String",
   nullptr
  };
  
  const char * const * GetHFAAuxMetaDataList()
  {
      return apszAuxMetadataItems;
  }
  
  /************************************************************************/
  /*                          HFAGetDictionary()                          */
  /************************************************************************/
  
  static char * HFAGetDictionary( HFAHandle hHFA )
  
  {
      int nDictMax = 100;
      char *pszDictionary = static_cast<char *>(CPLMalloc(nDictMax));
      int nDictSize = 0;
  
      if( VSIFSeekL(hHFA->fp, hHFA->nDictionaryPos, SEEK_SET) < 0 )
      {
          pszDictionary[nDictSize] = '\0';
          return pszDictionary;
      }
  
      while( true )
      {
          if( nDictSize >= nDictMax - 1 )
          {
              nDictMax = nDictSize * 2 + 100;
              pszDictionary =
                  static_cast<char *>(CPLRealloc(pszDictionary, nDictMax));
          }
  
          if( VSIFReadL(pszDictionary + nDictSize, 1, 1, hHFA->fp) < 1 ||
              pszDictionary[nDictSize] == '\0' ||
              (nDictSize > 2 && pszDictionary[nDictSize - 2] == ',' &&
               pszDictionary[nDictSize-1] == '.') )
              break;
  
          nDictSize++;
      }
  
      pszDictionary[nDictSize] = '\0';
  
      return pszDictionary;
  }
  
  /************************************************************************/
  /*                              HFAOpen()                               */
  /************************************************************************/
  
  HFAHandle HFAOpen( const char *pszFilename, const char *pszAccess )
  
  {
      VSILFILE *fp = VSIFOpenL(
          pszFilename,
          (EQUAL(pszAccess, "r") || EQUAL(pszAccess, "rb")) ? "rb" : "r+b");
  
      // Should this be changed to use some sort of CPLFOpen() which will
      // set the error?
      if( fp == nullptr )
      {
          CPLError(CE_Failure, CPLE_OpenFailed, "File open of %s failed.",
                   pszFilename);
  
          return nullptr;
      }
  
      // Read and verify the header.
      char szHeader[16] = {};
      if( VSIFReadL(szHeader, 16, 1, fp) < 1 )
      {
          CPLError(CE_Failure, CPLE_AppDefined,
                   "Attempt to read 16 byte header failed for\n%s.", pszFilename);
          CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
          return nullptr;
      }
  
      if( !STARTS_WITH_CI(szHeader, "EHFA_HEADER_TAG") )
      {
          CPLError(CE_Failure, CPLE_AppDefined,
                   "File %s is not an Imagine HFA file ... header wrong.",
                   pszFilename);
          CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
          return nullptr;
      }
  
      // Create the HFAInfo_t.
      HFAInfo_t *psInfo =
          static_cast<HFAInfo_t *>(CPLCalloc(sizeof(HFAInfo_t), 1));
  
      psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
      psInfo->pszPath = CPLStrdup(CPLGetPath(pszFilename));
      psInfo->fp = fp;
      if( EQUAL(pszAccess, "r") || EQUAL(pszAccess, "rb") )
          psInfo->eAccess = HFA_ReadOnly;
      else
          psInfo->eAccess = HFA_Update;
      psInfo->bTreeDirty = false;
  
      // Where is the header?
      GUInt32 nHeaderPos = 0;
      bool bRet = VSIFReadL(&nHeaderPos, sizeof(GInt32), 1, fp) > 0;
      HFAStandard(4, &nHeaderPos);
  
      // Read the header.
      bRet &= VSIFSeekL(fp, nHeaderPos, SEEK_SET) >= 0;
  
      bRet &= VSIFReadL(&(psInfo->nVersion), sizeof(GInt32), 1, fp) > 0;
      HFAStandard(4, &(psInfo->nVersion));
  
      bRet &= VSIFReadL(szHeader, 4, 1, fp) > 0;  // Skip freeList.
  
      bRet &= VSIFReadL(&(psInfo->nRootPos), sizeof(GInt32), 1, fp) > 0;
      HFAStandard(4, &(psInfo->nRootPos));
  
      bRet &= VSIFReadL(&(psInfo->nEntryHeaderLength), sizeof(GInt16), 1, fp) > 0;
      HFAStandard(2, &(psInfo->nEntryHeaderLength));
  
      bRet &= VSIFReadL(&(psInfo->nDictionaryPos), sizeof(GInt32), 1, fp) > 0;
      HFAStandard(4, &(psInfo->nDictionaryPos));
  
      // Collect file size.
      bRet &= VSIFSeekL(fp, 0, SEEK_END) >= 0;
      if( !bRet )
      {
          CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
          CPLFree(psInfo->pszFilename);
          CPLFree(psInfo->pszPath);
          CPLFree(psInfo);
          return nullptr;
      }
      psInfo->nEndOfFile = static_cast<GUInt32>(VSIFTellL(fp));
  
      // Instantiate the root entry.
      psInfo->poRoot = HFAEntry::New(psInfo, psInfo->nRootPos, nullptr, nullptr);
      if( psInfo->poRoot == nullptr )
      {
          CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
          CPLFree(psInfo->pszFilename);
          CPLFree(psInfo->pszPath);
          CPLFree(psInfo);
          return nullptr;
      }
  
      // Read the dictionary.
      psInfo->pszDictionary = HFAGetDictionary(psInfo);
      psInfo->poDictionary = new HFADictionary(psInfo->pszDictionary);
  
      // Collect band definitions.
      HFAParseBandInfo(psInfo);
  
      return psInfo;
  }
  
  /************************************************************************/
  /*                         HFACreateDependent()                         */
  /*                                                                      */
  /*      Create a .rrd file for the named file if it does not exist,     */
  /*      or return the existing dependent if it already exists.          */
  /************************************************************************/
  
  HFAInfo_t *HFACreateDependent( HFAInfo_t *psBase )
  
  {
      if( psBase->psDependent != nullptr )
          return psBase->psDependent;
  
      // Create desired RRD filename.
      const CPLString oBasename = CPLGetBasename(psBase->pszFilename);
      const CPLString oRRDFilename =
          CPLFormFilename(psBase->pszPath, oBasename, "rrd");
  
      // Does this file already exist?  If so, re-use it.
      VSILFILE *fp = VSIFOpenL(oRRDFilename, "rb");
      if( fp != nullptr )
      {
          CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
          psBase->psDependent = HFAOpen(oRRDFilename, "rb");
          // FIXME? this is not going to be reused but recreated...
      }
  
      // Otherwise create it now.
      HFAInfo_t *psDep = HFACreateLL(oRRDFilename);
      psBase->psDependent = psDep;
      if( psDep == nullptr )
          return nullptr;
  
  /* -------------------------------------------------------------------- */
  /*      Add the DependentFile node with the pointer back to the         */
  /*      parent.  When working from an .aux file we really want the      */
  /*      .rrd to point back to the original file, not the .aux file.     */
  /* -------------------------------------------------------------------- */
      HFAEntry *poEntry = psBase->poRoot->GetNamedChild("DependentFile");
      const char *pszDependentFile = nullptr;
      if( poEntry != nullptr )
          pszDependentFile = poEntry->GetStringField("dependent.string");
      if( pszDependentFile == nullptr )
          pszDependentFile = psBase->pszFilename;
  
      HFAEntry *poDF = HFAEntry::New(psDep, "DependentFile",
                                     "Eimg_DependentFile", psDep->poRoot);
  
      poDF->MakeData(static_cast<int>(strlen(pszDependentFile) + 50));
      poDF->SetPosition();
      poDF->SetStringField("dependent.string", pszDependentFile);
  
      return psDep;
  }
  
  /************************************************************************/
  /*                          HFAGetDependent()                           */
  /************************************************************************/
  
  HFAInfo_t *HFAGetDependent(HFAInfo_t *psBase, const char *pszFilename)
  
  {
      if( EQUAL(pszFilename, psBase->pszFilename) )
          return psBase;
  
      if( psBase->psDependent != nullptr )
      {
          if( EQUAL(pszFilename, psBase->psDependent->pszFilename) )
              return psBase->psDependent;
          else
              return nullptr;
      }
  
      // Try to open the dependent file.
      const char *pszMode = psBase->eAccess == HFA_Update ? "r+b" : "rb";
  
      char *pszDependent =
          CPLStrdup(CPLFormFilename(psBase->pszPath, pszFilename, nullptr));
  
      VSILFILE *fp = VSIFOpenL(pszDependent, pszMode);
      if( fp != nullptr )
      {
          CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
          psBase->psDependent = HFAOpen(pszDependent, pszMode);
      }
  
      CPLFree(pszDependent);
  
      return psBase->psDependent;
  }
  
  /************************************************************************/
  /*                          HFAParseBandInfo()                          */
  /*                                                                      */
  /*      This is used by HFAOpen() and HFACreate() to initialize the     */
  /*      band structures.                                                */
  /************************************************************************/
  
  CPLErr HFAParseBandInfo( HFAInfo_t *psInfo )
  
  {
      // Find the first band node.
      psInfo->nBands = 0;
      HFAEntry *poNode = psInfo->poRoot->GetChild();
      while( poNode != nullptr )
      {
          if( EQUAL(poNode->GetType(), "Eimg_Layer") &&
              poNode->GetIntField("width") > 0 &&
              poNode->GetIntField("height") > 0 )
          {
              if( psInfo->nBands == 0 )
              {
                  psInfo->nXSize = poNode->GetIntField("width");
                  psInfo->nYSize = poNode->GetIntField("height");
              }
              else if( poNode->GetIntField("width") != psInfo->nXSize ||
                       poNode->GetIntField("height") != psInfo->nYSize )
              {
                  return CE_Failure;
              }
  
              psInfo->papoBand = static_cast<HFABand **>(CPLRealloc(
                  psInfo->papoBand, sizeof(HFABand *) * (psInfo->nBands + 1)));
              psInfo->papoBand[psInfo->nBands] = new HFABand(psInfo, poNode);
              if( psInfo->papoBand[psInfo->nBands]->nWidth == 0 )
              {
                  delete psInfo->papoBand[psInfo->nBands];
                  return CE_Failure;
              }
              psInfo->nBands++;
          }
  
          poNode = poNode->GetNext();
      }
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                              HFAClose()                              */
  /************************************************************************/
  
  int HFAClose( HFAHandle hHFA )
  
  {
      if( hHFA->eAccess == HFA_Update && (hHFA->bTreeDirty ||
          (hHFA->poDictionary != nullptr &&
           hHFA->poDictionary->bDictionaryTextDirty)) )
          HFAFlush(hHFA);
  
      int nRet = 0;
      if( hHFA->psDependent != nullptr )
      {
          if( HFAClose(hHFA->psDependent) != 0 )
              nRet = -1;
      }
  
      delete hHFA->poRoot;
  
      if( VSIFCloseL(hHFA->fp) != 0 )
          nRet = -1;
  
      if( hHFA->poDictionary != nullptr )
          delete hHFA->poDictionary;
  
      CPLFree(hHFA->pszDictionary);
      CPLFree(hHFA->pszFilename);
      CPLFree(hHFA->pszIGEFilename);
      CPLFree(hHFA->pszPath);
  
      for( int i = 0; i < hHFA->nBands; i++ )
      {
          delete hHFA->papoBand[i];
      }
  
      CPLFree(hHFA->papoBand);
  
      if( hHFA->pProParameters != nullptr )
      {
          Eprj_ProParameters *psProParms = (Eprj_ProParameters *)
              hHFA->pProParameters;
  
          CPLFree(psProParms->proExeName);
          CPLFree(psProParms->proName);
          CPLFree(psProParms->proSpheroid.sphereName);
  
          CPLFree(psProParms);
      }
  
      if( hHFA->pDatum != nullptr )
      {
          CPLFree(((Eprj_Datum *)hHFA->pDatum)->datumname);
          CPLFree(((Eprj_Datum *)hHFA->pDatum)->gridname);
          CPLFree(hHFA->pDatum);
      }
  
      if( hHFA->pMapInfo != nullptr )
      {
          CPLFree(((Eprj_MapInfo *)hHFA->pMapInfo)->proName);
          CPLFree(((Eprj_MapInfo *)hHFA->pMapInfo)->units);
          CPLFree(hHFA->pMapInfo);
      }
  
      CPLFree(hHFA);
      return nRet;
  }
  
  /************************************************************************/
  /*                              HFARemove()                             */
  /*  Used from HFADelete() function.                                     */
  /************************************************************************/
  
  static CPLErr HFARemove( const char *pszFilename )
  
  {
      VSIStatBufL sStat;
  
      if( VSIStatL(pszFilename, &sStat) == 0 && VSI_ISREG(sStat.st_mode) )
      {
          if( VSIUnlink(pszFilename) == 0 )
              return CE_None;
          else
          {
              CPLError(CE_Failure, CPLE_AppDefined,
                       "Attempt to unlink %s failed.", pszFilename);
              return CE_Failure;
          }
      }
  
      CPLError(CE_Failure, CPLE_AppDefined, "Unable to delete %s, not a file.",
               pszFilename);
      return CE_Failure;
  }
  
  /************************************************************************/
  /*                              HFADelete()                             */
  /************************************************************************/
  
  CPLErr HFADelete( const char *pszFilename )
  
  {
      HFAInfo_t *psInfo = HFAOpen(pszFilename, "rb");
      HFAEntry *poDMS = nullptr;
      HFAEntry *poLayer = nullptr;
      HFAEntry *poNode = nullptr;
  
      if( psInfo != nullptr )
      {
          poNode = psInfo->poRoot->GetChild();
          while( (poNode != nullptr) && (poLayer == nullptr) )
          {
              if( EQUAL(poNode->GetType(), "Eimg_Layer") )
              {
                  poLayer = poNode;
              }
              poNode = poNode->GetNext();
          }
  
          if( poLayer != nullptr )
              poDMS = poLayer->GetNamedChild("ExternalRasterDMS");
  
          if( poDMS )
          {
              const char *pszRawFilename =
                  poDMS->GetStringField("fileName.string");
  
              if( pszRawFilename != nullptr )
                  HFARemove(
                      CPLFormFilename(psInfo->pszPath, pszRawFilename, nullptr));
          }
  
          CPL_IGNORE_RET_VAL(HFAClose(psInfo));
      }
      return HFARemove(pszFilename);
  }
  
  /************************************************************************/
  /*                          HFAGetRasterInfo()                          */
  /************************************************************************/
  
  CPLErr HFAGetRasterInfo( HFAHandle hHFA, int *pnXSize, int *pnYSize,
                           int *pnBands )
  
  {
      if( pnXSize != nullptr )
          *pnXSize = hHFA->nXSize;
      if( pnYSize != nullptr )
          *pnYSize = hHFA->nYSize;
      if( pnBands != nullptr )
          *pnBands = hHFA->nBands;
      return CE_None;
  }
  
  /************************************************************************/
  /*                           HFAGetBandInfo()                           */
  /************************************************************************/
  
  CPLErr HFAGetBandInfo( HFAHandle hHFA, int nBand, EPTType *peDataType,
                         int *pnBlockXSize, int *pnBlockYSize,
                         int *pnCompressionType )
  
  {
      if( nBand < 0 || nBand > hHFA->nBands )
      {
          CPLAssert(false);
          return CE_Failure;
      }
  
      HFABand *poBand = hHFA->papoBand[nBand - 1];
  
      if( peDataType != nullptr )
          *peDataType = poBand->eDataType;
  
      if( pnBlockXSize != nullptr )
          *pnBlockXSize = poBand->nBlockXSize;
  
      if( pnBlockYSize != nullptr )
          *pnBlockYSize = poBand->nBlockYSize;
  
      // Get compression code from RasterDMS.
      if( pnCompressionType != nullptr )
      {
          *pnCompressionType = 0;
  
          HFAEntry *poDMS = poBand->poNode->GetNamedChild("RasterDMS");
  
          if( poDMS != nullptr )
              *pnCompressionType = poDMS->GetIntField("compressionType");
      }
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                          HFAGetBandNoData()                          */
  /*                                                                      */
  /*      returns TRUE if value is set, otherwise FALSE.                  */
  /************************************************************************/
  
  int HFAGetBandNoData( HFAHandle hHFA, int nBand, double *pdfNoData )
  
  {
      if( nBand < 0 || nBand > hHFA->nBands )
      {
          CPLAssert(false);
          return CE_Failure;
      }
  
      HFABand *poBand = hHFA->papoBand[nBand - 1];
  
      if( !poBand->bNoDataSet && poBand->nOverviews > 0 )
      {
          poBand = poBand->papoOverviews[0];
          if( poBand == nullptr )
              return FALSE;
      }
  
      *pdfNoData = poBand->dfNoData;
      return poBand->bNoDataSet;
  }
  
  /************************************************************************/
  /*                          HFASetBandNoData()                          */
  /*                                                                      */
  /*      attempts to set a no-data value on the given band               */
  /************************************************************************/
  
  CPLErr HFASetBandNoData( HFAHandle hHFA, int nBand, double dfValue )
  
  {
      if( nBand < 0 || nBand > hHFA->nBands )
      {
          CPLAssert(false);
          return CE_Failure;
      }
  
      HFABand *poBand = hHFA->papoBand[nBand - 1];
  
      return poBand->SetNoDataValue( dfValue );
  }
  
  /************************************************************************/
  /*                        HFAGetOverviewCount()                         */
  /************************************************************************/
  
  int HFAGetOverviewCount( HFAHandle hHFA, int nBand )
  
  {
      if( nBand < 0 || nBand > hHFA->nBands )
      {
          CPLAssert(false);
          return CE_Failure;
      }
  
      HFABand *poBand = hHFA->papoBand[nBand - 1];
      poBand->LoadOverviews();
  
      return poBand->nOverviews;
  }
  
  /************************************************************************/
  /*                         HFAGetOverviewInfo()                         */
  /************************************************************************/
  
  CPLErr HFAGetOverviewInfo( HFAHandle hHFA, int nBand, int iOverview,
                             int *pnXSize, int *pnYSize,
                             int *pnBlockXSize, int *pnBlockYSize,
                             EPTType *peHFADataType )
  
  {
      if( nBand < 0 || nBand > hHFA->nBands )
      {
          CPLAssert(false);
          return CE_Failure;
      }
  
      HFABand *poBand = hHFA->papoBand[nBand - 1];
      poBand->LoadOverviews();
  
      if( iOverview < 0 || iOverview >= poBand->nOverviews )
      {
          CPLAssert(false);
          return CE_Failure;
      }
      poBand = poBand->papoOverviews[iOverview];
      if( poBand == nullptr )
      {
          return CE_Failure;
      }
  
      if( pnXSize != nullptr )
          *pnXSize = poBand->nWidth;
  
      if( pnYSize != nullptr )
          *pnYSize = poBand->nHeight;
  
      if( pnBlockXSize != nullptr )
          *pnBlockXSize = poBand->nBlockXSize;
  
      if( pnBlockYSize != nullptr )
          *pnBlockYSize = poBand->nBlockYSize;
  
      if( peHFADataType != nullptr )
          *peHFADataType = poBand->eDataType;
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                         HFAGetRasterBlock()                          */
  /************************************************************************/
  
  CPLErr HFAGetRasterBlock( HFAHandle hHFA, int nBand,
                            int nXBlock, int nYBlock, void *pData )
  
  {
      return HFAGetRasterBlockEx(hHFA, nBand, nXBlock, nYBlock, pData, -1);
  }
  
  /************************************************************************/
  /*                        HFAGetRasterBlockEx()                         */
  /************************************************************************/
  
  CPLErr HFAGetRasterBlockEx( HFAHandle hHFA, int nBand,
                              int nXBlock, int nYBlock, void *pData,
                              int nDataSize )
  
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return CE_Failure;
  
      return hHFA->papoBand[nBand - 1]->
          GetRasterBlock(nXBlock, nYBlock, pData,nDataSize);
  }
  
  /************************************************************************/
  /*                     HFAGetOverviewRasterBlock()                      */
  /************************************************************************/
  
  CPLErr HFAGetOverviewRasterBlock( HFAHandle hHFA, int nBand, int iOverview,
                                    int nXBlock, int nYBlock, void *pData )
  
  {
      return HFAGetOverviewRasterBlockEx(hHFA, nBand, iOverview, nXBlock, nYBlock,
                                         pData, -1);
  }
  
  /************************************************************************/
  /*                   HFAGetOverviewRasterBlockEx()                      */
  /************************************************************************/
  
  CPLErr HFAGetOverviewRasterBlockEx( HFAHandle hHFA, int nBand, int iOverview,
                                      int nXBlock, int nYBlock, void *pData,
                                      int nDataSize )
  
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return CE_Failure;
  
      if( iOverview < 0 || iOverview >= hHFA->papoBand[nBand - 1]->nOverviews )
          return CE_Failure;
  
      return hHFA->papoBand[nBand - 1]->papoOverviews[iOverview]->
          GetRasterBlock(nXBlock, nYBlock, pData, nDataSize);
  }
  
  /************************************************************************/
  /*                         HFASetRasterBlock()                          */
  /************************************************************************/
  
  CPLErr HFASetRasterBlock( HFAHandle hHFA, int nBand,
                            int nXBlock, int nYBlock, void *pData )
  
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return CE_Failure;
  
      return hHFA->papoBand[nBand - 1]->SetRasterBlock(nXBlock, nYBlock, pData);
  }
  
  /************************************************************************/
  /*                         HFASetRasterBlock()                          */
  /************************************************************************/
  
  CPLErr HFASetOverviewRasterBlock( HFAHandle hHFA, int nBand, int iOverview,
                                    int nXBlock, int nYBlock, void *pData )
  
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return CE_Failure;
  
      if( iOverview < 0 || iOverview >= hHFA->papoBand[nBand - 1]->nOverviews )
          return CE_Failure;
  
      return hHFA->papoBand[nBand - 1]->papoOverviews[iOverview]->
          SetRasterBlock(nXBlock, nYBlock, pData);
  }
  
  /************************************************************************/
  /*                         HFAGetBandName()                             */
  /************************************************************************/
  
  const char *HFAGetBandName( HFAHandle hHFA, int nBand )
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return "";
  
      return hHFA->papoBand[nBand - 1]->GetBandName();
  }
  
  /************************************************************************/
  /*                         HFASetBandName()                             */
  /************************************************************************/
  
  void HFASetBandName( HFAHandle hHFA, int nBand, const char *pszName )
  {
      if( nBand < 1 || nBand > hHFA->nBands )
        return;
  
      hHFA->papoBand[nBand - 1]->SetBandName(pszName);
  }
  
  /************************************************************************/
  /*                         HFAGetDataTypeBits()                         */
  /************************************************************************/
  
  int HFAGetDataTypeBits( EPTType eDataType )
  
  {
      switch( eDataType )
      {
      case EPT_u1:
          return 1;
  
      case EPT_u2:
          return 2;
  
      case EPT_u4:
          return 4;
  
      case EPT_u8:
      case EPT_s8:
          return 8;
  
      case EPT_u16:
      case EPT_s16:
          return 16;
  
      case EPT_u32:
      case EPT_s32:
      case EPT_f32:
          return 32;
  
      case EPT_f64:
      case EPT_c64:
          return 64;
  
      case EPT_c128:
          return 128;
      }
  
      CPLAssert(false);
      return 1;
  }
  
  /************************************************************************/
  /*                         HFAGetDataTypeName()                         */
  /************************************************************************/
  
  const char *HFAGetDataTypeName( EPTType eDataType )
  
  {
      switch( eDataType )
      {
      case EPT_u1:
          return "u1";
  
      case EPT_u2:
          return "u2";
  
      case EPT_u4:
          return "u4";
  
      case EPT_u8:
          return "u8";
  
      case EPT_s8:
          return "s8";
  
      case EPT_u16:
          return "u16";
  
      case EPT_s16:
          return "s16";
  
      case EPT_u32:
          return "u32";
  
      case EPT_s32:
          return "s32";
  
      case EPT_f32:
          return "f32";
  
      case EPT_f64:
          return "f64";
  
      case EPT_c64:
          return "c64";
  
      case EPT_c128:
          return "c128";
  
      default:
          CPLAssert(false);
          return "unknown";
      }
  }
  
  /************************************************************************/
  /*                           HFAGetMapInfo()                            */
  /************************************************************************/
  
  const Eprj_MapInfo *HFAGetMapInfo( HFAHandle hHFA )
  
  {
      if( hHFA->nBands < 1 )
          return nullptr;
  
      // Do we already have it?
      if( hHFA->pMapInfo != nullptr )
          return (Eprj_MapInfo *)hHFA->pMapInfo;
  
      // Get the HFA node.  If we don't find it under the usual name
      // we search for any node of the right type (#3338).
      HFAEntry *poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild("Map_Info");
      if( poMIEntry == nullptr )
      {
          for( HFAEntry *poChild = hHFA->papoBand[0]->poNode->GetChild();
               poChild != nullptr && poMIEntry == nullptr;
               poChild = poChild->GetNext() )
          {
              if( EQUAL(poChild->GetType(), "Eprj_MapInfo") )
                  poMIEntry = poChild;
          }
      }
  
      if( poMIEntry == nullptr )
      {
          return nullptr;
      }
  
      // Allocate the structure.
      Eprj_MapInfo *psMapInfo =
          static_cast<Eprj_MapInfo *>(CPLCalloc(sizeof(Eprj_MapInfo), 1));
  
      // Fetch the fields.
      psMapInfo->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
  
      psMapInfo->upperLeftCenter.x =
          poMIEntry->GetDoubleField("upperLeftCenter.x");
      psMapInfo->upperLeftCenter.y =
          poMIEntry->GetDoubleField("upperLeftCenter.y");
  
      psMapInfo->lowerRightCenter.x =
          poMIEntry->GetDoubleField("lowerRightCenter.x");
      psMapInfo->lowerRightCenter.y =
          poMIEntry->GetDoubleField("lowerRightCenter.y");
  
      CPLErr eErr = CE_None;
      psMapInfo->pixelSize.width =
          poMIEntry->GetDoubleField("pixelSize.width", &eErr);
      psMapInfo->pixelSize.height =
          poMIEntry->GetDoubleField("pixelSize.height", &eErr);
  
      // The following is basically a hack to get files with
      // non-standard MapInfo's that misname the pixelSize fields. (#3338)
      if( eErr != CE_None )
      {
          psMapInfo->pixelSize.width = poMIEntry->GetDoubleField("pixelSize.x");
          psMapInfo->pixelSize.height = poMIEntry->GetDoubleField("pixelSize.y");
      }
  
      psMapInfo->units = CPLStrdup(poMIEntry->GetStringField("units"));
  
      hHFA->pMapInfo = (void *)psMapInfo;
  
      return psMapInfo;
  }
  
  /************************************************************************/
  /*                        HFAInvGeoTransform()                          */
  /************************************************************************/
  
  static bool HFAInvGeoTransform( const double *gt_in, double *gt_out )
  
  {
      // Assume a 3rd row that is [1 0 0].
      // Compute determinate.
      const double det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
  
      if( fabs(det) < 1.0e-15 )
          return false;
  
      const double inv_det = 1.0 / det;
  
      // Compute adjoint, and divide by determinate.
      gt_out[1] =  gt_in[5] * inv_det;
      gt_out[4] = -gt_in[4] * inv_det;
  
      gt_out[2] = -gt_in[2] * inv_det;
      gt_out[5] =  gt_in[1] * inv_det;
  
      gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
      gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
  
      return true;
  }
  
  /************************************************************************/
  /*                         HFAGetGeoTransform()                         */
  /************************************************************************/
  
  int HFAGetGeoTransform( HFAHandle hHFA, double *padfGeoTransform )
  
  {
      const Eprj_MapInfo *psMapInfo = HFAGetMapInfo(hHFA);
  
      padfGeoTransform[0] = 0.0;
      padfGeoTransform[1] = 1.0;
      padfGeoTransform[2] = 0.0;
      padfGeoTransform[3] = 0.0;
      padfGeoTransform[4] = 0.0;
      padfGeoTransform[5] = 1.0;
  
      // Simple (north up) MapInfo approach.
      if( psMapInfo != nullptr )
      {
          padfGeoTransform[0] =
              psMapInfo->upperLeftCenter.x - psMapInfo->pixelSize.width * 0.5;
          padfGeoTransform[1] = psMapInfo->pixelSize.width;
          if( padfGeoTransform[1] == 0.0 )
              padfGeoTransform[1] = 1.0;
          padfGeoTransform[2] = 0.0;
          if( psMapInfo->upperLeftCenter.y >= psMapInfo->lowerRightCenter.y )
              padfGeoTransform[5] = -psMapInfo->pixelSize.height;
          else
              padfGeoTransform[5] = psMapInfo->pixelSize.height;
          if( padfGeoTransform[5] == 0.0 )
              padfGeoTransform[5] = 1.0;
  
          padfGeoTransform[3] =
              psMapInfo->upperLeftCenter.y - padfGeoTransform[5] * 0.5;
          padfGeoTransform[4] = 0.0;
  
          // Special logic to fixup odd angular units.
          if( EQUAL(psMapInfo->units, "ds") )
          {
              padfGeoTransform[0] /= 3600.0;
              padfGeoTransform[1] /= 3600.0;
              padfGeoTransform[2] /= 3600.0;
              padfGeoTransform[3] /= 3600.0;
              padfGeoTransform[4] /= 3600.0;
              padfGeoTransform[5] /= 3600.0;
          }
  
          return TRUE;
      }
  
      // Try for a MapToPixelXForm affine polynomial supporting
      // rotated and sheared affine transformations.
      if( hHFA->nBands == 0 )
          return FALSE;
  
      HFAEntry *poXForm0 =
          hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm0");
  
      if( poXForm0 == nullptr )
          return FALSE;
  
      if( poXForm0->GetIntField("order") != 1 ||
          poXForm0->GetIntField("numdimtransform") != 2 ||
          poXForm0->GetIntField("numdimpolynomial") != 2 ||
          poXForm0->GetIntField("termcount") != 3 )
          return FALSE;
  
      // Verify that there aren't any further xform steps.
      if( hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm1")
          != nullptr )
          return FALSE;
  
      // We should check that the exponent list is 0 0 1 0 0 1, but
      // we don't because we are lazy.
  
      // Fetch geotransform values.
      double adfXForm[6] = {
          poXForm0->GetDoubleField("polycoefvector[0]"),
          poXForm0->GetDoubleField("polycoefmtx[0]"),
          poXForm0->GetDoubleField("polycoefmtx[2]"),
          poXForm0->GetDoubleField("polycoefvector[1]"),
          poXForm0->GetDoubleField("polycoefmtx[1]"),
          poXForm0->GetDoubleField("polycoefmtx[3]")
      };
  
      // Invert.
  
      if( !HFAInvGeoTransform(adfXForm, padfGeoTransform) )
          memset(padfGeoTransform, 0, 6 * sizeof(double));
  
      // Adjust origin from center of top left pixel to top left corner
      // of top left pixel.
      padfGeoTransform[0] -= padfGeoTransform[1] * 0.5;
      padfGeoTransform[0] -= padfGeoTransform[2] * 0.5;
      padfGeoTransform[3] -= padfGeoTransform[4] * 0.5;
      padfGeoTransform[3] -= padfGeoTransform[5] * 0.5;
  
      return TRUE;
  }
  
  /************************************************************************/
  /*                           HFASetMapInfo()                            */
  /************************************************************************/
  
  CPLErr HFASetMapInfo( HFAHandle hHFA, const Eprj_MapInfo *poMapInfo )
  
  {
      // Loop over bands, setting information on each one.
      for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
      {
          // Create a new Map_Info if there isn't one present already.
          HFAEntry *poMIEntry =
              hHFA->papoBand[iBand]->poNode->GetNamedChild("Map_Info");
          if( poMIEntry == nullptr )
          {
              poMIEntry = HFAEntry::New(hHFA, "Map_Info", "Eprj_MapInfo",
                                        hHFA->papoBand[iBand]->poNode);
          }
  
          poMIEntry->MarkDirty();
  
          // Ensure we have enough space for all the data.
          // TODO(schwehr): Explain 48 and 40 constants.
          const int nSize = static_cast<int>(
              48 + 40
              + strlen(poMapInfo->proName) + 1
              + strlen(poMapInfo->units) + 1);
  
          GByte *pabyData = poMIEntry->MakeData(nSize);
          memset(pabyData, 0, nSize);
  
          poMIEntry->SetPosition();
  
          // Write the various fields.
          poMIEntry->SetStringField("proName", poMapInfo->proName);
  
          poMIEntry->SetDoubleField("upperLeftCenter.x",
                                    poMapInfo->upperLeftCenter.x);
          poMIEntry->SetDoubleField("upperLeftCenter.y",
                                    poMapInfo->upperLeftCenter.y);
  
          poMIEntry->SetDoubleField("lowerRightCenter.x",
                                    poMapInfo->lowerRightCenter.x);
          poMIEntry->SetDoubleField("lowerRightCenter.y",
                                    poMapInfo->lowerRightCenter.y);
  
          poMIEntry->SetDoubleField("pixelSize.width",
                                    poMapInfo->pixelSize.width);
          poMIEntry->SetDoubleField("pixelSize.height",
                                    poMapInfo->pixelSize.height);
  
          poMIEntry->SetStringField("units", poMapInfo->units);
      }
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                           HFAGetPEString()                           */
  /*                                                                      */
  /*      Some files have a ProjectionX node containing the ESRI style    */
  /*      PE_STRING.  This function allows fetching from it.              */
  /************************************************************************/
  
  char *HFAGetPEString( HFAHandle hHFA )
  
  {
      if( hHFA->nBands == 0 )
          return nullptr;
  
      // Get the HFA node.
      HFAEntry *poProX = hHFA->papoBand[0]->poNode->GetNamedChild("ProjectionX");
      if( poProX == nullptr )
          return nullptr;
  
      const char *pszType = poProX->GetStringField("projection.type.string");
      if( pszType == nullptr || !EQUAL(pszType, "PE_COORDSYS") )
          return nullptr;
  
      // Use a gross hack to scan ahead to the actual projection
      // string. We do it this way because we don't have general
      // handling for MIFObjects.
      GByte *pabyData = poProX->GetData();
      int nDataSize = poProX->GetDataSize();
  
      while( nDataSize > 10 &&
             !STARTS_WITH_CI((const char *)pabyData, "PE_COORDSYS,.") ) {
          pabyData++;
          nDataSize--;
      }
  
      if( nDataSize < 31 )
          return nullptr;
  
      // Skip ahead to the actual string.
      pabyData += 30;
      // nDataSize -= 30;
  
      return CPLStrdup((const char *)pabyData);
  }
  
  /************************************************************************/
  /*                           HFASetPEString()                           */
  /************************************************************************/
  
  CPLErr HFASetPEString( HFAHandle hHFA, const char *pszPEString )
  
  {
      if( !CPLTestBool(CPLGetConfigOption("HFA_WRITE_PE_STRING", "YES")) )
          return CE_None;
  
      // Loop over bands, setting information on each one.
      for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
      {
          // Verify we don't already have the node, since update-in-place
          // is likely to be more complicated.
          HFAEntry *poProX =
              hHFA->papoBand[iBand]->poNode->GetNamedChild("ProjectionX");
  
          // If we are setting an empty string then a missing entry is equivalent.
          if( strlen(pszPEString) == 0 && poProX == nullptr )
              continue;
  
          // Create the node.
          if( poProX == nullptr )
          {
              poProX = HFAEntry::New(hHFA, "ProjectionX", "Eprj_MapProjection842",
                                     hHFA->papoBand[iBand]->poNode);
              if( poProX->GetTypeObject() == nullptr )
                  return CE_Failure;
          }
  
          // Prepare the data area with some extra space just in case.
          GByte *pabyData =
              poProX->MakeData(static_cast<int>(700 + strlen(pszPEString)));
          if( !pabyData )
            return CE_Failure;
  
          memset(pabyData, 0, 250 + strlen(pszPEString));
  
          poProX->SetPosition();
  
          poProX->SetStringField("projection.type.string", "PE_COORDSYS");
          poProX->SetStringField("projection.MIFDictionary.string",
                                 "{0:pcstring,}Emif_String,{1:x{0:pcstring,}"
                                 "Emif_String,coordSys,}PE_COORDSYS,.");
  
          // Use a gross hack to scan ahead to the actual projection
          // string. We do it this way because we don't have general
          // handling for MIFObjects
          pabyData = poProX->GetData();
          int nDataSize = poProX->GetDataSize();
          GUInt32 iOffset = poProX->GetDataPos();
  
          while( nDataSize > 10 &&
                 !STARTS_WITH_CI((const char *) pabyData, "PE_COORDSYS,.") )
          {
              pabyData++;
              nDataSize--;
              iOffset++;
          }
  
          CPLAssert(nDataSize > static_cast<int>(strlen(pszPEString)) + 10);
  
          pabyData += 14;
          iOffset += 14;
  
          // Set the size and offset of the mifobject.
          iOffset += 8;
  
          GUInt32 nSize = static_cast<GUInt32>(strlen(pszPEString) + 9);
  
          HFAStandard(4, &nSize);
          memcpy(pabyData, &nSize, 4);
          pabyData += 4;
  
          HFAStandard(4, &iOffset);
          memcpy(pabyData, &iOffset, 4);
          pabyData += 4;
  
          // Set the size and offset of the string value.
          nSize = static_cast<GUInt32>(strlen(pszPEString) + 1);
  
          HFAStandard(4, &nSize);
          memcpy(pabyData, &nSize, 4);
          pabyData += 4;
  
          iOffset = 8;
          HFAStandard(4, &iOffset);
          memcpy(pabyData, &iOffset, 4);
          pabyData += 4;
  
          // Place the string itself.
          memcpy(pabyData, pszPEString, strlen(pszPEString) + 1);
  
          poProX->SetStringField("title.string", "PE");
      }
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                        HFAGetProParameters()                         */
  /************************************************************************/
  
  const Eprj_ProParameters *HFAGetProParameters( HFAHandle hHFA )
  
  {
      if( hHFA->nBands < 1 )
          return nullptr;
  
      // Do we already have it?
      if( hHFA->pProParameters != nullptr )
          return (Eprj_ProParameters *)hHFA->pProParameters;
  
      // Get the HFA node.
      HFAEntry *poMIEntry =
          hHFA->papoBand[0]->poNode->GetNamedChild("Projection");
      if( poMIEntry == nullptr )
          return nullptr;
  
      // Allocate the structure.
      Eprj_ProParameters *psProParms = static_cast<Eprj_ProParameters *>(
          CPLCalloc(sizeof(Eprj_ProParameters), 1));
  
      // Fetch the fields.
      const int proType = poMIEntry->GetIntField("proType");
      if( proType != EPRJ_INTERNAL && proType != EPRJ_EXTERNAL )
      {
          CPLError(CE_Failure, CPLE_AppDefined, "Wrong value for proType");
          CPLFree(psProParms);
          return nullptr;
      }
      psProParms->proType = static_cast<Eprj_ProType>(proType);
      psProParms->proNumber = poMIEntry->GetIntField("proNumber");
      psProParms->proExeName = CPLStrdup(poMIEntry->GetStringField("proExeName"));
      psProParms->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
      psProParms->proZone = poMIEntry->GetIntField("proZone");
  
      for( int i = 0; i < 15; i++ )
      {
          char szFieldName[40] = {};
  
          snprintf(szFieldName, sizeof(szFieldName), "proParams[%d]", i);
          psProParms->proParams[i] = poMIEntry->GetDoubleField(szFieldName);
      }
  
      psProParms->proSpheroid.sphereName =
          CPLStrdup(poMIEntry->GetStringField("proSpheroid.sphereName"));
      psProParms->proSpheroid.a = poMIEntry->GetDoubleField("proSpheroid.a");
      psProParms->proSpheroid.b = poMIEntry->GetDoubleField("proSpheroid.b");
      psProParms->proSpheroid.eSquared =
          poMIEntry->GetDoubleField("proSpheroid.eSquared");
      psProParms->proSpheroid.radius =
          poMIEntry->GetDoubleField("proSpheroid.radius");
  
      hHFA->pProParameters = (void *)psProParms;
  
      return psProParms;
  }
  
  /************************************************************************/
  /*                        HFASetProParameters()                         */
  /************************************************************************/
  
  CPLErr HFASetProParameters( HFAHandle hHFA, const Eprj_ProParameters *poPro )
  
  {
      // Loop over bands, setting information on each one.
      for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
      {
          // Create a new Projection if there isn't one present already.
          HFAEntry *poMIEntry =
              hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
          if( poMIEntry == nullptr )
          {
              poMIEntry = HFAEntry::New(hHFA, "Projection", "Eprj_ProParameters",
                                        hHFA->papoBand[iBand]->poNode);
          }
  
          poMIEntry->MarkDirty();
  
          // Ensure we have enough space for all the data.
          // TODO(schwehr): Explain all these constants.
          int nSize =
              static_cast<int>(34 + 15 * 8 + 8 + strlen(poPro->proName) + 1 + 32 +
                               8 + strlen(poPro->proSpheroid.sphereName) + 1);
  
          if( poPro->proExeName != nullptr )
              nSize += static_cast<int>(strlen(poPro->proExeName) + 1);
  
          GByte *pabyData = poMIEntry->MakeData(nSize);
          if( !pabyData )
              return CE_Failure;
  
          poMIEntry->SetPosition();
  
          // Initialize the whole thing to zeros for a clean start.
          memset(poMIEntry->GetData(), 0, poMIEntry->GetDataSize());
  
          // Write the various fields.
          poMIEntry->SetIntField("proType", poPro->proType);
  
          poMIEntry->SetIntField("proNumber", poPro->proNumber);
  
          poMIEntry->SetStringField("proExeName", poPro->proExeName);
          poMIEntry->SetStringField("proName", poPro->proName);
          poMIEntry->SetIntField("proZone", poPro->proZone);
          poMIEntry->SetDoubleField("proParams[0]", poPro->proParams[0]);
          poMIEntry->SetDoubleField("proParams[1]", poPro->proParams[1]);
          poMIEntry->SetDoubleField("proParams[2]", poPro->proParams[2]);
          poMIEntry->SetDoubleField("proParams[3]", poPro->proParams[3]);
          poMIEntry->SetDoubleField("proParams[4]", poPro->proParams[4]);
          poMIEntry->SetDoubleField("proParams[5]", poPro->proParams[5]);
          poMIEntry->SetDoubleField("proParams[6]", poPro->proParams[6]);
          poMIEntry->SetDoubleField("proParams[7]", poPro->proParams[7]);
          poMIEntry->SetDoubleField("proParams[8]", poPro->proParams[8]);
          poMIEntry->SetDoubleField("proParams[9]", poPro->proParams[9]);
          poMIEntry->SetDoubleField("proParams[10]", poPro->proParams[10]);
          poMIEntry->SetDoubleField("proParams[11]", poPro->proParams[11]);
          poMIEntry->SetDoubleField("proParams[12]", poPro->proParams[12]);
          poMIEntry->SetDoubleField("proParams[13]", poPro->proParams[13]);
          poMIEntry->SetDoubleField("proParams[14]", poPro->proParams[14]);
          poMIEntry->SetStringField("proSpheroid.sphereName",
                                    poPro->proSpheroid.sphereName);
          poMIEntry->SetDoubleField("proSpheroid.a", poPro->proSpheroid.a);
          poMIEntry->SetDoubleField("proSpheroid.b", poPro->proSpheroid.b);
          poMIEntry->SetDoubleField("proSpheroid.eSquared",
                                    poPro->proSpheroid.eSquared);
          poMIEntry->SetDoubleField("proSpheroid.radius",
                                    poPro->proSpheroid.radius);
      }
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                            HFAGetDatum()                             */
  /************************************************************************/
  
  const Eprj_Datum *HFAGetDatum( HFAHandle hHFA )
  
  {
      if( hHFA->nBands < 1 )
          return nullptr;
  
      // Do we already have it?
      if( hHFA->pDatum != nullptr )
          return (Eprj_Datum *)hHFA->pDatum;
  
      // Get the HFA node.
      HFAEntry *poMIEntry =
          hHFA->papoBand[0]->poNode->GetNamedChild("Projection.Datum");
      if( poMIEntry == nullptr )
          return nullptr;
  
      // Allocate the structure.
      Eprj_Datum *psDatum =
          static_cast<Eprj_Datum *>(CPLCalloc(sizeof(Eprj_Datum), 1));
  
      // Fetch the fields.
      psDatum->datumname = CPLStrdup(poMIEntry->GetStringField("datumname"));
      const int nDatumType = poMIEntry->GetIntField("type");
      if( nDatumType < 0 || nDatumType > EPRJ_DATUM_NONE )
      {
          CPLDebug("HFA", "Invalid value for datum type: %d", nDatumType);
          psDatum->type = EPRJ_DATUM_NONE;
      }
      else
          psDatum->type = static_cast<Eprj_DatumType>(nDatumType);
  
      for( int i = 0; i < 7; i++ )
      {
          char szFieldName[30] = {};
          snprintf(szFieldName, sizeof(szFieldName), "params[%d]", i);
          psDatum->params[i] = poMIEntry->GetDoubleField(szFieldName);
      }
  
      psDatum->gridname = CPLStrdup(poMIEntry->GetStringField("gridname"));
  
      hHFA->pDatum = (void *)psDatum;
  
      return psDatum;
  }
  
  /************************************************************************/
  /*                            HFASetDatum()                             */
  /************************************************************************/
  
  CPLErr HFASetDatum( HFAHandle hHFA, const Eprj_Datum *poDatum )
  
  {
      // Loop over bands, setting information on each one.
      for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
      {
          // Create a new Projection if there isn't one present already.
          HFAEntry *poProParms =
              hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
          if( poProParms == nullptr )
          {
              CPLError(CE_Failure, CPLE_AppDefined,
                       "Can't add Eprj_Datum with no Eprj_ProjParameters.");
              return CE_Failure;
          }
  
          HFAEntry *poDatumEntry = poProParms->GetNamedChild("Datum");
          if( poDatumEntry == nullptr )
          {
              poDatumEntry =
                  HFAEntry::New(hHFA, "Datum", "Eprj_Datum", poProParms);
          }
  
          poDatumEntry->MarkDirty();
  
          // Ensure we have enough space for all the data.
          // TODO(schwehr): Explain constants.
          int nSize =
              static_cast<int>(26 + strlen(poDatum->datumname) + 1 + 7 * 8);
  
          if( poDatum->gridname != nullptr )
              nSize += static_cast<int>(strlen(poDatum->gridname) + 1);
  
          GByte *pabyData = poDatumEntry->MakeData(nSize);
          if( !pabyData )
              return CE_Failure;
  
          poDatumEntry->SetPosition();
  
          // Initialize the whole thing to zeros for a clean start.
          memset(poDatumEntry->GetData(), 0, poDatumEntry->GetDataSize());
  
          // Write the various fields.
          poDatumEntry->SetStringField("datumname", poDatum->datumname);
          poDatumEntry->SetIntField("type", poDatum->type);
  
          poDatumEntry->SetDoubleField("params[0]", poDatum->params[0]);
          poDatumEntry->SetDoubleField("params[1]", poDatum->params[1]);
          poDatumEntry->SetDoubleField("params[2]", poDatum->params[2]);
          poDatumEntry->SetDoubleField("params[3]", poDatum->params[3]);
          poDatumEntry->SetDoubleField("params[4]", poDatum->params[4]);
          poDatumEntry->SetDoubleField("params[5]", poDatum->params[5]);
          poDatumEntry->SetDoubleField("params[6]", poDatum->params[6]);
  
          poDatumEntry->SetStringField("gridname", poDatum->gridname);
      }
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                             HFAGetPCT()                              */
  /*                                                                      */
  /*      Read the PCT from a band, if it has one.                        */
  /************************************************************************/
  
  CPLErr HFAGetPCT( HFAHandle hHFA, int nBand, int *pnColors,
                    double **ppadfRed, double **ppadfGreen,
                    double **ppadfBlue, double **ppadfAlpha,
                    double **ppadfBins )
  
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return CE_Failure;
  
      return hHFA->papoBand[nBand - 1]->GetPCT(pnColors, ppadfRed, ppadfGreen,
                                               ppadfBlue, ppadfAlpha, ppadfBins);
  }
  
  /************************************************************************/
  /*                             HFASetPCT()                              */
  /*                                                                      */
  /*      Set the PCT on a band.                                          */
  /************************************************************************/
  
  CPLErr HFASetPCT( HFAHandle hHFA, int nBand, int nColors,
                    double *padfRed, double *padfGreen, double *padfBlue,
                    double *padfAlpha )
  
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return CE_Failure;
  
      return hHFA->papoBand[nBand - 1]->SetPCT(nColors, padfRed, padfGreen,
                                               padfBlue, padfAlpha);
  }
  
  /************************************************************************/
  /*                          HFAGetDataRange()                           */
  /************************************************************************/
  
  CPLErr HFAGetDataRange( HFAHandle hHFA, int nBand,
                          double *pdfMin, double *pdfMax )
  
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return CE_Failure;
  
      HFAEntry *poBinInfo =
          hHFA->papoBand[nBand-1]->poNode->GetNamedChild("Statistics" );
  
      if( poBinInfo == nullptr )
          return CE_Failure;
  
      *pdfMin = poBinInfo->GetDoubleField("minimum");
      *pdfMax = poBinInfo->GetDoubleField("maximum");
  
      if( *pdfMax > *pdfMin )
          return CE_None;
      else
          return CE_Failure;
  }
  
  /************************************************************************/
  /*                            HFADumpNode()                             */
  /************************************************************************/
  
  static void HFADumpNode( HFAEntry *poEntry, int nIndent, bool bVerbose,
                           FILE *fp )
  
  {
      std::string osSpaces(nIndent * 2, ' ');
  
      fprintf(fp, "%s%s(%s) @ %u + %u @ %u\n", osSpaces.c_str(),
              poEntry->GetName(), poEntry->GetType(),
              poEntry->GetFilePos(),
              poEntry->GetDataSize(), poEntry->GetDataPos());
  
      if( bVerbose )
      {
          osSpaces += "+ ";
          poEntry->DumpFieldValues(fp, osSpaces.c_str());
          fprintf(fp, "\n");
      }
  
      if( poEntry->GetChild() != nullptr )
          HFADumpNode(poEntry->GetChild(), nIndent+1, bVerbose, fp);
  
      if( poEntry->GetNext() != nullptr )
          HFADumpNode(poEntry->GetNext(), nIndent, bVerbose, fp);
  }
  
  /************************************************************************/
  /*                            HFADumpTree()                             */
  /*                                                                      */
  /*      Dump the tree of information in a HFA file.                     */
  /************************************************************************/
  
  void HFADumpTree( HFAHandle hHFA, FILE *fpOut )
  
  {
      HFADumpNode(hHFA->poRoot, 0, true, fpOut);
  }
  
  /************************************************************************/
  /*                         HFADumpDictionary()                          */
  /*                                                                      */
  /*      Dump the dictionary (in raw, and parsed form) to the named      */
  /*      device.                                                         */
  /************************************************************************/
  
  void HFADumpDictionary( HFAHandle hHFA, FILE *fpOut )
  
  {
      fprintf(fpOut, "%s\n", hHFA->pszDictionary);
  
      hHFA->poDictionary->Dump(fpOut);
  }
  
  /************************************************************************/
  /*                            HFAStandard()                             */
  /*                                                                      */
  /*      Swap byte order on MSB systems.                                 */
  /************************************************************************/
  
  #ifdef CPL_MSB
  void HFAStandard( int nBytes, void *pData )
  
  {
      GByte *pabyData = static_cast<GByte *>(pData);
  
      for( int i = nBytes / 2 - 1; i >= 0; i-- )
      {
          GByte byTemp = pabyData[i];
          pabyData[i] = pabyData[nBytes - i - 1];
          pabyData[nBytes - i - 1] = byTemp;
      }
  }
  #endif
  
  /* ==================================================================== */
  /*      Default data dictionary.  Emitted verbatim into the imagine     */
  /*      file.                                                           */
  /* ==================================================================== */
  
  static const char * const aszDefaultDD[] = {
  "{1:lversion,1:LfreeList,1:LrootEntryPtr,1:sentryHeaderLength,1:LdictionaryPtr,}Ehfa_File,{1:Lnext,1:Lprev,1:Lparent,1:Lchild,1:Ldata,1:ldataSize,64:cname,32:ctype,1:tmodTime,}Ehfa_Entry,{16:clabel,1:LheaderPtr,}Ehfa_HeaderTag,{1:LfreeList,1:lfreeSize,}Ehfa_FreeListNode,{1:lsize,1:Lptr,}Ehfa_Data,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft of real-valued data,layerType,",
  "1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft of real-valued data,layerType,1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer_SubSample,{1:e2:raster,vector,type,1:LdictionaryPtr,}Ehfa_Layer,{1:LspaceUsedForRasterData,}ImgFormatInfo831,{1:sfileCode,1:Loffset,1:lsize,1:e2:false,true,logvalid,",
  "1:e2:no compression,ESRI GRID compression,compressionType,}Edms_VirtualBlockInfo,{1:lmin,1:lmax,}Edms_FreeIDList,{1:lnumvirtualblocks,1:lnumobjectsperblock,1:lnextobjectnum,1:e2:no compression,RLC compression,compressionType,0:poEdms_VirtualBlockInfo,blockinfo,0:poEdms_FreeIDList,freelist,1:tmodTime,}Edms_State,{0:pcstring,}Emif_String,{1:oEmif_String,fileName,2:LlayerStackValidFlagsOffset,2:LlayerStackDataOffset,1:LlayerStackCount,1:LlayerStackIndex,}ImgExternalRaster,{1:oEmif_String,algorithm,0:poEmif_String,nameList,}Eimg_RRDNamesList,{1:oEmif_String,projection,1:oEmif_String,units,}Eimg_MapInformation,",
  "{1:oEmif_String,dependent,}Eimg_DependentFile,{1:oEmif_String,ImageLayerName,}Eimg_DependentLayerName,{1:lnumrows,1:lnumcolumns,1:e13:EGDA_TYPE_U1,EGDA_TYPE_U2,EGDA_TYPE_U4,EGDA_TYPE_U8,EGDA_TYPE_S8,EGDA_TYPE_U16,EGDA_TYPE_S16,EGDA_TYPE_U32,EGDA_TYPE_S32,EGDA_TYPE_F32,EGDA_TYPE_F64,EGDA_TYPE_C64,EGDA_TYPE_C128,datatype,1:e4:EGDA_SCALAR_OBJECT,EGDA_TABLE_OBJECT,EGDA_MATRIX_OBJECT,EGDA_RASTER_OBJECT,objecttype,}Egda_BaseData,{1:*bvalueBD,}Eimg_NonInitializedValue,{1:dx,1:dy,}Eprj_Coordinate,{1:dwidth,1:dheight,}Eprj_Size,{0:pcproName,1:*oEprj_Coordinate,upperLeftCenter,",
  "1:*oEprj_Coordinate,lowerRightCenter,1:*oEprj_Size,pixelSize,0:pcunits,}Eprj_MapInfo,{0:pcdatumname,1:e3:EPRJ_DATUM_PARAMETRIC,EPRJ_DATUM_GRID,EPRJ_DATUM_REGRESSION,type,0:pdparams,0:pcgridname,}Eprj_Datum,{0:pcsphereName,1:da,1:db,1:deSquared,1:dradius,}Eprj_Spheroid,{1:e2:EPRJ_INTERNAL,EPRJ_EXTERNAL,proType,1:lproNumber,0:pcproExeName,0:pcproName,1:lproZone,0:pdproParams,1:*oEprj_Spheroid,proSpheroid,}Eprj_ProParameters,{1:dminimum,1:dmaximum,1:dmean,1:dmedian,1:dmode,1:dstddev,}Esta_Statistics,{1:lnumBins,1:e4:direct,linear,logarithmic,explicit,binFunctionType,1:dminLimit,1:dmaxLimit,1:*bbinLimits,}Edsc_BinFunction,{0:poEmif_String,LayerNames,1:*bExcludedValues,1:oEmif_String,AOIname,",
  "1:lSkipFactorX,1:lSkipFactorY,1:*oEdsc_BinFunction,BinFunction,}Eimg_StatisticsParameters830,{1:lnumrows,}Edsc_Table,{1:lnumRows,1:LcolumnDataPtr,1:e4:integer,real,complex,string,dataType,1:lmaxNumChars,}Edsc_Column,{1:lposition,0:pcname,1:e2:EMSC_FALSE,EMSC_TRUE,editable,1:e3:LEFT,CENTER,RIGHT,alignment,0:pcformat,1:e3:DEFAULT,APPLY,AUTO-APPLY,formulamode,0:pcformula,1:dcolumnwidth,0:pcunits,1:e5:NO_COLOR,RED,GREEN,BLUE,COLOR,colorflag,0:pcgreenname,0:pcbluename,}Eded_ColumnAttributes_1,{1:lversion,1:lnumobjects,1:e2:EAOI_UNION,EAOI_INTERSECTION,operation,}Eaoi_AreaOfInterest,",
  "{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,MIFDictionary,0:pCMIFObject,}Emif_MIFObject,",
  "{1:x{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,MIFDictionary,0:pCMIFObject,}Emif_MIFObject,projection,1:x{0:pcstring,}Emif_String,title,}Eprj_MapProjection842,",
  "{0:poEmif_String,titleList,}Exfr_GenericXFormHeader,{1:lorder,1:lnumdimtransform,1:lnumdimpolynomial,1:ltermcount,0:plexponentlist,1:*bpolycoefmtx,1:*bpolycoefvector,}Efga_Polynomial,",
  ".",
  nullptr
  };
  
  /************************************************************************/
  /*                            HFACreateLL()                             */
  /*                                                                      */
  /*      Low level creation of an Imagine file.  Writes out the          */
  /*      Ehfa_HeaderTag, dictionary and Ehfa_File.                       */
  /************************************************************************/
  
  HFAHandle HFACreateLL( const char *pszFilename )
  
  {
      // Create the file in the file system.
      VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
      if( fp == nullptr )
      {
          CPLError(CE_Failure, CPLE_OpenFailed, "Creation of file %s failed.",
                   pszFilename);
          return nullptr;
      }
  
      // Create the HFAInfo_t.
      HFAInfo_t *psInfo =
          static_cast<HFAInfo_t *>(CPLCalloc(sizeof(HFAInfo_t), 1));
  
      psInfo->fp = fp;
      psInfo->eAccess = HFA_Update;
      psInfo->nXSize = 0;
      psInfo->nYSize = 0;
      psInfo->nBands = 0;
      psInfo->papoBand = nullptr;
      psInfo->pMapInfo = nullptr;
      psInfo->pDatum = nullptr;
      psInfo->pProParameters = nullptr;
      psInfo->bTreeDirty = false;
      psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
      psInfo->pszPath = CPLStrdup(CPLGetPath(pszFilename));
  
      // Write out the Ehfa_HeaderTag.
      bool bRet = VSIFWriteL((void *)"EHFA_HEADER_TAG", 1, 16, fp) > 0;
  
      GInt32 nHeaderPos = 20;
      HFAStandard(4, &nHeaderPos);
      bRet &= VSIFWriteL(&nHeaderPos, 4, 1, fp) > 0;
  
      // Write the Ehfa_File node, locked in at offset 20.
      GInt32 nVersion = 1;
      GInt32 nFreeList = 0;
      GInt32 nRootEntry = 0;
      GInt16 nEntryHeaderLength = 128;
      GInt32 nDictionaryPtr = 38;
  
      psInfo->nEntryHeaderLength = nEntryHeaderLength;
      psInfo->nRootPos = 0;
      psInfo->nDictionaryPos = nDictionaryPtr;
      psInfo->nVersion = nVersion;
  
      HFAStandard(4, &nVersion);
      HFAStandard(4, &nFreeList);
      HFAStandard(4, &nRootEntry);
      HFAStandard(2, &nEntryHeaderLength);
      HFAStandard(4, &nDictionaryPtr);
  
      bRet &= VSIFWriteL(&nVersion, 4, 1, fp) > 0;
      bRet &= VSIFWriteL(&nFreeList, 4, 1, fp) > 0;
      bRet &= VSIFWriteL(&nRootEntry, 4, 1, fp) > 0;
      bRet &= VSIFWriteL(&nEntryHeaderLength, 2, 1, fp) > 0;
      bRet &= VSIFWriteL(&nDictionaryPtr, 4, 1, fp) > 0;
  
      // Write the dictionary, locked in at location 38.  Note that
      // we jump through a bunch of hoops to operate on the
      // dictionary in chunks because some compiles (such as VC++)
      // don't allow particularly large static strings.
      int nDictLen = 0;
  
      for( int iChunk = 0; aszDefaultDD[iChunk] != nullptr; iChunk++ )
          nDictLen += static_cast<int>(strlen(aszDefaultDD[iChunk]));
  
      psInfo->pszDictionary = static_cast<char *>(CPLMalloc(nDictLen + 1));
      psInfo->pszDictionary[0] = '\0';
  
      for( int iChunk = 0; aszDefaultDD[iChunk] != nullptr; iChunk++ )
          strcat(psInfo->pszDictionary, aszDefaultDD[iChunk]);
  
      bRet &= VSIFWriteL((void *)psInfo->pszDictionary,
                         strlen(psInfo->pszDictionary) + 1, 1, fp) > 0;
      if( !bRet )
      {
          CPL_IGNORE_RET_VAL(HFAClose(psInfo));
          return nullptr;
      }
  
      psInfo->poDictionary = new HFADictionary(psInfo->pszDictionary);
  
      psInfo->nEndOfFile = static_cast<GUInt32>(VSIFTellL(fp));
  
      // Create a root entry.
      psInfo->poRoot = new HFAEntry(psInfo, "root", "root", nullptr);
  
      // If an .ige or .rrd file exists with the same base name,
      // delete them.  (#1784)
      CPLString osExtension = CPLGetExtension(pszFilename);
      if( !EQUAL(osExtension, "rrd") && !EQUAL(osExtension, "aux") )
      {
          CPLString osPath = CPLGetPath(pszFilename);
          CPLString osBasename = CPLGetBasename(pszFilename);
          VSIStatBufL sStatBuf;
          CPLString osSupFile = CPLFormCIFilename(osPath, osBasename, "rrd");
  
          if( VSIStatL(osSupFile, &sStatBuf) == 0 )
              VSIUnlink(osSupFile);
  
          osSupFile = CPLFormCIFilename(osPath, osBasename, "ige");
  
          if( VSIStatL(osSupFile, &sStatBuf) == 0 )
              VSIUnlink(osSupFile);
      }
  
      return psInfo;
  }
  
  /************************************************************************/
  /*                          HFAAllocateSpace()                          */
  /*                                                                      */
  /*      Return an area in the file to the caller to write the           */
  /*      requested number of bytes.  Currently this is always at the     */
  /*      end of the file, but eventually we might actually keep track    */
  /*      of free space.  The HFAInfo_t's concept of file size is         */
  /*      updated, even if nothing ever gets written to this region.      */
  /*                                                                      */
  /*      Returns the offset to the requested space, or zero one          */
  /*      failure.                                                        */
  /************************************************************************/
  
  GUInt32 HFAAllocateSpace( HFAInfo_t *psInfo, GUInt32 nBytes )
  
  {
      // TODO(schwehr): Check if this will wrap over 2GB limit.
  
      psInfo->nEndOfFile += nBytes;
      return psInfo->nEndOfFile - nBytes;
  }
  
  /************************************************************************/
  /*                              HFAFlush()                              */
  /*                                                                      */
  /*      Write out any dirty tree information to disk, putting the       */
  /*      disk file in a consistent state.                                */
  /************************************************************************/
  
  CPLErr HFAFlush( HFAHandle hHFA )
  
  {
      if( !hHFA->bTreeDirty && !hHFA->poDictionary->bDictionaryTextDirty )
          return CE_None;
  
      CPLAssert(hHFA->poRoot != nullptr);
  
      // Flush HFAEntry tree to disk.
      if( hHFA->bTreeDirty )
      {
          const CPLErr eErr = hHFA->poRoot->FlushToDisk();
          if( eErr != CE_None )
              return eErr;
  
          hHFA->bTreeDirty = false;
      }
  
      // Flush Dictionary to disk.
      GUInt32 nNewDictionaryPos = hHFA->nDictionaryPos;
      bool bRet = true;
      if( hHFA->poDictionary->bDictionaryTextDirty )
      {
          bRet &= VSIFSeekL(hHFA->fp, 0, SEEK_END) >= 0;
          nNewDictionaryPos = static_cast<GUInt32>(VSIFTellL(hHFA->fp));
          bRet &=
              VSIFWriteL(hHFA->poDictionary->osDictionaryText.c_str(),
                         strlen(hHFA->poDictionary->osDictionaryText.c_str()) + 1,
                         1, hHFA->fp) > 0;
          hHFA->poDictionary->bDictionaryTextDirty = false;
      }
  
      // Do we need to update the Ehfa_File pointer to the root node?
      if( hHFA->nRootPos != hHFA->poRoot->GetFilePos() ||
          nNewDictionaryPos != hHFA->nDictionaryPos )
      {
          GUInt32 nHeaderPos = 0;
  
          bRet &= VSIFSeekL(hHFA->fp, 16, SEEK_SET) >= 0;
          bRet &= VSIFReadL(&nHeaderPos, sizeof(GInt32), 1, hHFA->fp) > 0;
          HFAStandard(4, &nHeaderPos);
  
          GUInt32 nOffset = hHFA->poRoot->GetFilePos();
          hHFA->nRootPos = nOffset;
          HFAStandard(4, &nOffset);
          bRet &= VSIFSeekL(hHFA->fp, nHeaderPos + 8, SEEK_SET) >= 0;
          bRet &= VSIFWriteL(&nOffset, 4, 1, hHFA->fp) > 0;
  
          nOffset = nNewDictionaryPos;
          hHFA->nDictionaryPos = nNewDictionaryPos;
          HFAStandard(4, &nOffset);
          bRet &= VSIFSeekL(hHFA->fp, nHeaderPos + 14, SEEK_SET) >= 0;
          bRet &= VSIFWriteL(&nOffset, 4, 1, hHFA->fp) > 0;
      }
  
      return bRet ? CE_None : CE_Failure;
  }
  
  /************************************************************************/
  /*                           HFACreateLayer()                           */
  /*                                                                      */
  /*      Create a layer object, and corresponding RasterDMS.             */
  /*      Suitable for use with primary layers, and overviews.            */
  /************************************************************************/
  
  int
  HFACreateLayer( HFAHandle psInfo, HFAEntry *poParent,
                  const char *pszLayerName,
                  int bOverview, int nBlockSize,
                  int bCreateCompressed, int bCreateLargeRaster,
                  int bDependentLayer,
                  int nXSize, int nYSize, EPTType eDataType,
                  char ** /* papszOptions */,
                  // These are only related to external (large) files.
                  GIntBig nStackValidFlagsOffset,
                  GIntBig nStackDataOffset,
                  int nStackCount, int nStackIndex )
  
  {
      const char *pszLayerType =
          bOverview ? "Eimg_Layer_SubSample" : "Eimg_Layer";
  
      if( nBlockSize <= 0 )
      {
          CPLError(CE_Failure, CPLE_IllegalArg,
                   "HFACreateLayer: nBlockXSize < 0");
          return FALSE;
      }
  
      // Work out some details about the tiling scheme.
      const int nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
      const int nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
      const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
      const int nBytesPerBlock =
          (nBlockSize * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
  
      // Create the Eimg_Layer for the band.
      HFAEntry *poEimg_Layer =
          HFAEntry::New(psInfo, pszLayerName, pszLayerType, poParent);
  
      poEimg_Layer->SetIntField("width", nXSize);
      poEimg_Layer->SetIntField("height", nYSize);
      poEimg_Layer->SetStringField("layerType", "athematic");
      poEimg_Layer->SetIntField("pixelType", eDataType);
      poEimg_Layer->SetIntField("blockWidth", nBlockSize);
      poEimg_Layer->SetIntField("blockHeight", nBlockSize);
  
      // Create the RasterDMS (block list).  This is a complex type
      // with pointers, and variable size.  We set the superstructure
      // ourselves rather than trying to have the HFA type management
      // system do it for us (since this would be hard to implement).
      if( !bCreateLargeRaster && !bDependentLayer )
      {
          HFAEntry *poEdms_State =
              HFAEntry::New(psInfo, "RasterDMS", "Edms_State", poEimg_Layer);
  
          // TODO(schwehr): Explain constants.
          const int nDmsSize = 14 * nBlocks + 38;
          GByte *pabyData = poEdms_State->MakeData(nDmsSize);
  
          // Set some simple values.
          poEdms_State->SetIntField("numvirtualblocks", nBlocks);
          poEdms_State->SetIntField("numobjectsperblock",
                                    nBlockSize * nBlockSize);
          poEdms_State->SetIntField("nextobjectnum",
                                    nBlockSize * nBlockSize * nBlocks);
  
          // Is file compressed or not?
          if( bCreateCompressed )
          {
              poEdms_State->SetStringField("compressionType", "RLC compression");
          }
          else
          {
              poEdms_State->SetStringField("compressionType", "no compression");
          }
  
          // We need to hardcode file offset into the data, so locate it now.
          poEdms_State->SetPosition();
  
          // Set block info headers.
  
          // Blockinfo count.
          GUInt32 nValue = nBlocks;
          HFAStandard(4, &nValue);
          memcpy(pabyData + 14, &nValue, 4);
  
          // Blockinfo position.
          nValue = poEdms_State->GetDataPos() + 22;
          HFAStandard(4, &nValue);
          memcpy(pabyData + 18, &nValue, 4);
  
          // Set each blockinfo.
          for( int iBlock = 0; iBlock < nBlocks; iBlock++ )
          {
              int nOffset = 22 + 14 * iBlock;
  
              // fileCode.
              GInt16 nValue16 = 0;
              HFAStandard(2, &nValue16);
              memcpy(pabyData + nOffset, &nValue16, 2);
  
              // Offset.
              if( bCreateCompressed )
              {
                  // Flag it with zero offset. Allocate space when we compress it.
                  nValue = 0;
              }
              else
              {
                  nValue = HFAAllocateSpace(psInfo, nBytesPerBlock);
              }
              HFAStandard(4, &nValue);
              memcpy(pabyData + nOffset + 2, &nValue, 4);
  
              // Size.
              if( bCreateCompressed )
              {
                  // Flag with zero size. Don't know until we compress it.
                  nValue = 0;
              }
              else
              {
                  nValue = nBytesPerBlock;
              }
              HFAStandard(4, &nValue);
              memcpy(pabyData + nOffset + 6, &nValue, 4);
  
              // logValid (false).
              nValue16 = 0;
              HFAStandard(2, &nValue16);
              memcpy(pabyData + nOffset + 10, &nValue16, 2);
  
              // compressionType.
              if( bCreateCompressed )
                  nValue16 = 1;
              else
                  nValue16 = 0;
  
              HFAStandard(2, &nValue16);
              memcpy(pabyData + nOffset + 12, &nValue16, 2);
          }
      }
  
      // Create ExternalRasterDMS object.
      else if( bCreateLargeRaster )
      {
          HFAEntry *poEdms_State = HFAEntry::New(
              psInfo, "ExternalRasterDMS", "ImgExternalRaster", poEimg_Layer);
          poEdms_State->MakeData(
              static_cast<int>(8 + strlen(psInfo->pszIGEFilename) + 1 + 6 * 4));
  
          poEdms_State->SetStringField("fileName.string", psInfo->pszIGEFilename);
  
          poEdms_State->SetIntField(
              "layerStackValidFlagsOffset[0]",
              static_cast<int>(nStackValidFlagsOffset & 0xFFFFFFFF));
          poEdms_State->SetIntField(
              "layerStackValidFlagsOffset[1]",
              static_cast<int>(nStackValidFlagsOffset >> 32));
  
          poEdms_State->SetIntField(
              "layerStackDataOffset[0]",
              static_cast<int>(nStackDataOffset & 0xFFFFFFFF));
          poEdms_State->SetIntField("layerStackDataOffset[1]",
                                    static_cast<int>(nStackDataOffset >> 32));
          poEdms_State->SetIntField("layerStackCount", nStackCount);
          poEdms_State->SetIntField("layerStackIndex", nStackIndex);
      }
      // Dependent...
      else if( bDependentLayer )
      {
          HFAEntry *poDepLayerName =
              HFAEntry::New(psInfo, "DependentLayerName",
                            "Eimg_DependentLayerName", poEimg_Layer);
          poDepLayerName->MakeData(
              static_cast<int>(8 + strlen(pszLayerName) + 2));
  
          poDepLayerName->SetStringField("ImageLayerName.string", pszLayerName);
      }
  
      // Create the Ehfa_Layer.
      char chBandType = '\0';
  
      if(eDataType == EPT_u1)
          chBandType = '1';
      else if(eDataType == EPT_u2)
          chBandType = '2';
      else if(eDataType == EPT_u4)
          chBandType = '4';
      else if(eDataType == EPT_u8)
          chBandType = 'c';
      else if(eDataType == EPT_s8)
          chBandType = 'C';
      else if(eDataType == EPT_u16)
          chBandType = 's';
      else if(eDataType == EPT_s16)
          chBandType = 'S';
      else if(eDataType == EPT_u32)
          // For some reason erdas imagine expects an L for unsigned 32 bit ints
          // otherwise it gives strange "out of memory errors".
          chBandType = 'L';
      else if(eDataType == EPT_s32)
          chBandType = 'L';
      else if(eDataType == EPT_f32)
          chBandType = 'f';
      else if(eDataType == EPT_f64)
          chBandType = 'd';
      else if(eDataType == EPT_c64)
          chBandType = 'm';
      else if(eDataType == EPT_c128)
          chBandType = 'M';
      else
      {
          CPLAssert(false);
          chBandType = 'c';
      }
  
      // The first value in the entry below gives the number of pixels
      // within a block.
      char szLDict[128] = {};
      snprintf(szLDict, sizeof(szLDict), "{%d:%cdata,}RasterDMS,.",
               nBlockSize * nBlockSize, chBandType);
  
      HFAEntry *poEhfa_Layer =
          HFAEntry::New(psInfo, "Ehfa_Layer", "Ehfa_Layer", poEimg_Layer);
      poEhfa_Layer->MakeData();
      poEhfa_Layer->SetPosition();
      const GUInt32 nLDict =
          HFAAllocateSpace(psInfo, static_cast<GUInt32>(strlen(szLDict) + 1));
  
      poEhfa_Layer->SetStringField("type", "raster");
      poEhfa_Layer->SetIntField("dictionaryPtr", nLDict);
  
      bool bRet = VSIFSeekL(psInfo->fp, nLDict, SEEK_SET) >= 0;
      bRet &= VSIFWriteL((void *)szLDict, strlen(szLDict) + 1, 1, psInfo->fp) > 0;
  
      return bRet;
  }
  
  /************************************************************************/
  /*                             HFACreate()                              */
  /************************************************************************/
  
  HFAHandle HFACreate( const char *pszFilename,
                       int nXSize, int nYSize, int nBands,
                       EPTType eDataType, char **papszOptions )
  
  {
      int nBlockSize = 64;
      const char *pszValue = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
  
      if( pszValue != nullptr )
      {
          nBlockSize = atoi(pszValue);
          // Check for sane values.
          if( nBlockSize == 0 ||
              ((( nBlockSize < 32 ) || (nBlockSize > 2048)) &&
              !CPLTestBool(CPLGetConfigOption("FORCE_BLOCKSIZE", "NO"))) )
          {
              if( nBlockSize != 0 )
                  CPLError(CE_Warning, CPLE_AppDefined,
                           "Forcing BLOCKSIZE to %d", 64);
              nBlockSize = 64;
          }
      }
      bool bCreateLargeRaster = CPLFetchBool(papszOptions, "USE_SPILL", false);
      bool bCreateCompressed = CPLFetchBool(papszOptions, "COMPRESS", false) ||
                               CPLFetchBool(papszOptions, "COMPRESSED", false);
      const bool bCreateAux = CPLFetchBool(papszOptions, "AUX", false);
  
      char *pszFullFilename = nullptr;
      char *pszRawFilename = nullptr;
  
      // Work out some details about the tiling scheme.
      const int nBlocksPerRow = DIV_ROUND_UP(nXSize, nBlockSize);
      const int nBlocksPerColumn = DIV_ROUND_UP(nYSize, nBlockSize);
      if( nBlocksPerRow > INT_MAX / nBlocksPerColumn )
      {
          CPLError(CE_Failure, CPLE_NotSupported, "Too many blocks");
          return nullptr;
      }
      const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
      const GInt64 nBytesPerBlock64 =
          (static_cast<GInt64>(nBlockSize) * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
      if( nBytesPerBlock64 > INT_MAX )
      {
          CPLError(CE_Failure, CPLE_NotSupported, "Too large block");
          return nullptr;
      }
      const int nBytesPerBlock = static_cast<int>(nBytesPerBlock64);
  
      // Create the low level structure.
      HFAHandle psInfo = HFACreateLL(pszFilename);
      if( psInfo == nullptr )
          return nullptr;
  
      // Create the DependentFile node if requested.
      const char *pszDependentFile =
          CSLFetchNameValue(papszOptions, "DEPENDENT_FILE");
  
      if( pszDependentFile != nullptr )
      {
          HFAEntry *poDF = HFAEntry::New(psInfo, "DependentFile",
                                         "Eimg_DependentFile", psInfo->poRoot);
  
          poDF->MakeData(static_cast<int>(strlen(pszDependentFile) + 50));
          poDF->SetPosition();
          poDF->SetStringField("dependent.string", pszDependentFile);
      }
  
      CPLDebug("HFACreate", "Blocks per row %d, blocks per column %d, "
               "total number of blocks %d, bytes per block %d.",
               nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock);
  
      // Check whether we should create external large file with
      // image.  We create a spill file if the amount of imagery is
      // close to 2GB.  We don't check the amount of auxiliary
      // information, so in theory if there were an awful lot of
      // non-imagery data our approximate size could be smaller than
      // the file will actually we be.  We leave room for 10MB of
      // auxiliary data.
      // We can also force spill file creation using option
      // SPILL_FILE=YES.
      const double dfApproxSize =
          static_cast<double>(nBytesPerBlock) * static_cast<double>(nBlocks) *
          static_cast<double>(nBands) + 10000000.0;
  
      if( dfApproxSize > 2147483648.0 && !bCreateAux )
          bCreateLargeRaster = true;
  
      // Erdas Imagine creates this entry even if an external spill file is used.
      if( !bCreateAux )
      {
          HFAEntry *poImgFormat = HFAEntry::New(
              psInfo, "IMGFormatInfo", "ImgFormatInfo831", psInfo->poRoot);
          poImgFormat->MakeData();
          if( bCreateLargeRaster )
          {
              poImgFormat->SetIntField("spaceUsedForRasterData", 0);
              // Can't be compressed if we are creating a spillfile.
              bCreateCompressed = false;
          }
          else
          {
              poImgFormat->SetIntField("spaceUsedForRasterData",
                                       nBytesPerBlock * nBlocks * nBands);
          }
      }
  
      // Create external file and write its header.
      GIntBig nValidFlagsOffset = 0;
      GIntBig nDataOffset = 0;
  
      if( bCreateLargeRaster )
      {
          if( !HFACreateSpillStack(psInfo, nXSize, nYSize, nBands,
                                   nBlockSize, eDataType,
                                   &nValidFlagsOffset, &nDataOffset) )
          {
              CPLFree(pszRawFilename);
              CPLFree(pszFullFilename);
              return nullptr;
          }
      }
  
      // Create each band (layer).
      for( int iBand = 0; iBand < nBands; iBand++ )
      {
          char szName[128] = {};
  
          snprintf(szName, sizeof(szName), "Layer_%d", iBand + 1);
  
          if( !HFACreateLayer(psInfo, psInfo->poRoot, szName, FALSE, nBlockSize,
                              bCreateCompressed, bCreateLargeRaster, bCreateAux,
                              nXSize, nYSize, eDataType, papszOptions,
                              nValidFlagsOffset, nDataOffset,
                              nBands, iBand) )
          {
              CPL_IGNORE_RET_VAL(HFAClose(psInfo));
              return nullptr;
          }
      }
  
      // Initialize the band information.
      HFAParseBandInfo(psInfo);
  
      return psInfo;
  }
  
  /************************************************************************/
  /*                         HFACreateOverview()                          */
  /*                                                                      */
  /*      Create an overview layer object for a band.                     */
  /************************************************************************/
  
  int HFACreateOverview( HFAHandle hHFA, int nBand, int nOverviewLevel,
                         const char *pszResampling )
  
  {
      if( nBand < 1 || nBand > hHFA->nBands )
          return -1;
  
      HFABand *poBand = hHFA->papoBand[nBand - 1];
      return poBand->CreateOverview(nOverviewLevel, pszResampling);
  }
  
  /************************************************************************/
  /*                           HFAGetMetadata()                           */
  /*                                                                      */
  /*      Read metadata structured in a table called GDAL_MetaData.       */
  /************************************************************************/
  
  char **HFAGetMetadata( HFAHandle hHFA, int nBand )
  
  {
      HFAEntry *poTable = nullptr;
  
      if( nBand > 0 && nBand <= hHFA->nBands )
          poTable = hHFA->papoBand[nBand - 1]->poNode->GetChild();
      else if( nBand == 0 )
          poTable = hHFA->poRoot->GetChild();
      else
          return nullptr;
  
      for( ;
           poTable != nullptr && !EQUAL(poTable->GetName(), "GDAL_MetaData");
           poTable = poTable->GetNext() ) {}
  
      if( poTable == nullptr || !EQUAL(poTable->GetType(), "Edsc_Table") )
          return nullptr;
  
      if( poTable->GetIntField("numRows") != 1 )
      {
          CPLDebug("HFADataset", "GDAL_MetaData.numRows = %d, expected 1!",
                   poTable->GetIntField("numRows"));
          return nullptr;
      }
  
      // Loop over each column.  Each column will be one metadata
      // entry, with the title being the key, and the row value being
      // the value.  There is only ever one row in GDAL_MetaData tables.
      char **papszMD = nullptr;
  
      for( HFAEntry *poColumn = poTable->GetChild();
           poColumn != nullptr;
           poColumn = poColumn->GetNext() )
      {
          // Skip the #Bin_Function# entry.
          if( STARTS_WITH_CI(poColumn->GetName(), "#") )
              continue;
  
          const char *pszValue = poColumn->GetStringField("dataType");
          if( pszValue == nullptr || !EQUAL(pszValue, "string") )
              continue;
  
          const int columnDataPtr = poColumn->GetIntField("columnDataPtr");
          if( columnDataPtr <= 0 )
              continue;
  
          // Read up to nMaxNumChars bytes from the indicated location.
          // allocate required space temporarily nMaxNumChars should have been
          // set by GDAL originally so we should trust it, but who knows.
          const int nMaxNumChars = poColumn->GetIntField("maxNumChars");
  
          if( nMaxNumChars <= 0 )
          {
              papszMD = CSLSetNameValue(papszMD, poColumn->GetName(), "");
          }
          else
          {
              char *pszMDValue =
                  static_cast<char *>(VSI_MALLOC_VERBOSE(nMaxNumChars));
              if( pszMDValue == nullptr )
              {
                  continue;
              }
  
              if( VSIFSeekL(hHFA->fp, columnDataPtr, SEEK_SET) != 0 )
              {
                  CPLFree(pszMDValue);
                  continue;
              }
  
              const int nMDBytes = static_cast<int>(
                  VSIFReadL(pszMDValue, 1, nMaxNumChars, hHFA->fp));
              if( nMDBytes == 0 )
              {
                  CPLFree(pszMDValue);
                  continue;
              }
  
              pszMDValue[nMaxNumChars - 1] = '\0';
  
              papszMD = CSLSetNameValue(papszMD, poColumn->GetName(), pszMDValue);
              CPLFree(pszMDValue);
          }
      }
  
      return papszMD;
  }
  
  /************************************************************************/
  /*                         HFASetGDALMetadata()                         */
  /*                                                                      */
  /*      This function is used to set metadata in a table called         */
  /*      GDAL_MetaData.  It is called by HFASetMetadata() for all        */
  /*      metadata items that aren't some specific supported              */
  /*      information (like histogram or stats info).                     */
  /************************************************************************/
  
  static CPLErr
  HFASetGDALMetadata( HFAHandle hHFA, int nBand, char **papszMD )
  
  {
      if( papszMD == nullptr )
          return CE_None;
  
      HFAEntry *poNode = nullptr;
  
      if( nBand > 0 && nBand <= hHFA->nBands )
          poNode = hHFA->papoBand[nBand - 1]->poNode;
      else if( nBand == 0 )
          poNode = hHFA->poRoot;
      else
          return CE_Failure;
  
      // Create the Descriptor table.
      // Check we have no table with this name already.
      HFAEntry *poEdsc_Table = poNode->GetNamedChild("GDAL_MetaData");
  
      if( poEdsc_Table == nullptr || !EQUAL(poEdsc_Table->GetType(), "Edsc_Table") )
          poEdsc_Table =
              HFAEntry::New(hHFA, "GDAL_MetaData", "Edsc_Table", poNode);
  
      poEdsc_Table->SetIntField("numrows", 1);
  
      // Create the Binning function node.  Do we really need this though?
      // Check it doesn't exist already.
      HFAEntry *poEdsc_BinFunction =
          poEdsc_Table->GetNamedChild("#Bin_Function#");
  
      if( poEdsc_BinFunction == nullptr ||
          !EQUAL(poEdsc_BinFunction->GetType(), "Edsc_BinFunction") )
          poEdsc_BinFunction = HFAEntry::New(hHFA, "#Bin_Function#",
                                             "Edsc_BinFunction", poEdsc_Table);
  
      // Because of the BaseData we have to hardcode the size.
      poEdsc_BinFunction->MakeData( 30 );
  
      poEdsc_BinFunction->SetIntField("numBins", 1);
      poEdsc_BinFunction->SetStringField("binFunction", "direct");
      poEdsc_BinFunction->SetDoubleField("minLimit", 0.0);
      poEdsc_BinFunction->SetDoubleField("maxLimit", 0.0);
  
      // Process each metadata item as a separate column.
      bool bRet = true;
      for( int iColumn = 0; papszMD[iColumn] != nullptr; iColumn++ )
      {
          char *pszKey = nullptr;
          const char *pszValue = CPLParseNameValue(papszMD[iColumn], &pszKey);
          if( pszValue == nullptr )
              continue;
  
          // Create the Edsc_Column.
          // Check it doesn't exist already.
          HFAEntry *poEdsc_Column = poEdsc_Table->GetNamedChild(pszKey);
  
          if( poEdsc_Column == nullptr ||
              !EQUAL(poEdsc_Column->GetType(), "Edsc_Column") )
              poEdsc_Column =
                  HFAEntry::New(hHFA, pszKey, "Edsc_Column", poEdsc_Table);
  
          poEdsc_Column->SetIntField("numRows", 1);
          poEdsc_Column->SetStringField("dataType", "string");
          poEdsc_Column->SetIntField("maxNumChars",
                                     static_cast<GUInt32>(strlen(pszValue) + 1));
  
          // Write the data out.
          const int nOffset =
              HFAAllocateSpace(hHFA, static_cast<GUInt32>(strlen(pszValue) + 1));
  
          poEdsc_Column->SetIntField("columnDataPtr", nOffset);
  
          bRet &= VSIFSeekL(hHFA->fp, nOffset, SEEK_SET) >= 0;
          bRet &=
              VSIFWriteL((void *)pszValue, strlen(pszValue) + 1, 1, hHFA->fp) > 0;
  
          CPLFree(pszKey);
      }
  
      return bRet ? CE_None : CE_Failure;
  }
  
  /************************************************************************/
  /*                           HFASetMetadata()                           */
  /************************************************************************/
  
  CPLErr HFASetMetadata( HFAHandle hHFA, int nBand, char **papszMD )
  
  {
      char **papszGDALMD = nullptr;
  
      if( CSLCount(papszMD) == 0 )
          return CE_None;
  
      HFAEntry *poNode = nullptr;
  
      if( nBand > 0 && nBand <= hHFA->nBands )
          poNode = hHFA->papoBand[nBand - 1]->poNode;
      else if( nBand == 0 )
          poNode = hHFA->poRoot;
      else
          return CE_Failure;
  #ifdef DEBUG
      // To please Clang Static Analyzer (CSA).
      if( poNode == nullptr )
      {
          CPLAssert(false);
          return CE_Failure;
      }
  #endif
  
      // Check if the Metadata is an "known" entity which should be
      // stored in a better place.
      char *pszBinValues = nullptr;
      bool bCreatedHistogramParameters = false;
      bool bCreatedStatistics = false;
      const char *const *pszAuxMetaData = GetHFAAuxMetaDataList();
      // Check each metadata item.
      for( int iColumn = 0; papszMD[iColumn] != nullptr; iColumn++ )
      {
          char *pszKey = nullptr;
          const char *pszValue = CPLParseNameValue(papszMD[iColumn], &pszKey);
          if( pszValue == nullptr )
              continue;
  
          // Know look if its known.
          int i = 0;  // Used after for.
          for( ; pszAuxMetaData[i] != nullptr; i += 4 )
          {
              if( EQUALN(pszAuxMetaData[i + 2], pszKey, strlen(pszKey)) )
                  break;
          }
          if( pszAuxMetaData[i] != nullptr )
          {
              // Found one, get the right entry.
              HFAEntry *poEntry = nullptr;
  
              if( strlen(pszAuxMetaData[i]) > 0 )
                  poEntry = poNode->GetNamedChild(pszAuxMetaData[i]);
              else
                  poEntry = poNode;
  
              if( poEntry == nullptr && strlen(pszAuxMetaData[i + 3]) > 0 )
              {
                  // Child does not yet exist --> create it,
                  poEntry = HFAEntry::New(hHFA, pszAuxMetaData[i],
                                          pszAuxMetaData[i + 3], poNode);
  
                  if( STARTS_WITH_CI(pszAuxMetaData[i], "Statistics") )
                      bCreatedStatistics = true;
  
                  if( STARTS_WITH_CI(pszAuxMetaData[i], "HistogramParameters") )
                  {
                      // A bit nasty.  Need to set the string field for the object
                      // first because the SetStringField sets the count for the
                      // object BinFunction to the length of the string.
                      poEntry->MakeData(70);
                      poEntry->SetStringField("BinFunction.binFunctionType",
                                              "direct");
  
                      bCreatedHistogramParameters = true;
                  }
              }
              if( poEntry == nullptr )
              {
                  CPLFree(pszKey);
                  continue;
              }
  
              const char *pszFieldName = pszAuxMetaData[i + 1] + 1;
              switch( pszAuxMetaData[i + 1][0] )
              {
              case 'd':
              {
                  double dfValue = CPLAtof(pszValue);
                  poEntry->SetDoubleField(pszFieldName, dfValue);
              }
              break;
              case 'i':
              case 'l':
              {
                  int nValue = atoi(pszValue);
                  poEntry->SetIntField(pszFieldName, nValue);
              }
              break;
              case 's':
              case 'e':
              {
                  poEntry->SetStringField(pszFieldName, pszValue);
              }
              break;
              default:
                  CPLAssert(false);
              }
          }
          else if( STARTS_WITH_CI(pszKey, "STATISTICS_HISTOBINVALUES") )
          {
              CPLFree(pszBinValues);
              pszBinValues = CPLStrdup(pszValue);
          }
          else
          {
              papszGDALMD = CSLAddString(papszGDALMD, papszMD[iColumn]);
          }
  
          CPLFree(pszKey);
      }
  
      // Special case to write out the histogram.
      bool bRet = true;
      if( pszBinValues != nullptr )
      {
          HFAEntry *poEntry = poNode->GetNamedChild("HistogramParameters");
          if( poEntry != nullptr && bCreatedHistogramParameters )
          {
              // If this node exists we have added Histogram data -- complete with
              // some defaults.
              poEntry->SetIntField("SkipFactorX", 1);
              poEntry->SetIntField("SkipFactorY", 1);
  
              const int nNumBins = poEntry->GetIntField("BinFunction.numBins");
              const double dMinLimit =
                  poEntry->GetDoubleField("BinFunction.minLimit");
              const double dMaxLimit =
                  poEntry->GetDoubleField("BinFunction.maxLimit");
  
              // Fill the descriptor table - check it isn't there already.
              poEntry = poNode->GetNamedChild("Descriptor_Table");
              if( poEntry == nullptr || !EQUAL(poEntry->GetType(), "Edsc_Table") )
                  poEntry =
                      HFAEntry::New(hHFA, "Descriptor_Table",
                                    "Edsc_Table", poNode);
  
              poEntry->SetIntField("numRows", nNumBins);
  
              // Bin function.
              HFAEntry *poBinFunc = poEntry->GetNamedChild("#Bin_Function#");
              if( poBinFunc == nullptr ||
                  !EQUAL(poBinFunc->GetType(), "Edsc_BinFunction") )
                  poBinFunc =
                      HFAEntry::New(hHFA, "#Bin_Function#",
                                    "Edsc_BinFunction", poEntry);
  
              poBinFunc->MakeData(30);
              poBinFunc->SetIntField("numBins", nNumBins);
              poBinFunc->SetDoubleField("minLimit", dMinLimit);
              poBinFunc->SetDoubleField("maxLimit", dMaxLimit);
              // Direct for thematic layers, linear otherwise.
              if( STARTS_WITH_CI(poNode->GetStringField("layerType"),
                                 "thematic") )
                  poBinFunc->SetStringField("binFunctionType", "direct");
              else
                  poBinFunc->SetStringField("binFunctionType", "linear");
  
              // We need a child named histogram.
              HFAEntry *poHisto = poEntry->GetNamedChild("Histogram");
              if( poHisto == nullptr || !EQUAL(poHisto->GetType(), "Edsc_Column") )
                  poHisto =
                      HFAEntry::New(hHFA, "Histogram", "Edsc_Column", poEntry);
  
              poHisto->SetIntField("numRows", nNumBins);
              // Allocate space for the bin values.
              GUInt32 nOffset = HFAAllocateSpace(hHFA, nNumBins * 8);
              poHisto->SetIntField("columnDataPtr", nOffset);
              poHisto->SetStringField("dataType", "real");
              poHisto->SetIntField("maxNumChars", 0);
              // Write out histogram data.
              char *pszWork = pszBinValues;
              for( int nBin = 0; nBin < nNumBins; ++nBin )
              {
                  char *pszEnd = strchr(pszWork, '|');
                  if( pszEnd != nullptr )
                  {
                      *pszEnd = 0;
                      bRet &=
                          VSIFSeekL(hHFA->fp, nOffset + 8 * nBin, SEEK_SET) >= 0;
                      double nValue = CPLAtof(pszWork);
                      HFAStandard(8, &nValue);
  
                      bRet &= VSIFWriteL((void *)&nValue, 8, 1, hHFA->fp) > 0;
                      pszWork = pszEnd + 1;
                  }
              }
          }
          else if( poEntry != nullptr )
          {
              // In this case, there are HistogramParameters present, but we did
              // not create them. However, we might be modifying them, in the case
              // where the data has changed and the histogram counts need to be
              // updated. It could be worse than that, but that is all we are
              // going to cope with for now.  We are assuming that we did not
              // change any of the other stuff, like skip factors and so
              // forth. The main need for this case is for programs (such as
              // Imagine itself) which will happily modify the pixel values
              // without re-calculating the histogram counts.
              int nNumBins = poEntry->GetIntField("BinFunction.numBins");
              HFAEntry *poEntryDescrTbl =
                  poNode->GetNamedChild("Descriptor_Table");
              HFAEntry *poHisto = nullptr;
              if( poEntryDescrTbl != nullptr)
              {
                  poHisto = poEntryDescrTbl->GetNamedChild("Histogram");
              }
              if( poHisto != nullptr )
              {
                  int nOffset = poHisto->GetIntField("columnDataPtr");
                  // Write out histogram data.
                  char *pszWork = pszBinValues;
  
                  // Check whether histogram counts were written as int or double
                  bool bCountIsInt = true;
                  const char *pszDataType = poHisto->GetStringField("dataType");
                  if( STARTS_WITH_CI(pszDataType, "real") )
                  {
                      bCountIsInt = false;
                  }
                  for( int nBin = 0; nBin < nNumBins; ++nBin )
                  {
                      char *pszEnd = strchr(pszWork, '|');
                      if( pszEnd != nullptr )
                      {
                          *pszEnd = 0;
                          if( bCountIsInt )
                          {
                              // Histogram counts were written as ints, so
                              // re-write them the same way.
                              bRet &= VSIFSeekL(hHFA->fp, nOffset + 4 * nBin,
                                                SEEK_SET) >= 0;
                              int nValue = atoi(pszWork);
                              HFAStandard(4, &nValue);
                              bRet &=
                                  VSIFWriteL((void *)&nValue, 4, 1, hHFA->fp) > 0;
                          }
                          else
                          {
                              // Histogram were written as doubles, as is now the
                              // default behaviour.
                              bRet &= VSIFSeekL(hHFA->fp, nOffset + 8 * nBin,
                                                SEEK_SET) >= 0;
                              double nValue = CPLAtof(pszWork);
                              HFAStandard(8, &nValue);
                              bRet &=
                                  VSIFWriteL((void *)&nValue, 8, 1, hHFA->fp) > 0;
                          }
                          pszWork = pszEnd + 1;
                      }
                  }
              }
          }
          CPLFree(pszBinValues);
      }
  
      // If we created a statistics node then try to create a
      // StatisticsParameters node too.
      if( bCreatedStatistics )
      {
          HFAEntry *poEntry =
              HFAEntry::New(hHFA, "StatisticsParameters",
                            "Eimg_StatisticsParameters830", poNode);
  
          poEntry->MakeData(70);
          // poEntry->SetStringField( "BinFunction.binFunctionType", "linear" );
  
          poEntry->SetIntField("SkipFactorX", 1);
          poEntry->SetIntField("SkipFactorY", 1);
      }
  
      // Write out metadata items without a special place.
      if( bRet && CSLCount(papszGDALMD) != 0 )
      {
          CPLErr eErr = HFASetGDALMetadata(hHFA, nBand, papszGDALMD);
  
          CSLDestroy(papszGDALMD);
          return eErr;
      }
      else
      {
          CSLDestroy(papszGDALMD);
          return CE_Failure;
      }
  }
  
  /************************************************************************/
  /*                         HFAGetIGEFilename()                          */
  /*                                                                      */
  /*      Returns the .ige filename if one is associated with this        */
  /*      object.  For files not newly created we need to scan the        */
  /*      bands for spill files.  Presumably there will only be one.      */
  /*                                                                      */
  /*      NOTE: Returns full path, not just the filename portion.         */
  /************************************************************************/
  
  const char *HFAGetIGEFilename( HFAHandle hHFA )
  
  {
      if( hHFA->pszIGEFilename == nullptr )
      {
          std::vector<HFAEntry *> apoDMSList =
              hHFA->poRoot->FindChildren(nullptr, "ImgExternalRaster");
  
          HFAEntry *poDMS = apoDMSList.empty() ? nullptr : apoDMSList[0];
  
          // Get the IGE filename from if we have an ExternalRasterDMS.
          if( poDMS )
          {
              const char *pszRawFilename =
                  poDMS->GetStringField("fileName.string");
  
              if( pszRawFilename != nullptr )
              {
                  VSIStatBufL sStatBuf;
                  CPLString osFullFilename =
                      CPLFormFilename(hHFA->pszPath, pszRawFilename, nullptr);
  
                  if( VSIStatL(osFullFilename, &sStatBuf) != 0 )
                  {
                      const CPLString osExtension =
                          CPLGetExtension(pszRawFilename);
                      const CPLString osBasename =
                          CPLGetBasename(hHFA->pszFilename);
                      osFullFilename =
                          CPLFormFilename(hHFA->pszPath, osBasename, osExtension);
  
                      if( VSIStatL(osFullFilename, &sStatBuf) == 0 )
                          hHFA->pszIGEFilename = CPLStrdup(
                              CPLFormFilename(nullptr, osBasename, osExtension));
                      else
                          hHFA->pszIGEFilename = CPLStrdup(pszRawFilename);
                  }
                  else
                  {
                      hHFA->pszIGEFilename = CPLStrdup(pszRawFilename);
                  }
              }
          }
      }
  
      // Return the full filename.
      if( hHFA->pszIGEFilename )
          return CPLFormFilename(hHFA->pszPath, hHFA->pszIGEFilename, nullptr);
  
      return nullptr;
  }
  
  /************************************************************************/
  /*                        HFACreateSpillStack()                         */
  /*                                                                      */
  /*      Create a new stack of raster layers in the spill (.ige)         */
  /*      file.  Create the spill file if it didn't exist before.         */
  /************************************************************************/
  
  bool HFACreateSpillStack( HFAInfo_t *psInfo, int nXSize, int nYSize,
                            int nLayers, int nBlockSize, EPTType eDataType,
                            GIntBig *pnValidFlagsOffset,
                            GIntBig *pnDataOffset )
  
  {
      // Form .ige filename.
      if( nBlockSize <= 0 )
      {
          CPLError(CE_Failure, CPLE_IllegalArg,
                   "HFACreateSpillStack: nBlockXSize < 0");
          return false;
      }
  
      if( psInfo->pszIGEFilename == nullptr )
      {
          if( EQUAL(CPLGetExtension(psInfo->pszFilename), "rrd") )
              psInfo->pszIGEFilename =
                  CPLStrdup( CPLResetExtension(psInfo->pszFilename, "rde") );
          else if( EQUAL(CPLGetExtension(psInfo->pszFilename), "aux") )
              psInfo->pszIGEFilename =
                  CPLStrdup(CPLResetExtension(psInfo->pszFilename, "axe"));
          else
              psInfo->pszIGEFilename =
                  CPLStrdup(CPLResetExtension(psInfo->pszFilename, "ige"));
      }
  
      char *pszFullFilename = CPLStrdup(
          CPLFormFilename(psInfo->pszPath, psInfo->pszIGEFilename, nullptr));
  
      // Try and open it.  If we fail, create it and write the magic header.
      static const char *const pszMagick = "ERDAS_IMG_EXTERNAL_RASTER";
  
      bool bRet = true;
      VSILFILE *fpVSIL = VSIFOpenL(pszFullFilename, "r+b");
      if( fpVSIL == nullptr )
      {
          fpVSIL = VSIFOpenL(pszFullFilename, "w+");
          if( fpVSIL == nullptr )
          {
              CPLError(CE_Failure, CPLE_OpenFailed,
                       "Failed to create spill file %s.\n%s",
                       psInfo->pszIGEFilename, VSIStrerror(errno));
              return false;
          }
  
          bRet &=
              VSIFWriteL((void *)pszMagick, strlen(pszMagick) + 1, 1, fpVSIL) > 0;
      }
  
      CPLFree(pszFullFilename);
  
      // Work out some details about the tiling scheme.
      const int nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
      const int nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
      // const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
      const int nBytesPerBlock =
          (nBlockSize * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
  
      const int nBytesPerRow = (nBlocksPerRow + 7) / 8;
      const int nBlockMapSize = nBytesPerRow * nBlocksPerColumn;
      // const int iFlagsSize = nBlockMapSize + 20;
  
      // Write stack prefix information.
      bRet &= VSIFSeekL(fpVSIL, 0, SEEK_END) >= 0;
  
      GByte bUnknown = 1;
      bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
  
      GInt32 nValue32 = nLayers;
      HFAStandard(4, &nValue32);
      bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
      nValue32 = nXSize;
      HFAStandard(4, &nValue32);
      bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
      nValue32 = nYSize;
      HFAStandard(4, &nValue32);
      bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
      nValue32 = nBlockSize;
      HFAStandard(4, &nValue32);
      bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
      bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
      bUnknown = 3;
      bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
      bUnknown = 0;
      bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
  
      // Write out ValidFlags section(s).
      *pnValidFlagsOffset = VSIFTellL(fpVSIL);
  
      unsigned char *pabyBlockMap =
          static_cast<unsigned char *>(VSI_MALLOC_VERBOSE(nBlockMapSize));
      if( pabyBlockMap == nullptr )
      {
          CPL_IGNORE_RET_VAL(VSIFCloseL(fpVSIL));
          return false;
      }
  
      memset(pabyBlockMap, 0xff, nBlockMapSize);
      for( int iBand = 0; iBand < nLayers; iBand++ )
      {
          nValue32 = 1;  // Unknown
          HFAStandard(4, &nValue32);
          bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
          nValue32 = 0;  // Unknown
          bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
          nValue32 = nBlocksPerColumn;
          HFAStandard(4, &nValue32);
          bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
          nValue32 = nBlocksPerRow;
          HFAStandard(4, &nValue32);
          bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
          nValue32 = 0x30000;  // Unknown
          HFAStandard(4, &nValue32);
          bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
  
          const int iRemainder = nBlocksPerRow % 8;
          CPLDebug("HFACreate",
                   "Block map size %d, bytes per row %d, remainder %d.",
                   nBlockMapSize, nBytesPerRow, iRemainder);
          if( iRemainder )
          {
              for( int i = nBytesPerRow - 1; i < nBlockMapSize; i += nBytesPerRow )
                  pabyBlockMap[i] = static_cast<GByte>((1 << iRemainder) - 1);
          }
  
          bRet &= VSIFWriteL(pabyBlockMap, nBlockMapSize, 1, fpVSIL) > 0;
      }
      CPLFree(pabyBlockMap);
      pabyBlockMap = nullptr;
  
      // Extend the file to account for all the imagery space.
      const GIntBig nTileDataSize = static_cast<GIntBig>(nBytesPerBlock) *
                                    nBlocksPerRow * nBlocksPerColumn * nLayers;
  
      *pnDataOffset = VSIFTellL(fpVSIL);
  
      if( !bRet ||
          VSIFTruncateL(fpVSIL, nTileDataSize + *pnDataOffset) != 0 )
      {
          CPLError(CE_Failure, CPLE_FileIO,
                   "Failed to extend %s to full size (" CPL_FRMT_GIB " bytes), "
                   "likely out of disk space.\n%s",
                   psInfo->pszIGEFilename,
                   nTileDataSize + *pnDataOffset,
                   VSIStrerror(errno));
  
          CPL_IGNORE_RET_VAL(VSIFCloseL( fpVSIL ));
          return false;
      }
  
      if( VSIFCloseL(fpVSIL) != 0 )
          return false;
  
      return true;
  }
  
  /************************************************************************/
  /*                       HFAReadAndValidatePoly()                       */
  /************************************************************************/
  
  static bool HFAReadAndValidatePoly( HFAEntry *poTarget,
                                      const char *pszName,
                                      Efga_Polynomial *psRetPoly )
  
  {
      memset(psRetPoly, 0, sizeof(Efga_Polynomial));
  
      CPLString osFldName;
      osFldName.Printf("%sorder", pszName);
      psRetPoly->order = poTarget->GetIntField(osFldName);
  
      if( psRetPoly->order < 1 || psRetPoly->order > 3 )
          return false;
  
      // Validate that things are in a "well known" form.
      osFldName.Printf("%snumdimtransform", pszName);
      const int numdimtransform = poTarget->GetIntField(osFldName);
  
      osFldName.Printf("%snumdimpolynomial", pszName);
      const int numdimpolynomial = poTarget->GetIntField(osFldName);
  
      osFldName.Printf("%stermcount", pszName);
      const int termcount = poTarget->GetIntField(osFldName);
  
      if( numdimtransform != 2 || numdimpolynomial != 2 )
          return false;
  
      if( (psRetPoly->order == 1 && termcount != 3) ||
          (psRetPoly->order == 2 && termcount != 6) ||
          (psRetPoly->order == 3 && termcount != 10) )
          return false;
  
      // We don't check the exponent organization for now.  Hopefully
      // it is always standard.
  
      // Get coefficients.
      for( int i = 0; i < termcount*2 - 2; i++ )
      {
          osFldName.Printf("%spolycoefmtx[%d]", pszName, i);
          psRetPoly->polycoefmtx[i] = poTarget->GetDoubleField(osFldName);
      }
  
      for( int i = 0; i < 2; i++ )
      {
          osFldName.Printf("%spolycoefvector[%d]", pszName, i);
          psRetPoly->polycoefvector[i] = poTarget->GetDoubleField(osFldName);
      }
  
      return true;
  }
  
  /************************************************************************/
  /*                         HFAReadXFormStack()                          */
  /************************************************************************/
  
  int HFAReadXFormStack( HFAHandle hHFA,
                         Efga_Polynomial **ppasPolyListForward,
                         Efga_Polynomial **ppasPolyListReverse )
  
  {
      if( hHFA->nBands == 0 )
          return 0;
  
      // Get the HFA node.
      HFAEntry *poXFormHeader =
          hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm");
      if( poXFormHeader == nullptr )
          return 0;
  
      // Loop over children, collecting XForms.
      int nStepCount = 0;
      *ppasPolyListForward = nullptr;
      *ppasPolyListReverse = nullptr;
  
      for( HFAEntry *poXForm = poXFormHeader->GetChild();
           poXForm != nullptr;
           poXForm = poXForm->GetNext() )
      {
          bool bSuccess = false;
          Efga_Polynomial sForward;
          Efga_Polynomial sReverse;
          memset(&sForward, 0, sizeof(sForward));
          memset(&sReverse, 0, sizeof(sReverse));
  
          if( EQUAL(poXForm->GetType(), "Efga_Polynomial") )
          {
              bSuccess = HFAReadAndValidatePoly(poXForm, "", &sForward);
  
              if( bSuccess )
              {
                  double adfGT[6] = {
                      sForward.polycoefvector[0],
                      sForward.polycoefmtx[0],
                      sForward.polycoefmtx[2],
                      sForward.polycoefvector[1],
                      sForward.polycoefmtx[1],
                      sForward.polycoefmtx[3]
                  };
  
                  double adfInvGT[6] = {};
                  bSuccess = HFAInvGeoTransform(adfGT, adfInvGT);
                  if( !bSuccess )
                      memset(adfInvGT, 0, sizeof(adfInvGT));
  
                  sReverse.order = sForward.order;
                  sReverse.polycoefvector[0] = adfInvGT[0];
                  sReverse.polycoefmtx[0] = adfInvGT[1];
                  sReverse.polycoefmtx[2] = adfInvGT[2];
                  sReverse.polycoefvector[1] = adfInvGT[3];
                  sReverse.polycoefmtx[1] = adfInvGT[4];
                  sReverse.polycoefmtx[3] = adfInvGT[5];
              }
          }
          else if( EQUAL(poXForm->GetType(), "GM_PolyPair") )
          {
              bSuccess =
                  HFAReadAndValidatePoly(poXForm, "forward.", &sForward) &&
                  HFAReadAndValidatePoly(poXForm, "reverse.", &sReverse);
          }
  
          if( bSuccess )
          {
              nStepCount++;
              *ppasPolyListForward = static_cast<Efga_Polynomial *>(CPLRealloc(
                  *ppasPolyListForward, sizeof(Efga_Polynomial) * nStepCount));
              memcpy(*ppasPolyListForward + nStepCount - 1, &sForward,
                     sizeof(sForward));
  
              *ppasPolyListReverse = static_cast<Efga_Polynomial *>(CPLRealloc(
                  *ppasPolyListReverse, sizeof(Efga_Polynomial) * nStepCount));
              memcpy(*ppasPolyListReverse + nStepCount - 1, &sReverse,
                     sizeof(sReverse));
          }
      }
  
      return nStepCount;
  }
  
  /************************************************************************/
  /*                       HFAEvaluateXFormStack()                        */
  /************************************************************************/
  
  int HFAEvaluateXFormStack( int nStepCount, int bForward,
                             Efga_Polynomial *pasPolyList,
                             double *pdfX, double *pdfY )
  
  {
      for( int iStep = 0; iStep < nStepCount; iStep++ )
      {
          const Efga_Polynomial *psStep =
              bForward ? pasPolyList + iStep
                       : pasPolyList + nStepCount - iStep - 1;
  
          double dfXOut = 0.0;
          double dfYOut = 0.0;
          if( psStep->order == 1 )
          {
              dfXOut = psStep->polycoefvector[0] +
                       psStep->polycoefmtx[0] * *pdfX +
                       psStep->polycoefmtx[2] * *pdfY;
  
              dfYOut = psStep->polycoefvector[1] +
                       psStep->polycoefmtx[1] * *pdfX +
                       psStep->polycoefmtx[3] * *pdfY;
  
              *pdfX = dfXOut;
              *pdfY = dfYOut;
          }
          else if( psStep->order == 2 )
          {
              dfXOut = psStep->polycoefvector[0] +
                       psStep->polycoefmtx[0] * *pdfX +
                       psStep->polycoefmtx[2] * *pdfY +
                       psStep->polycoefmtx[4] * *pdfX * *pdfX +
                       psStep->polycoefmtx[6] * *pdfX * *pdfY +
                       psStep->polycoefmtx[8] * *pdfY * *pdfY;
              dfYOut = psStep->polycoefvector[1] +
                       psStep->polycoefmtx[1] * *pdfX +
                       psStep->polycoefmtx[3] * *pdfY +
                       psStep->polycoefmtx[5] * *pdfX * *pdfX +
                       psStep->polycoefmtx[7] * *pdfX * *pdfY +
                       psStep->polycoefmtx[9] * *pdfY * *pdfY;
  
              *pdfX = dfXOut;
              *pdfY = dfYOut;
          }
          else if( psStep->order == 3 )
          {
              dfXOut = psStep->polycoefvector[0] +
                       psStep->polycoefmtx[0] * *pdfX +
                       psStep->polycoefmtx[2] * *pdfY +
                       psStep->polycoefmtx[4] * *pdfX * *pdfX +
                       psStep->polycoefmtx[6] * *pdfX * *pdfY +
                       psStep->polycoefmtx[8] * *pdfY * *pdfY +
                       psStep->polycoefmtx[10] * *pdfX * *pdfX * *pdfX +
                       psStep->polycoefmtx[12] * *pdfX * *pdfX * *pdfY +
                       psStep->polycoefmtx[14] * *pdfX * *pdfY * *pdfY +
                       psStep->polycoefmtx[16] * *pdfY * *pdfY * *pdfY;
              dfYOut = psStep->polycoefvector[1] +
                       psStep->polycoefmtx[1] * *pdfX +
                       psStep->polycoefmtx[3] * *pdfY +
                       psStep->polycoefmtx[5] * *pdfX * *pdfX +
                       psStep->polycoefmtx[7] * *pdfX * *pdfY +
                       psStep->polycoefmtx[9] * *pdfY * *pdfY +
                       psStep->polycoefmtx[11] * *pdfX * *pdfX * *pdfX +
                       psStep->polycoefmtx[13] * *pdfX * *pdfX * *pdfY +
                       psStep->polycoefmtx[15] * *pdfX * *pdfY * *pdfY +
                       psStep->polycoefmtx[17] * *pdfY * *pdfY * *pdfY;
  
              *pdfX = dfXOut;
              *pdfY = dfYOut;
          }
          else
              return FALSE;
      }
  
      return TRUE;
  }
  
  /************************************************************************/
  /*                         HFAWriteXFormStack()                         */
  /************************************************************************/
  
  CPLErr HFAWriteXFormStack( HFAHandle hHFA, int nBand, int nXFormCount,
                             Efga_Polynomial **ppasPolyListForward,
                             Efga_Polynomial **ppasPolyListReverse )
  
  {
      if( nXFormCount == 0 )
          return CE_None;
  
      if( ppasPolyListForward[0]->order != 1 )
      {
          CPLError(
              CE_Failure, CPLE_AppDefined,
              "For now HFAWriteXFormStack() only supports order 1 polynomials");
          return CE_Failure;
      }
  
      if( nBand < 0 || nBand > hHFA->nBands )
          return CE_Failure;
  
      // If no band number is provided, operate on all bands.
      if( nBand == 0 )
      {
          for( nBand = 1; nBand <= hHFA->nBands; nBand++ )
          {
              CPLErr eErr =
                  HFAWriteXFormStack(hHFA, nBand, nXFormCount,
                                     ppasPolyListForward,
                                     ppasPolyListReverse);
              if( eErr != CE_None )
                  return eErr;
          }
  
          return CE_None;
      }
  
      // Fetch our band node.
      HFAEntry *poBandNode = hHFA->papoBand[nBand - 1]->poNode;
      HFAEntry *poXFormHeader = poBandNode->GetNamedChild("MapToPixelXForm");
      if( poXFormHeader == nullptr )
      {
          poXFormHeader = HFAEntry::New(hHFA, "MapToPixelXForm",
                                        "Exfr_GenericXFormHeader", poBandNode);
          poXFormHeader->MakeData(23);
          poXFormHeader->SetPosition();
          poXFormHeader->SetStringField("titleList.string", "Affine");
      }
  
      // Loop over XForms.
      for( int iXForm = 0; iXForm < nXFormCount; iXForm++ )
      {
          Efga_Polynomial *psForward = *ppasPolyListForward + iXForm;
          CPLString osXFormName;
          osXFormName.Printf("XForm%d", iXForm);
  
          HFAEntry *poXForm = poXFormHeader->GetNamedChild(osXFormName);
  
          if( poXForm == nullptr )
          {
              poXForm = HFAEntry::New(hHFA, osXFormName, "Efga_Polynomial",
                                      poXFormHeader);
              poXForm->MakeData(136);
              poXForm->SetPosition();
          }
  
          poXForm->SetIntField("order", 1);
          poXForm->SetIntField("numdimtransform", 2);
          poXForm->SetIntField("numdimpolynomial", 2);
          poXForm->SetIntField("termcount", 3);
          poXForm->SetIntField("exponentlist[0]", 0);
          poXForm->SetIntField("exponentlist[1]", 0);
          poXForm->SetIntField("exponentlist[2]", 1);
          poXForm->SetIntField("exponentlist[3]", 0);
          poXForm->SetIntField("exponentlist[4]", 0);
          poXForm->SetIntField("exponentlist[5]", 1);
  
          poXForm->SetIntField("polycoefmtx[-3]", EPT_f64);
          poXForm->SetIntField("polycoefmtx[-2]", 2);
          poXForm->SetIntField("polycoefmtx[-1]", 2);
          poXForm->SetDoubleField("polycoefmtx[0]", psForward->polycoefmtx[0]);
          poXForm->SetDoubleField("polycoefmtx[1]", psForward->polycoefmtx[1]);
          poXForm->SetDoubleField("polycoefmtx[2]", psForward->polycoefmtx[2]);
          poXForm->SetDoubleField("polycoefmtx[3]", psForward->polycoefmtx[3]);
  
          poXForm->SetIntField("polycoefvector[-3]", EPT_f64);
          poXForm->SetIntField("polycoefvector[-2]", 1);
          poXForm->SetIntField("polycoefvector[-1]", 2);
          poXForm->SetDoubleField("polycoefvector[0]",
                                  psForward->polycoefvector[0]);
          poXForm->SetDoubleField("polycoefvector[1]",
                                  psForward->polycoefvector[1]);
      }
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                         HFAReadCameraModel()                         */
  /************************************************************************/
  
  char **HFAReadCameraModel( HFAHandle hHFA )
  
  {
      if( hHFA->nBands == 0 )
          return nullptr;
  
      // Get the camera model node, and confirm its type.
      HFAEntry *poXForm =
          hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm0");
      if( poXForm == nullptr )
          return nullptr;
  
      if( !EQUAL(poXForm->GetType(), "Camera_ModelX") )
          return nullptr;
  
      // Convert the values to metadata.
      char **papszMD = nullptr;
      static const char * const apszFields[] = {
          "direction", "refType", "demsource", "PhotoDirection", "RotationSystem",
          "demfilename", "demzunits",
          "forSrcAffine[0]", "forSrcAffine[1]", "forSrcAffine[2]",
          "forSrcAffine[3]", "forSrcAffine[4]", "forSrcAffine[5]",
          "forDstAffine[0]", "forDstAffine[1]", "forDstAffine[2]",
          "forDstAffine[3]", "forDstAffine[4]", "forDstAffine[5]",
          "invSrcAffine[0]", "invSrcAffine[1]", "invSrcAffine[2]",
          "invSrcAffine[3]", "invSrcAffine[4]", "invSrcAffine[5]",
          "invDstAffine[0]", "invDstAffine[1]", "invDstAffine[2]",
          "invDstAffine[3]", "invDstAffine[4]", "invDstAffine[5]",
          "z_mean", "lat0", "lon0",
          "coeffs[0]", "coeffs[1]", "coeffs[2]",
          "coeffs[3]", "coeffs[4]", "coeffs[5]",
          "coeffs[6]", "coeffs[7]", "coeffs[8]",
          "LensDistortion[0]", "LensDistortion[1]", "LensDistortion[2]",
          nullptr };
  
      const char *pszValue = nullptr;
      for( int i = 0; apszFields[i] != nullptr; i++ )
      {
          pszValue = poXForm->GetStringField(apszFields[i]);
          if( pszValue == nullptr )
              pszValue = "";
  
          papszMD = CSLSetNameValue(papszMD, apszFields[i], pszValue);
      }
  
      // Create a pseudo-entry for the MIFObject with the outputProjection.
      HFAEntry *poProjInfo =
          HFAEntry::BuildEntryFromMIFObject(poXForm, "outputProjection");
      if( poProjInfo )
      {
          // Fetch the datum.
          Eprj_Datum sDatum;
  
          memset(&sDatum, 0, sizeof(sDatum));
  
          sDatum.datumname =
              (char *)poProjInfo->GetStringField("earthModel.datum.datumname");
  
          const int nDatumType = poProjInfo->GetIntField("earthModel.datum.type");
          if( nDatumType < 0 || nDatumType > EPRJ_DATUM_NONE )
          {
              CPLDebug("HFA", "Invalid value for datum type: %d", nDatumType);
              sDatum.type = EPRJ_DATUM_NONE;
          }
          else
          {
              sDatum.type = static_cast<Eprj_DatumType>(nDatumType);
          }
  
          for( int i = 0; i < 7; i++ )
          {
              char szFieldName[60] = {};
  
              snprintf(szFieldName, sizeof(szFieldName),
                       "earthModel.datum.params[%d]", i);
              sDatum.params[i] = poProjInfo->GetDoubleField(szFieldName);
          }
  
          sDatum.gridname =
              (char *)poProjInfo->GetStringField("earthModel.datum.gridname");
  
          // Fetch the projection parameters.
          Eprj_ProParameters sPro;
  
          memset(&sPro, 0, sizeof(sPro));
  
          sPro.proType =
              (Eprj_ProType)poProjInfo->GetIntField("projectionObject.proType");
          sPro.proNumber = poProjInfo->GetIntField("projectionObject.proNumber");
          sPro.proExeName =
              (char *)poProjInfo->GetStringField("projectionObject.proExeName");
          sPro.proName =
              (char *)poProjInfo->GetStringField("projectionObject.proName");
          sPro.proZone = poProjInfo->GetIntField("projectionObject.proZone");
  
          for( int i = 0; i < 15; i++ )
          {
              char szFieldName[40] = {};
  
              snprintf(szFieldName, sizeof(szFieldName),
                       "projectionObject.proParams[%d]", i);
              sPro.proParams[i] = poProjInfo->GetDoubleField(szFieldName);
          }
  
          // Fetch the spheroid.
          sPro.proSpheroid.sphereName = (char *)poProjInfo->GetStringField(
              "earthModel.proSpheroid.sphereName");
          sPro.proSpheroid.a =
              poProjInfo->GetDoubleField("earthModel.proSpheroid.a");
          sPro.proSpheroid.b =
              poProjInfo->GetDoubleField("earthModel.proSpheroid.b");
          sPro.proSpheroid.eSquared =
              poProjInfo->GetDoubleField("earthModel.proSpheroid.eSquared");
          sPro.proSpheroid.radius =
              poProjInfo->GetDoubleField("earthModel.proSpheroid.radius");
  
          // Fetch the projection info.
          // poProjInfo->DumpFieldValues( stdout, "" );
  
          char *pszProjection = HFAPCSStructToWKT(&sDatum, &sPro, nullptr, nullptr);
  
          if( pszProjection )
          {
              papszMD =
                  CSLSetNameValue(papszMD, "outputProjection", pszProjection);
              CPLFree(pszProjection);
          }
  
          delete poProjInfo;
      }
  
      // Fetch the horizontal units.
      pszValue = poXForm->GetStringField("outputHorizontalUnits.string");
      if( pszValue == nullptr )
          pszValue = "";
  
      papszMD = CSLSetNameValue(papszMD, "outputHorizontalUnits", pszValue);
  
      // Fetch the elevationinfo.
      HFAEntry *poElevInfo =
          HFAEntry::BuildEntryFromMIFObject(poXForm, "outputElevationInfo");
      if( poElevInfo )
      {
          // poElevInfo->DumpFieldValues( stdout, "" );
  
          if( poElevInfo->GetDataSize() != 0 )
          {
              static const char *const apszEFields[] = {
                  "verticalDatum.datumname",
                  "verticalDatum.type",
                  "elevationUnit",
                  "elevationType",
                  nullptr };
  
              for( int i = 0; apszEFields[i] != nullptr; i++ )
              {
                  pszValue = poElevInfo->GetStringField(apszEFields[i]);
                  if( pszValue == nullptr )
                      pszValue = "";
  
                  papszMD = CSLSetNameValue(papszMD, apszEFields[i], pszValue);
              }
          }
  
          delete poElevInfo;
      }
  
      return papszMD;
  }
  
  /************************************************************************/
  /*                         HFAReadElevationUnit()                       */
  /************************************************************************/
  
  const char *HFAReadElevationUnit( HFAHandle hHFA, int iBand )
  {
      if( hHFA->nBands <= iBand )
          return nullptr;
  
      HFABand *poBand(hHFA->papoBand[iBand]);
      if( poBand == nullptr || poBand->poNode == nullptr )
      {
          return nullptr;
      }
      HFAEntry *poElevInfo = poBand->poNode->GetNamedChild("Elevation_Info");
      if( poElevInfo == nullptr )
      {
          return nullptr;
      }
      return poElevInfo->GetStringField("elevationUnit");
  }
  
  /************************************************************************/
  /*                         HFASetGeoTransform()                         */
  /*                                                                      */
  /*      Set a MapInformation and XForm block.  Allows for rotated       */
  /*      and shared geotransforms.                                       */
  /************************************************************************/
  
  CPLErr HFASetGeoTransform( HFAHandle hHFA,
                             const char *pszProName,
                             const char *pszUnits,
                             double *padfGeoTransform )
  
  {
      // Write MapInformation.
      for( int nBand = 1; nBand <= hHFA->nBands; nBand++ )
      {
          HFAEntry *poBandNode = hHFA->papoBand[nBand - 1]->poNode;
  
          HFAEntry *poMI = poBandNode->GetNamedChild("MapInformation");
          if( poMI == nullptr )
          {
              poMI = HFAEntry::New(hHFA, "MapInformation", "Eimg_MapInformation",
                                   poBandNode);
              poMI->MakeData(
                  static_cast<int>(18 + strlen(pszProName) + strlen(pszUnits)));
              poMI->SetPosition();
          }
  
          poMI->SetStringField("projection.string", pszProName);
          poMI->SetStringField("units.string", pszUnits);
      }
  
      // Write XForm.
      double adfAdjTransform[6] = {};
  
      // Offset by half pixel.
  
      memcpy(adfAdjTransform, padfGeoTransform, sizeof(double) * 6);
      adfAdjTransform[0] += adfAdjTransform[1] * 0.5;
      adfAdjTransform[0] += adfAdjTransform[2] * 0.5;
      adfAdjTransform[3] += adfAdjTransform[4] * 0.5;
      adfAdjTransform[3] += adfAdjTransform[5] * 0.5;
  
      // Invert.
      double adfRevTransform[6] = {};
      if( !HFAInvGeoTransform( adfAdjTransform, adfRevTransform ) )
          memset(adfRevTransform, 0, sizeof(adfRevTransform));
  
      // Assign to polynomial object.
  
      Efga_Polynomial sForward;
      memset(&sForward, 0, sizeof(sForward));
      Efga_Polynomial *psForward = &sForward;
      sForward.order = 1;
      sForward.polycoefvector[0] = adfRevTransform[0];
      sForward.polycoefmtx[0] = adfRevTransform[1];
      sForward.polycoefmtx[1] = adfRevTransform[4];
      sForward.polycoefvector[1] = adfRevTransform[3];
      sForward.polycoefmtx[2] = adfRevTransform[2];
      sForward.polycoefmtx[3] = adfRevTransform[5];
  
      Efga_Polynomial sReverse = sForward;
      Efga_Polynomial *psReverse = &sReverse;
  
      return HFAWriteXFormStack(hHFA, 0, 1, &psForward, &psReverse);
  }
  
  /************************************************************************/
  /*                        HFARenameReferences()                         */
  /*                                                                      */
  /*      Rename references in this .img file from the old basename to    */
  /*      a new basename.  This should be passed on to .aux and .rrd      */
  /*      files and should include references to .aux, .rrd and .ige.     */
  /************************************************************************/
  
  CPLErr HFARenameReferences( HFAHandle hHFA,
                              const char *pszNewBase,
                              const char *pszOldBase )
  
  {
      // Handle RRDNamesList updates.
      std::vector<HFAEntry *> apoNodeList =
          hHFA->poRoot->FindChildren("RRDNamesList", nullptr);
  
      for( size_t iNode = 0; iNode < apoNodeList.size(); iNode++ )
      {
          HFAEntry *poRRDNL = apoNodeList[iNode];
          std::vector<CPLString> aosNL;
  
          // Collect all the existing names.
          const int nNameCount = poRRDNL->GetFieldCount("nameList");
  
          CPLString osAlgorithm = poRRDNL->GetStringField("algorithm.string");
          for( int i = 0; i < nNameCount; i++ )
          {
              CPLString osFN;
              osFN.Printf("nameList[%d].string", i);
              aosNL.push_back(poRRDNL->GetStringField(osFN));
          }
  
          // Adjust the names to the new form.
          for( int i = 0; i < nNameCount; i++ )
          {
              if( strncmp(aosNL[i], pszOldBase, strlen(pszOldBase)) == 0 )
              {
                  CPLString osNew = pszNewBase;
                  osNew += aosNL[i].c_str() + strlen(pszOldBase);
                  aosNL[i] = osNew;
              }
          }
  
          // Try to make sure the RRDNamesList is big enough to hold the
          // adjusted name list.
          if( strlen(pszNewBase) > strlen(pszOldBase) )
          {
              CPLDebug("HFA", "Growing RRDNamesList to hold new names");
              poRRDNL->MakeData(static_cast<int>(
                  poRRDNL->GetDataSize() +
                  nNameCount * (strlen(pszNewBase) - strlen(pszOldBase))));
          }
  
          // Initialize the whole thing to zeros for a clean start.
          memset(poRRDNL->GetData(), 0, poRRDNL->GetDataSize());
  
          // Write the updates back to the file.
          poRRDNL->SetStringField("algorithm.string", osAlgorithm);
          for( int i = 0; i < nNameCount; i++ )
          {
              CPLString osFN;
              osFN.Printf("nameList[%d].string", i);
              poRRDNL->SetStringField(osFN, aosNL[i]);
          }
      }
  
      // Spill file references.
      apoNodeList =
          hHFA->poRoot->FindChildren("ExternalRasterDMS", "ImgExternalRaster");
  
      for( size_t iNode = 0; iNode < apoNodeList.size(); iNode++ )
      {
          HFAEntry *poERDMS = apoNodeList[iNode];
  
          if( poERDMS == nullptr )
              continue;
  
          // Fetch all existing values.
          CPLString osFileName = poERDMS->GetStringField("fileName.string");
  
          GInt32 anValidFlagsOffset[2] = {
              poERDMS->GetIntField("layerStackValidFlagsOffset[0]"),
              poERDMS->GetIntField("layerStackValidFlagsOffset[1]")
          };
  
          GInt32 anStackDataOffset[2] = {
              poERDMS->GetIntField("layerStackDataOffset[0]"),
              poERDMS->GetIntField("layerStackDataOffset[1]")
          };
  
          const GInt32 nStackCount = poERDMS->GetIntField("layerStackCount");
          const GInt32 nStackIndex = poERDMS->GetIntField("layerStackIndex");
  
          // Update the filename.
          if( strncmp(osFileName, pszOldBase, strlen(pszOldBase)) == 0 )
          {
              CPLString osNew = pszNewBase;
              osNew += osFileName.c_str() + strlen(pszOldBase);
              osFileName = osNew;
          }
  
          // Grow the node if needed.
          if( strlen(pszNewBase) > strlen(pszOldBase) )
          {
              CPLDebug("HFA", "Growing ExternalRasterDMS to hold new names");
              poERDMS->MakeData(
                  static_cast<int>(poERDMS->GetDataSize() +
                                   (strlen(pszNewBase) - strlen(pszOldBase))));
          }
  
          // Initialize the whole thing to zeros for a clean start.
          memset(poERDMS->GetData(), 0, poERDMS->GetDataSize());
  
          // Write it all out again, this may change the size of the node.
          poERDMS->SetStringField("fileName.string", osFileName);
          poERDMS->SetIntField("layerStackValidFlagsOffset[0]",
                               anValidFlagsOffset[0]);
          poERDMS->SetIntField("layerStackValidFlagsOffset[1]",
                               anValidFlagsOffset[1]);
  
          poERDMS->SetIntField("layerStackDataOffset[0]", anStackDataOffset[0]);
          poERDMS->SetIntField("layerStackDataOffset[1]", anStackDataOffset[1]);
  
          poERDMS->SetIntField("layerStackCount", nStackCount);
          poERDMS->SetIntField("layerStackIndex", nStackIndex);
      }
  
      // DependentFile.
      apoNodeList =
          hHFA->poRoot->FindChildren("DependentFile", "Eimg_DependentFile");
  
      for( size_t iNode = 0; iNode < apoNodeList.size(); iNode++ )
      {
          CPLString osFileName =
              apoNodeList[iNode]->GetStringField("dependent.string");
  
          // Grow the node if needed.
          if( strlen(pszNewBase) > strlen(pszOldBase) )
          {
              CPLDebug("HFA", "Growing DependentFile to hold new names");
              apoNodeList[iNode]->MakeData(
                  static_cast<int>(apoNodeList[iNode]->GetDataSize() +
                                   (strlen(pszNewBase) - strlen(pszOldBase))));
          }
  
          // Update the filename.
          if( strncmp(osFileName, pszOldBase, strlen(pszOldBase)) == 0 )
          {
              CPLString osNew = pszNewBase;
              osNew += osFileName.c_str() + strlen(pszOldBase);
              osFileName = osNew;
          }
  
          apoNodeList[iNode]->SetStringField("dependent.string", osFileName);
      }
  
      return CE_None;
  }
  