LCOV - code coverage report
Current view: directory - frmts/hf2 - hf2dataset.cpp (source / functions) Found Hit Coverage
Test: gdal_filtered.info Lines: 507 438 86.4 %
Date: 2012-04-28 Functions: 22 15 68.2 %

       1                 : /******************************************************************************
       2                 :  * $Id: hf2dataset.cpp 21680 2011-02-11 21:12:07Z warmerdam $
       3                 :  *
       4                 :  * Project:  HF2 driver
       5                 :  * Purpose:  GDALDataset driver for HF2/HFZ dataset.
       6                 :  * Author:   Even Rouault, <even dot rouault at mines dash paris dot org>
       7                 :  *
       8                 :  ******************************************************************************
       9                 :  * Copyright (c) 2010, Even Rouault
      10                 :  *
      11                 :  * Permission is hereby granted, free of charge, to any person obtaining a
      12                 :  * copy of this software and associated documentation files (the "Software"),
      13                 :  * to deal in the Software without restriction, including without limitation
      14                 :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      15                 :  * and/or sell copies of the Software, and to permit persons to whom the
      16                 :  * Software is furnished to do so, subject to the following conditions:
      17                 :  *
      18                 :  * The above copyright notice and this permission notice shall be included
      19                 :  * in all copies or substantial portions of the Software.
      20                 :  *
      21                 :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      22                 :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      23                 :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      24                 :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      25                 :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      26                 :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      27                 :  * DEALINGS IN THE SOFTWARE.
      28                 :  ****************************************************************************/
      29                 : 
      30                 : #include "cpl_string.h"
      31                 : #include "gdal_pam.h"
      32                 : #include "ogr_spatialref.h"
      33                 : 
      34                 : CPL_CVSID("$Id: hf2dataset.cpp 21680 2011-02-11 21:12:07Z warmerdam $");
      35                 : 
      36                 : CPL_C_START
      37                 : void    GDALRegister_HF2(void);
      38                 : CPL_C_END
      39                 : 
      40                 : /************************************************************************/
      41                 : /* ==================================================================== */
      42                 : /*                              HF2Dataset                              */
      43                 : /* ==================================================================== */
      44                 : /************************************************************************/
      45                 : 
      46                 : class HF2RasterBand;
      47                 : 
      48                 : class HF2Dataset : public GDALPamDataset
      49                 : {
      50                 :     friend class HF2RasterBand;
      51                 :     
      52                 :     VSILFILE   *fp;
      53                 :     double      adfGeoTransform[6];
      54                 :     char       *pszWKT;
      55                 :     vsi_l_offset    *panBlockOffset;
      56                 : 
      57                 :     int         nTileSize;
      58                 :     int         bHasLoaderBlockMap;
      59                 :     int         LoadBlockMap();
      60                 : 
      61                 :   public:
      62                 :                  HF2Dataset();
      63                 :     virtual     ~HF2Dataset();
      64                 :     
      65                 :     virtual CPLErr GetGeoTransform( double * );
      66                 :     virtual const char* GetProjectionRef();
      67                 :     
      68                 :     static GDALDataset *Open( GDALOpenInfo * );
      69                 :     static int          Identify( GDALOpenInfo * );
      70                 :     static GDALDataset *CreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
      71                 :                                     int bStrict, char ** papszOptions, 
      72                 :                                     GDALProgressFunc pfnProgress, void * pProgressData );
      73                 : };
      74                 : 
      75                 : /************************************************************************/
      76                 : /* ==================================================================== */
      77                 : /*                            HF2RasterBand                             */
      78                 : /* ==================================================================== */
      79                 : /************************************************************************/
      80                 : 
      81                 : class HF2RasterBand : public GDALPamRasterBand
      82                 : {
      83                 :     friend class HF2Dataset;
      84                 : 
      85                 :     float*  pafBlockData;
      86                 :     int     nLastBlockYOff;
      87                 : 
      88                 :   public:
      89                 : 
      90                 :                 HF2RasterBand( HF2Dataset *, int, GDALDataType );
      91                 :                ~HF2RasterBand();
      92                 : 
      93                 :     virtual CPLErr IReadBlock( int, int, void * );
      94                 : };
      95                 : 
      96                 : 
      97                 : /************************************************************************/
      98                 : /*                           HF2RasterBand()                            */
      99                 : /************************************************************************/
     100                 : 
     101              46 : HF2RasterBand::HF2RasterBand( HF2Dataset *poDS, int nBand, GDALDataType eDT )
     102                 : 
     103                 : {
     104              46 :     this->poDS = poDS;
     105              46 :     this->nBand = nBand;
     106                 : 
     107              46 :     eDataType = eDT;
     108                 : 
     109              46 :     nBlockXSize = poDS->nTileSize;
     110              46 :     nBlockYSize = 1;
     111                 : 
     112              46 :     pafBlockData = NULL;
     113              46 :     nLastBlockYOff = -1;
     114              46 : }
     115                 : 
     116                 : /************************************************************************/
     117                 : /*                          ~HF2RasterBand()                            */
     118                 : /************************************************************************/
     119                 : 
     120              46 : HF2RasterBand::~HF2RasterBand()
     121                 : {
     122              46 :     CPLFree(pafBlockData);
     123              46 : }
     124                 : 
     125                 : /************************************************************************/
     126                 : /*                             IReadBlock()                             */
     127                 : /************************************************************************/
     128                 : 
     129            1528 : CPLErr HF2RasterBand::IReadBlock( int nBlockXOff, int nLineYOff,
     130                 :                                   void * pImage )
     131                 : 
     132                 : {
     133            1528 :     HF2Dataset *poGDS = (HF2Dataset *) poDS;
     134                 : 
     135            1528 :     int nXBlocks = (nRasterXSize + nBlockXSize - 1) / nBlockXSize;
     136            1528 :     int nYBlocks = (nRasterYSize + nBlockXSize - 1) / nBlockXSize;
     137                 :     
     138            1528 :     if (!poGDS->LoadBlockMap())
     139               0 :         return CE_Failure;
     140                 :     
     141            1528 :     if (pafBlockData == NULL)
     142                 :     {
     143              16 :         pafBlockData = (float*)VSIMalloc3(nXBlocks * sizeof(float), poGDS->nTileSize, poGDS->nTileSize);
     144              16 :         if (pafBlockData == NULL)
     145               0 :             return CE_Failure;
     146                 :     }
     147                 :     
     148            1528 :     nLineYOff = nRasterYSize - 1 - nLineYOff;
     149                 : 
     150            1528 :     int nBlockYOff = nLineYOff / nBlockXSize;
     151            1528 :     int nYOffInTile = nLineYOff % nBlockXSize;
     152                 : 
     153            1528 :     if (nBlockYOff != nLastBlockYOff)
     154                 :     {
     155              20 :         nLastBlockYOff = nBlockYOff;
     156                 : 
     157              20 :         memset(pafBlockData, 0, nXBlocks * sizeof(float) * nBlockXSize * nBlockXSize);
     158                 : 
     159                 :         /* 4 * nBlockXSize is the upper bound */
     160              20 :         void* pabyData = CPLMalloc( 4 * nBlockXSize );
     161                 : 
     162                 :         int nxoff;
     163              48 :         for(nxoff = 0; nxoff < nXBlocks; nxoff++)
     164                 :         {
     165              28 :             VSIFSeekL(poGDS->fp, poGDS->panBlockOffset[(nYBlocks - 1 - nBlockYOff) * nXBlocks + nxoff], SEEK_SET);
     166                 :             float fScale, fOff;
     167              28 :             VSIFReadL(&fScale, 4, 1, poGDS->fp);
     168              28 :             VSIFReadL(&fOff, 4, 1, poGDS->fp);
     169                 :             CPL_LSBPTR32(&fScale);
     170                 :             CPL_LSBPTR32(&fOff);
     171                 :     
     172              28 :             int nTileWidth = MIN(nBlockXSize, nRasterXSize - nxoff * nBlockXSize);
     173              28 :             int nTileHeight = MIN(nBlockXSize, nRasterYSize - nBlockYOff * nBlockXSize);
     174                 :             
     175                 :             int j;
     176            1556 :             for(j=0;j<nTileHeight;j++)
     177                 :             {
     178                 :                 GByte nWordSize;
     179            1528 :                 VSIFReadL(&nWordSize, 1, 1, poGDS->fp);
     180            1528 :                 if (nWordSize != 1 && nWordSize != 2 && nWordSize != 4)
     181                 :                 {
     182                 :                     CPLError(CE_Failure, CPLE_AppDefined,
     183               0 :                              "Unexpected word size : %d", (int)nWordSize);
     184               0 :                     break;
     185                 :                 }
     186                 : 
     187                 :                 GInt32 nVal;
     188            1528 :                 VSIFReadL(&nVal, 4, 1, poGDS->fp);
     189                 :                 CPL_LSBPTR32(&nVal);
     190            1528 :                 VSIFReadL(pabyData, nWordSize * (nTileWidth - 1), 1, poGDS->fp);
     191                 : #if defined(CPL_MSB)
     192                 :                 if (nWordSize > 1)
     193                 :                     GDALSwapWords(pabyData, nWordSize, nTileWidth - 1, nWordSize);
     194                 : #endif
     195                 : 
     196            1528 :                 pafBlockData[nxoff * nBlockXSize * nBlockXSize + j * nBlockXSize + 0] = nVal * fScale + fOff;
     197                 :                 int i;
     198          223368 :                 for(i=1;i<nTileWidth;i++)
     199                 :                 {
     200          221840 :                     if (nWordSize == 1)
     201           77840 :                         nVal += ((char*)pabyData)[i-1];
     202          144000 :                     else if (nWordSize == 2)
     203          144000 :                         nVal += ((GInt16*)pabyData)[i-1];
     204                 :                     else
     205               0 :                         nVal += ((GInt32*)pabyData)[i-1];
     206          221840 :                     pafBlockData[nxoff * nBlockXSize * nBlockXSize + j * nBlockXSize + i] = nVal * fScale + fOff;
     207                 :                 }
     208                 :             }
     209                 :         }
     210                 : 
     211              20 :         CPLFree(pabyData);
     212                 :     }
     213                 : 
     214            1528 :     int nTileWidth = MIN(nBlockXSize, nRasterXSize - nBlockXOff * nBlockXSize);
     215                 :     memcpy(pImage, pafBlockData + nBlockXOff * nBlockXSize * nBlockXSize +
     216                 :                                   nYOffInTile * nBlockXSize,
     217            1528 :            nTileWidth * sizeof(float));
     218                 : 
     219            1528 :     return CE_None;
     220                 : }
     221                 : 
     222                 : /************************************************************************/
     223                 : /*                            ~HF2Dataset()                            */
     224                 : /************************************************************************/
     225                 : 
     226              46 : HF2Dataset::HF2Dataset()
     227                 : {
     228              46 :     fp = NULL;
     229              46 :     pszWKT = NULL;
     230              46 :     panBlockOffset = NULL;
     231              46 :     adfGeoTransform[0] = 0;
     232              46 :     adfGeoTransform[1] = 1;
     233              46 :     adfGeoTransform[2] = 0;
     234              46 :     adfGeoTransform[3] = 0;
     235              46 :     adfGeoTransform[4] = 0;
     236              46 :     adfGeoTransform[5] = 1;
     237              46 :     bHasLoaderBlockMap = FALSE;
     238              46 :     nTileSize = 0;
     239              46 : }
     240                 : 
     241                 : /************************************************************************/
     242                 : /*                            ~HF2Dataset()                            */
     243                 : /************************************************************************/
     244                 : 
     245              46 : HF2Dataset::~HF2Dataset()
     246                 : 
     247                 : {
     248              46 :     FlushCache();
     249              46 :     CPLFree(pszWKT);
     250              46 :     CPLFree(panBlockOffset);
     251              46 :     if (fp)
     252              46 :         VSIFCloseL(fp);
     253              46 : }
     254                 : 
     255                 : /************************************************************************/
     256                 : /*                            LoadBlockMap()                            */
     257                 : /************************************************************************/
     258                 : 
     259            1528 : int HF2Dataset::LoadBlockMap()
     260                 : {
     261            1528 :     if (bHasLoaderBlockMap)
     262            1512 :         return panBlockOffset != NULL;
     263                 : 
     264              16 :     bHasLoaderBlockMap = TRUE;
     265                 : 
     266              16 :     int nXBlocks = (nRasterXSize + nTileSize - 1) / nTileSize;
     267              16 :     int nYBlocks = (nRasterYSize + nTileSize - 1) / nTileSize;
     268              16 :     panBlockOffset = (vsi_l_offset*) VSIMalloc3(sizeof(vsi_l_offset), nXBlocks, nYBlocks);
     269              16 :     if (panBlockOffset == NULL)
     270                 :     {
     271               0 :         return FALSE;
     272                 :     }
     273                 :     int i, j;
     274              36 :     for(j = 0; j < nYBlocks; j++)
     275                 :     {
     276              48 :         for(i = 0; i < nXBlocks; i++)
     277                 :         {
     278              28 :             vsi_l_offset nOff = VSIFTellL(fp);
     279              28 :             panBlockOffset[(nYBlocks - 1 - j) * nXBlocks + i] = nOff;
     280                 :             //VSIFSeekL(fp, 4 + 4, SEEK_CUR);
     281                 :             float fScale, fOff;
     282              28 :             VSIFReadL(&fScale, 4, 1, fp);
     283              28 :             VSIFReadL(&fOff, 4, 1, fp);
     284                 :             CPL_LSBPTR32(&fScale);
     285                 :             CPL_LSBPTR32(&fOff);
     286                 :             //printf("fScale = %f, fOff = %f\n", fScale, fOff);
     287                 :             int k;
     288              28 :             int nCols = MIN(nTileSize, nRasterXSize - nTileSize *i);
     289              28 :             int nLines = MIN(nTileSize, nRasterYSize - nTileSize *j);
     290            3112 :             for(k = 0; k < nLines; k++)
     291                 :             {
     292                 :                 GByte nWordSize;
     293            1528 :                 VSIFReadL(&nWordSize, 1, 1, fp);
     294                 :                 //printf("nWordSize=%d\n", nWordSize);
     295            1528 :                 if (nWordSize == 1 || nWordSize == 2 || nWordSize == 4)
     296            1528 :                     VSIFSeekL(fp, 4 + nWordSize * (nCols - 1), SEEK_CUR);
     297                 :                 else
     298                 :                 {
     299                 :                     CPLError(CE_Failure, CPLE_AppDefined,
     300                 :                             "Got unexpected byte depth (%d) for block (%d, %d) line %d",
     301               0 :                             (int)nWordSize, i, j, k);
     302               0 :                     VSIFree(panBlockOffset);
     303               0 :                     panBlockOffset = NULL;
     304               0 :                     return FALSE;
     305                 :                 }
     306                 :             }
     307                 :         }
     308                 :     }
     309                 : 
     310              16 :     return TRUE;
     311                 : }
     312                 : 
     313                 : /************************************************************************/
     314                 : /*                     GetProjectionRef()                               */
     315                 : /************************************************************************/
     316                 : 
     317               0 : const char* HF2Dataset::GetProjectionRef()
     318                 : {
     319               0 :     if (pszWKT)
     320               0 :         return pszWKT;
     321               0 :     return GDALPamDataset::GetProjectionRef();
     322                 : }
     323                 : 
     324                 : /************************************************************************/
     325                 : /*                             Identify()                               */
     326                 : /************************************************************************/
     327                 : 
     328           21868 : int HF2Dataset::Identify( GDALOpenInfo * poOpenInfo)
     329                 : {
     330                 : 
     331           21868 :     GDALOpenInfo* poOpenInfoToDelete = NULL;
     332                 :     /*  GZipped .hf2 files are common, so automagically open them */
     333                 :     /*  if the /vsigzip/ has not been explicitely passed */
     334           21868 :     CPLString osFilename(poOpenInfo->pszFilename);
     335           21868 :     if ((EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hfz") ||
     336                 :         (strlen(poOpenInfo->pszFilename) > 6 &&
     337                 :          EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6, "hf2.gz"))) &&
     338                 :          !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
     339                 :     {
     340               6 :         osFilename = "/vsigzip/";
     341               6 :         osFilename += poOpenInfo->pszFilename;
     342                 :         poOpenInfo = poOpenInfoToDelete =
     343                 :                 new GDALOpenInfo(osFilename.c_str(), GA_ReadOnly,
     344               6 :                                  poOpenInfo->papszSiblingFiles);
     345                 :     }
     346                 : 
     347           21868 :     if (poOpenInfo->nHeaderBytes < 28)
     348                 :     {
     349           21140 :         delete poOpenInfoToDelete;
     350           21140 :         return FALSE;
     351                 :     }
     352                 : 
     353             728 :     if (memcmp(poOpenInfo->pabyHeader, "HF2\0\0\0\0", 6) != 0)
     354                 :     {
     355             682 :         delete poOpenInfoToDelete;
     356             682 :         return FALSE;
     357                 :     }
     358                 : 
     359              46 :     delete poOpenInfoToDelete;
     360              46 :     return TRUE;
     361                 : }
     362                 : 
     363                 : /************************************************************************/
     364                 : /*                                Open()                                */
     365                 : /************************************************************************/
     366                 : 
     367            2974 : GDALDataset *HF2Dataset::Open( GDALOpenInfo * poOpenInfo )
     368                 : 
     369                 : {
     370            2974 :     CPLString osOriginalFilename(poOpenInfo->pszFilename);
     371                 : 
     372            2974 :     if (!Identify(poOpenInfo))
     373            2928 :         return NULL;
     374                 : 
     375              46 :     GDALOpenInfo* poOpenInfoToDelete = NULL;
     376                 :     /*  GZipped .hf2 files are common, so automagically open them */
     377                 :     /*  if the /vsigzip/ has not been explicitely passed */
     378              46 :     CPLString osFilename(poOpenInfo->pszFilename);
     379              46 :     if ((EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hfz") ||
     380                 :         (strlen(poOpenInfo->pszFilename) > 6 &&
     381                 :          EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6, "hf2.gz"))) &&
     382                 :          !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
     383                 :     {
     384               4 :         osFilename = "/vsigzip/";
     385               4 :         osFilename += poOpenInfo->pszFilename;
     386                 :         poOpenInfo = poOpenInfoToDelete =
     387                 :                 new GDALOpenInfo(osFilename.c_str(), GA_ReadOnly,
     388               4 :                                  poOpenInfo->papszSiblingFiles);
     389                 :     }
     390                 : 
     391                 : /* -------------------------------------------------------------------- */
     392                 : /*      Parse header                                                    */
     393                 : /* -------------------------------------------------------------------- */
     394                 : 
     395                 :     int nXSize, nYSize;
     396              46 :     memcpy(&nXSize, poOpenInfo->pabyHeader + 6, 4);
     397                 :     CPL_LSBPTR32(&nXSize);
     398              46 :     memcpy(&nYSize, poOpenInfo->pabyHeader + 10, 4);
     399                 :     CPL_LSBPTR32(&nYSize);
     400                 : 
     401                 :     GUInt16 nTileSize;
     402              46 :     memcpy(&nTileSize, poOpenInfo->pabyHeader + 14, 2);
     403                 :     CPL_LSBPTR16(&nTileSize);
     404                 : 
     405                 :     float fVertPres, fHorizScale;
     406              46 :     memcpy(&fVertPres, poOpenInfo->pabyHeader + 16, 4);
     407                 :     CPL_LSBPTR32(&fVertPres);
     408              46 :     memcpy(&fHorizScale, poOpenInfo->pabyHeader + 20, 4);
     409                 :     CPL_LSBPTR32(&fHorizScale);
     410                 : 
     411                 :     GUInt32 nExtendedHeaderLen;
     412              46 :     memcpy(&nExtendedHeaderLen, poOpenInfo->pabyHeader + 24, 4);
     413                 :     CPL_LSBPTR32(&nExtendedHeaderLen);
     414                 : 
     415              46 :     delete poOpenInfoToDelete;
     416              46 :     poOpenInfoToDelete = NULL;
     417                 : 
     418              46 :     if (nTileSize < 8)
     419               0 :         return NULL;
     420              46 :     if (nXSize <= 0 || nXSize > INT_MAX - nTileSize ||
     421                 :         nYSize <= 0 || nYSize > INT_MAX - nTileSize)
     422               0 :         return NULL;
     423                 :     /* To avoid later potential int overflows */
     424              46 :     if (nExtendedHeaderLen > 1024 * 65536)
     425               0 :         return NULL;
     426                 : 
     427              46 :     if (!GDALCheckDatasetDimensions(nXSize, nYSize))
     428                 :     {
     429               0 :         return NULL;
     430                 :     }
     431                 : 
     432                 : /* -------------------------------------------------------------------- */
     433                 : /*      Parse extended blocks                                           */
     434                 : /* -------------------------------------------------------------------- */
     435                 : 
     436              46 :     VSILFILE* fp = VSIFOpenL(osFilename.c_str(), "rb");
     437              46 :     if (fp == NULL)
     438               0 :         return NULL;
     439                 : 
     440              46 :     VSIFSeekL(fp, 28, SEEK_SET);
     441                 : 
     442              46 :     int bHasExtent = FALSE;
     443              46 :     double dfMinX = 0, dfMaxX = 0, dfMinY = 0, dfMaxY = 0;
     444              46 :     int bHasUTMZone = FALSE;
     445              46 :     GInt16 nUTMZone = 0;
     446              46 :     int bHasEPSGDatumCode = FALSE;
     447              46 :     GInt16 nEPSGDatumCode = 0;
     448              46 :     int bHasEPSGCode = FALSE;
     449              46 :     GInt16 nEPSGCode = 0;
     450              46 :     int bHasRelativePrecision = FALSE;
     451              46 :     float fRelativePrecision = 0;
     452                 :     char szApplicationName[256];
     453              46 :     szApplicationName[0] = 0;
     454                 : 
     455              46 :     GUInt32 nExtendedHeaderOff = 0;
     456             208 :     while(nExtendedHeaderOff < nExtendedHeaderLen)
     457                 :     {
     458                 :         char pabyBlockHeader[24];
     459             116 :         VSIFReadL(pabyBlockHeader, 24, 1, fp);
     460                 : 
     461                 :         char szBlockName[16 + 1];
     462             116 :         memcpy(szBlockName, pabyBlockHeader + 4, 16);
     463             116 :         szBlockName[16] = 0;
     464                 :         GUInt32 nBlockSize;
     465             116 :         memcpy(&nBlockSize, pabyBlockHeader + 20, 4);
     466                 :         CPL_LSBPTR32(&nBlockSize);
     467             116 :         if (nBlockSize > 65536)
     468               0 :             break;
     469                 : 
     470             116 :         nExtendedHeaderOff += 24 + nBlockSize;
     471                 : 
     472             162 :         if (strcmp(szBlockName, "georef-extents") == 0 &&
     473                 :             nBlockSize == 34)
     474                 :         {
     475                 :             char pabyBlockData[34];
     476              46 :             VSIFReadL(pabyBlockData, 34, 1, fp);
     477                 : 
     478              46 :             memcpy(&dfMinX, pabyBlockData + 2, 8);
     479                 :             CPL_LSBPTR64(&dfMinX);
     480              46 :             memcpy(&dfMaxX, pabyBlockData + 2 + 8, 8);
     481                 :             CPL_LSBPTR64(&dfMaxX);
     482              46 :             memcpy(&dfMinY, pabyBlockData + 2 + 8 + 8, 8);
     483                 :             CPL_LSBPTR64(&dfMinY);
     484              46 :             memcpy(&dfMaxY, pabyBlockData + 2 + 8 + 8 + 8, 8);
     485                 :             CPL_LSBPTR64(&dfMaxY);
     486                 : 
     487              46 :             bHasExtent = TRUE;
     488                 :         }
     489              88 :         else if (strcmp(szBlockName, "georef-utm") == 0 &&
     490                 :                 nBlockSize == 2)
     491                 :         {
     492              18 :             VSIFReadL(&nUTMZone, 2, 1, fp);
     493                 :             CPL_LSBPTR16(&nUTMZone);
     494              18 :             CPLDebug("HF2", "UTM Zone = %d", nUTMZone);
     495                 : 
     496              18 :             bHasUTMZone = TRUE;
     497                 :         }
     498              92 :         else if (strcmp(szBlockName, "georef-datum") == 0 &&
     499                 :                  nBlockSize == 2)
     500                 :         {
     501              40 :             VSIFReadL(&nEPSGDatumCode, 2, 1, fp);
     502                 :             CPL_LSBPTR16(&nEPSGDatumCode);
     503              40 :             CPLDebug("HF2", "EPSG Datum Code = %d", nEPSGDatumCode);
     504                 : 
     505              40 :             bHasEPSGDatumCode = TRUE;
     506                 :         }
     507              24 :         else if (strcmp(szBlockName, "georef-epsg-prj") == 0 &&
     508                 :                  nBlockSize == 2)
     509                 :         {
     510              12 :             VSIFReadL(&nEPSGCode, 2, 1, fp);
     511                 :             CPL_LSBPTR16(&nEPSGCode);
     512              12 :             CPLDebug("HF2", "EPSG Code = %d", nEPSGCode);
     513                 : 
     514              12 :             bHasEPSGCode = TRUE;
     515                 :         }
     516               0 :         else if (strcmp(szBlockName, "precis-rel") == 0 &&
     517                 :                  nBlockSize == 4)
     518                 :         {
     519               0 :             VSIFReadL(&fRelativePrecision, 4, 1, fp);
     520                 :             CPL_LSBPTR32(&fRelativePrecision);
     521                 : 
     522               0 :             bHasRelativePrecision = TRUE;
     523                 :         }
     524               0 :         else if (strcmp(szBlockName, "app-name") == 0 &&
     525                 :                  nBlockSize < 256)
     526                 :         {
     527               0 :             VSIFReadL(szApplicationName, nBlockSize, 1, fp);
     528               0 :             szApplicationName[nBlockSize] = 0;
     529                 :         }
     530                 :         else
     531                 :         {
     532               0 :             CPLDebug("HF2", "Skipping block %s", szBlockName);
     533               0 :             VSIFSeekL(fp, nBlockSize, SEEK_CUR);
     534                 :         }
     535                 :     }
     536                 : 
     537                 : /* -------------------------------------------------------------------- */
     538                 : /*      Create a corresponding GDALDataset.                             */
     539                 : /* -------------------------------------------------------------------- */
     540                 :     HF2Dataset         *poDS;
     541                 : 
     542              46 :     poDS = new HF2Dataset();
     543              46 :     poDS->fp = fp;
     544              46 :     poDS->nRasterXSize = nXSize;
     545              46 :     poDS->nRasterYSize = nYSize;
     546              46 :     poDS->nTileSize = nTileSize;
     547              46 :     CPLDebug("HF2", "nXSize = %d, nYSize = %d, nTileSize = %d", nXSize, nYSize, nTileSize);
     548              46 :     if (bHasExtent)
     549                 :     {
     550              46 :         poDS->adfGeoTransform[0] = dfMinX;
     551              46 :         poDS->adfGeoTransform[3] = dfMaxY;
     552              46 :         poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nXSize;
     553              46 :         poDS->adfGeoTransform[5] = -(dfMaxY - dfMinY) / nYSize;
     554                 :     }
     555                 :     else
     556                 :     {
     557               0 :         poDS->adfGeoTransform[1] = fHorizScale;
     558               0 :         poDS->adfGeoTransform[5] = fHorizScale;
     559                 :     }
     560                 : 
     561              46 :     if (bHasEPSGCode)
     562                 :     {
     563              12 :         OGRSpatialReference oSRS;
     564              12 :         if (oSRS.importFromEPSG(nEPSGCode) == OGRERR_NONE)
     565              12 :             oSRS.exportToWkt(&poDS->pszWKT);
     566                 :     }
     567                 :     else
     568                 :     {
     569              34 :         int bHasSRS = FALSE;
     570              34 :         OGRSpatialReference oSRS;
     571              34 :         oSRS.SetGeogCS("unknown", "unknown", "unknown", SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING); 
     572              34 :         if (bHasEPSGDatumCode)
     573                 :         {
     574              56 :             if (nEPSGDatumCode == 23 || nEPSGDatumCode == 6326)
     575                 :             {
     576              28 :                 bHasSRS = TRUE;
     577              28 :                 oSRS.SetWellKnownGeogCS("WGS84");
     578                 :             }
     579               0 :             else if (nEPSGDatumCode >= 6000)
     580                 :             {
     581                 :                 char szName[32];
     582               0 :                 sprintf( szName, "EPSG:%d", nEPSGDatumCode-2000 );
     583               0 :                 oSRS.SetWellKnownGeogCS( szName );
     584               0 :                 bHasSRS = TRUE;
     585                 :             }
     586                 :         }
     587                 : 
     588              34 :         if (bHasUTMZone && ABS(nUTMZone) >= 1 && ABS(nUTMZone) <= 60)
     589                 :         {
     590               6 :             bHasSRS = TRUE;
     591               6 :             oSRS.SetUTM(ABS(nUTMZone), nUTMZone > 0);
     592                 :         }
     593              34 :         if (bHasSRS)
     594              34 :             oSRS.exportToWkt(&poDS->pszWKT);
     595                 :     }
     596                 : 
     597                 : /* -------------------------------------------------------------------- */
     598                 : /*      Create band information objects.                                */
     599                 : /* -------------------------------------------------------------------- */
     600              46 :     poDS->nBands = 1;
     601                 :     int i;
     602              92 :     for( i = 0; i < poDS->nBands; i++ )
     603                 :     {
     604              46 :         poDS->SetBand( i+1, new HF2RasterBand( poDS, i+1, GDT_Float32 ) );
     605              46 :         poDS->GetRasterBand(i+1)->SetUnitType("m");
     606                 :     }
     607                 : 
     608              46 :     if (szApplicationName[0] != '\0')
     609               0 :         poDS->SetMetadataItem("APPLICATION_NAME", szApplicationName);
     610              46 :     poDS->SetMetadataItem("VERTICAL_PRECISION", CPLString().Printf("%f", fVertPres));
     611              46 :     if (bHasRelativePrecision)
     612               0 :         poDS->SetMetadataItem("RELATIVE_VERTICAL_PRECISION", CPLString().Printf("%f", fRelativePrecision));
     613                 : 
     614                 : /* -------------------------------------------------------------------- */
     615                 : /*      Initialize any PAM information.                                 */
     616                 : /* -------------------------------------------------------------------- */
     617              46 :     poDS->SetDescription( osOriginalFilename.c_str() );
     618              46 :     poDS->TryLoadXML();
     619                 : 
     620                 : /* -------------------------------------------------------------------- */
     621                 : /*      Support overviews.                                              */
     622                 : /* -------------------------------------------------------------------- */
     623              46 :     poDS->oOvManager.Initialize( poDS, osOriginalFilename.c_str() );
     624              46 :     return( poDS );
     625                 : }
     626                 : 
     627                 : /************************************************************************/
     628                 : /*                          GetGeoTransform()                           */
     629                 : /************************************************************************/
     630                 : 
     631               2 : CPLErr HF2Dataset::GetGeoTransform( double * padfTransform )
     632                 : 
     633                 : {
     634               2 :     memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
     635                 : 
     636               2 :     return( CE_None );
     637                 : }
     638                 : 
     639                 : 
     640           72128 : static void WriteShort(VSILFILE* fp, GInt16 val)
     641                 : {
     642                 :     CPL_LSBPTR16(&val);
     643           72128 :     VSIFWriteL(&val, 2, 1, fp);
     644           72128 : }
     645                 : 
     646            1142 : static void WriteInt(VSILFILE* fp, GInt32 val)
     647                 : {
     648                 :     CPL_LSBPTR32(&val);
     649            1142 :     VSIFWriteL(&val, 4, 1, fp);
     650            1142 : }
     651                 : 
     652             132 : static void WriteFloat(VSILFILE* fp, float val)
     653                 : {
     654                 :     CPL_LSBPTR32(&val);
     655             132 :     VSIFWriteL(&val, 4, 1, fp);
     656             132 : }
     657                 : 
     658             120 : static void WriteDouble(VSILFILE* fp, double val)
     659                 : {
     660                 :     CPL_LSBPTR64(&val);
     661             120 :     VSIFWriteL(&val, 8, 1, fp);
     662             120 : }
     663                 : 
     664                 : /************************************************************************/
     665                 : /*                             CreateCopy()                             */
     666                 : /************************************************************************/
     667                 : 
     668              44 : GDALDataset* HF2Dataset::CreateCopy( const char * pszFilename,
     669                 :                                      GDALDataset *poSrcDS, 
     670                 :                                      int bStrict, char ** papszOptions, 
     671                 :                                      GDALProgressFunc pfnProgress,
     672                 :                                      void * pProgressData )
     673                 : {
     674                 : /* -------------------------------------------------------------------- */
     675                 : /*      Some some rudimentary checks                                    */
     676                 : /* -------------------------------------------------------------------- */
     677              44 :     int nBands = poSrcDS->GetRasterCount();
     678              44 :     if (nBands == 0)
     679                 :     {
     680                 :         CPLError( CE_Failure, CPLE_NotSupported, 
     681               2 :                   "HF2 driver does not support source dataset with zero band.\n");
     682               2 :         return NULL;
     683                 :     }
     684                 : 
     685              42 :     if (nBands != 1)
     686                 :     {
     687                 :         CPLError( (bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported, 
     688               8 :                   "HF2 driver only uses the first band of the dataset.\n");
     689               8 :         if (bStrict)
     690               8 :             return NULL;
     691                 :     }
     692                 : 
     693              34 :     if( pfnProgress && !pfnProgress( 0.0, NULL, pProgressData ) )
     694               0 :         return NULL;
     695                 : 
     696                 : /* -------------------------------------------------------------------- */
     697                 : /*      Get source dataset info                                         */
     698                 : /* -------------------------------------------------------------------- */
     699                 : 
     700              34 :     int nXSize = poSrcDS->GetRasterXSize();
     701              34 :     int nYSize = poSrcDS->GetRasterYSize();
     702                 :     double adfGeoTransform[6];
     703              34 :     poSrcDS->GetGeoTransform(adfGeoTransform);
     704              34 :     int bHasGeoTransform = !(adfGeoTransform[0] == 0 &&
     705               0 :                              adfGeoTransform[1] == 1 &&
     706               0 :                              adfGeoTransform[2] == 0 &&
     707               0 :                              adfGeoTransform[3] == 0 &&
     708               0 :                              adfGeoTransform[4] == 0 &&
     709              34 :                              adfGeoTransform[5] == 1);
     710              34 :     if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
     711                 :     {
     712                 :         CPLError( CE_Failure, CPLE_NotSupported, 
     713               0 :                   "HF2 driver does not support CreateCopy() from skewed or rotated dataset.\n");
     714               0 :         return NULL;
     715                 :     }
     716                 : 
     717              34 :     GDALDataType eSrcDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
     718                 :     GDALDataType eReqDT;
     719              34 :     float fVertPres = (float) 0.01;
     720              48 :     if (eSrcDT == GDT_Byte || eSrcDT == GDT_Int16)
     721                 :     {
     722              14 :         fVertPres = 1;
     723              14 :         eReqDT = GDT_Int16;
     724                 :     }
     725                 :     else
     726              20 :         eReqDT = GDT_Float32;
     727                 : 
     728                 : /* -------------------------------------------------------------------- */
     729                 : /*      Read creation options                                           */
     730                 : /* -------------------------------------------------------------------- */
     731              34 :     const char* pszCompressed = CSLFetchNameValue(papszOptions, "COMPRESS");
     732              34 :     int bCompress = FALSE;
     733              34 :     if (pszCompressed)
     734               2 :         bCompress = CSLTestBoolean(pszCompressed);
     735                 :     
     736              34 :     const char* pszVerticalPrecision = CSLFetchNameValue(papszOptions, "VERTICAL_PRECISION");
     737              34 :     if (pszVerticalPrecision)
     738                 :     {
     739               0 :         fVertPres = (float) CPLAtofM(pszVerticalPrecision);
     740               0 :         if (fVertPres <= 0)
     741                 :         {
     742                 :             CPLError(CE_Warning, CPLE_AppDefined,
     743               0 :                      "Unsupported value for VERTICAL_PRECISION. Defaulting to 0.01");
     744               0 :             fVertPres = (float) 0.01;
     745                 :         }
     746               0 :         if (eReqDT == GDT_Int16 && fVertPres > 1)
     747               0 :             eReqDT = GDT_Float32;
     748                 :     }
     749                 : 
     750              34 :     const char* pszBlockSize = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
     751              34 :     int nTileSize = 256;
     752              34 :     if (pszBlockSize)
     753                 :     {
     754               2 :         nTileSize = atoi(pszBlockSize);
     755               2 :         if (nTileSize < 8 || nTileSize > 4096)
     756                 :         {
     757                 :             CPLError(CE_Warning, CPLE_AppDefined,
     758               0 :                      "Unsupported value for BLOCKSIZE. Defaulting to 256");
     759               0 :             nTileSize = 256;
     760                 :         }
     761                 :     }
     762                 : 
     763                 : /* -------------------------------------------------------------------- */
     764                 : /*      Parse source dataset georeferencing info                        */
     765                 : /* -------------------------------------------------------------------- */
     766                 : 
     767              34 :     int nExtendedHeaderLen = 0;
     768              34 :     if (bHasGeoTransform)
     769              34 :         nExtendedHeaderLen += 58;
     770              34 :     const char* pszProjectionRef = poSrcDS->GetProjectionRef();
     771              34 :     int nDatumCode = -2;
     772              34 :     int nUTMZone = 0;
     773              34 :     int bNorth = FALSE;
     774              34 :     int nEPSGCode = 0;
     775              34 :     int nExtentUnits = 1;
     776              34 :     if (pszProjectionRef != NULL && pszProjectionRef[0] != '\0')
     777                 :     {
     778              34 :         OGRSpatialReference oSRS;
     779              34 :         char* pszTemp = (char*) pszProjectionRef;
     780              34 :         if (oSRS.importFromWkt(&pszTemp) == OGRERR_NONE)
     781                 :         {
     782              34 :             const char* pszValue = NULL;
     783              34 :             if( oSRS.GetAuthorityName( "GEOGCS|DATUM" ) != NULL
     784                 :                 && EQUAL(oSRS.GetAuthorityName( "GEOGCS|DATUM" ),"EPSG") )
     785              26 :                 nDatumCode = atoi(oSRS.GetAuthorityCode( "GEOGCS|DATUM" ));
     786               8 :             else if ((pszValue = oSRS.GetAttrValue("GEOGCS|DATUM")) != NULL)
     787                 :             {
     788               8 :                 if (strstr(pszValue, "WGS") && strstr(pszValue, "84"))
     789               6 :                     nDatumCode = 6326;
     790                 :             }
     791                 : 
     792              34 :             nUTMZone = oSRS.GetUTMZone(&bNorth);
     793                 :         }
     794              34 :         if( oSRS.GetAuthorityName( "PROJCS" ) != NULL
     795                 :             && EQUAL(oSRS.GetAuthorityName( "PROJCS" ),"EPSG") )
     796               4 :             nEPSGCode = atoi(oSRS.GetAuthorityCode( "PROJCS" ));
     797                 : 
     798              34 :         if( oSRS.IsGeographic() )
     799                 :         {
     800              28 :             nExtentUnits = 0;
     801                 :         }
     802                 :         else
     803                 :         {
     804               6 :             double dfLinear = oSRS.GetLinearUnits();
     805                 : 
     806               6 :             if( ABS(dfLinear - 0.3048) < 0.0000001 )
     807               0 :                 nExtentUnits = 2;
     808               6 :             else if( ABS(dfLinear - atof(SRS_UL_US_FOOT_CONV)) < 0.00000001 )
     809               0 :                 nExtentUnits = 3;
     810                 :             else
     811               6 :                 nExtentUnits = 1;
     812              34 :         }
     813                 :     }
     814              34 :     if (nDatumCode != -2)
     815              32 :         nExtendedHeaderLen += 26;
     816              34 :     if (nUTMZone != 0)
     817               6 :         nExtendedHeaderLen += 26;
     818              34 :     if (nEPSGCode)
     819               4 :         nExtendedHeaderLen += 26;
     820                 : 
     821                 : /* -------------------------------------------------------------------- */
     822                 : /*      Create target file                                              */
     823                 : /* -------------------------------------------------------------------- */
     824                 : 
     825              34 :     CPLString osFilename;
     826              34 :     if (bCompress)
     827                 :     {
     828               2 :         osFilename = "/vsigzip/";
     829               2 :         osFilename += pszFilename;
     830                 :     }
     831                 :     else
     832              32 :         osFilename = pszFilename;
     833              34 :     VSILFILE* fp = VSIFOpenL(osFilename.c_str(), "wb");
     834              34 :     if (fp == NULL)
     835                 :     {
     836                 :         CPLError( CE_Failure, CPLE_AppDefined, 
     837               4 :                   "Cannot create %s", pszFilename );
     838               4 :         return NULL;
     839                 :     }
     840                 : 
     841                 : /* -------------------------------------------------------------------- */
     842                 : /*      Write header                                                    */
     843                 : /* -------------------------------------------------------------------- */
     844                 : 
     845              30 :     VSIFWriteL("HF2\0", 4, 1, fp);
     846              30 :     WriteShort(fp, 0);
     847              30 :     WriteInt(fp, nXSize);
     848              30 :     WriteInt(fp, nYSize);
     849              30 :     WriteShort(fp, (GInt16) nTileSize);
     850              30 :     WriteFloat(fp, fVertPres);
     851              30 :     float fHorizScale = (float) ((fabs(adfGeoTransform[1]) + fabs(adfGeoTransform[5])) / 2);
     852              30 :     WriteFloat(fp, fHorizScale);
     853              30 :     WriteInt(fp, nExtendedHeaderLen);
     854                 : 
     855                 : /* -------------------------------------------------------------------- */
     856                 : /*      Write extended header                                           */
     857                 : /* -------------------------------------------------------------------- */
     858                 : 
     859                 :     char szBlockName[16 + 1];
     860              30 :     if (bHasGeoTransform)
     861                 :     {
     862              30 :         VSIFWriteL("bin\0", 4, 1, fp);
     863              30 :         memset(szBlockName, 0, 16 + 1);
     864              30 :         strcpy(szBlockName, "georef-extents");
     865              30 :         VSIFWriteL(szBlockName, 16, 1, fp);
     866              30 :         WriteInt(fp, 34);
     867              30 :         WriteShort(fp, (GInt16) nExtentUnits);
     868              30 :         WriteDouble(fp, adfGeoTransform[0]);
     869              30 :         WriteDouble(fp, adfGeoTransform[0] + nXSize * adfGeoTransform[1]);
     870              30 :         WriteDouble(fp, adfGeoTransform[3] + nYSize * adfGeoTransform[5]);
     871              30 :         WriteDouble(fp, adfGeoTransform[3]);
     872                 :     }
     873              30 :     if (nUTMZone != 0)
     874                 :     {
     875               6 :         VSIFWriteL("bin\0", 4, 1, fp);
     876               6 :         memset(szBlockName, 0, 16 + 1);
     877               6 :         strcpy(szBlockName, "georef-utm");
     878               6 :         VSIFWriteL(szBlockName, 16, 1, fp);
     879               6 :         WriteInt(fp, 2);
     880               6 :         WriteShort(fp, (GInt16) ((bNorth) ? nUTMZone : -nUTMZone));
     881                 :     }
     882              30 :     if (nDatumCode != -2)
     883                 :     {
     884              28 :         VSIFWriteL("bin\0", 4, 1, fp);
     885              28 :         memset(szBlockName, 0, 16 + 1);
     886              28 :         strcpy(szBlockName, "georef-datum");
     887              28 :         VSIFWriteL(szBlockName, 16, 1, fp);
     888              28 :         WriteInt(fp, 2);
     889              28 :         WriteShort(fp, (GInt16) nDatumCode);
     890                 :     }
     891              30 :     if (nEPSGCode != 0)
     892                 :     {
     893               4 :         VSIFWriteL("bin\0", 4, 1, fp);
     894               4 :         memset(szBlockName, 0, 16 + 1);
     895               4 :         strcpy(szBlockName, "georef-epsg-prj");
     896               4 :         VSIFWriteL(szBlockName, 16, 1, fp);
     897               4 :         WriteInt(fp, 2);
     898               4 :         WriteShort(fp, (GInt16) nEPSGCode);
     899                 :     }
     900                 : 
     901                 : /* -------------------------------------------------------------------- */
     902                 : /*      Copy imagery                                                    */
     903                 : /* -------------------------------------------------------------------- */
     904              30 :     int nXBlocks = (nXSize + nTileSize - 1) / nTileSize;
     905              30 :     int nYBlocks = (nYSize + nTileSize - 1) / nTileSize;
     906                 : 
     907              30 :     void* pTileBuffer = (void*) VSIMalloc(nTileSize * nTileSize * (GDALGetDataTypeSize(eReqDT) / 8));
     908              30 :     if (pTileBuffer == NULL)
     909                 :     {
     910               0 :         CPLError( CE_Failure, CPLE_OutOfMemory, "Out of memory");
     911               0 :         VSIFCloseL(fp);
     912               0 :         return NULL;
     913                 :     }
     914                 : 
     915                 :     int i, j, k, l;
     916              30 :     CPLErr eErr = CE_None;
     917              62 :     for(j=0;j<nYBlocks && eErr == CE_None;j++)
     918                 :     {
     919              68 :         for(i=0;i<nXBlocks && eErr == CE_None;i++)
     920                 :         {
     921              36 :             int nReqXSize = MIN(nTileSize, nXSize - i * nTileSize);
     922              36 :             int nReqYSize = MIN(nTileSize, nYSize - j * nTileSize);
     923                 :             eErr = poSrcDS->GetRasterBand(1)->RasterIO(GF_Read,
     924                 :                                                 i * nTileSize, MAX(0, nYSize - (j + 1) * nTileSize),
     925                 :                                                 nReqXSize, nReqYSize,
     926                 :                                                 pTileBuffer, nReqXSize, nReqYSize,
     927              36 :                                                 eReqDT, 0, 0);
     928              36 :             if (eErr != CE_None)
     929               0 :                 break;
     930                 : 
     931              36 :             if (eReqDT == GDT_Int16)
     932                 :             {
     933              16 :                 WriteFloat(fp, 1); /* scale */
     934              16 :                 WriteFloat(fp, 0); /* offset */
     935             418 :                 for(k=0;k<nReqYSize;k++)
     936                 :                 {
     937             402 :                     int nLastVal = ((short*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
     938             402 :                     GByte nWordSize = 1;
     939           31282 :                     for(l=1;l<nReqXSize;l++)
     940                 :                     {
     941           30880 :                         int nVal = ((short*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + l];
     942           30880 :                         int nDiff = nVal - nLastVal;
     943           30880 :                         if (nDiff < -32768 || nDiff > 32767)
     944                 :                         {
     945               0 :                             nWordSize = 4;
     946               0 :                             break;
     947                 :                         }
     948           30880 :                         if (nDiff < -128 || nDiff > 127)
     949               0 :                             nWordSize = 2;
     950           30880 :                         nLastVal = nVal;
     951                 :                     }
     952                 : 
     953             402 :                     VSIFWriteL(&nWordSize, 1, 1, fp);
     954             402 :                     nLastVal = ((short*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
     955             402 :                     WriteInt(fp, nLastVal);
     956           31282 :                     for(l=1;l<nReqXSize;l++)
     957                 :                     {
     958           30880 :                         int nVal = ((short*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + l];
     959           30880 :                         int nDiff = nVal - nLastVal;
     960           30880 :                         if (nWordSize == 1)
     961                 :                         {
     962           30880 :                             CPLAssert(nDiff >= -128 && nDiff <= 127);
     963           30880 :                             char chDiff = (char)nDiff;
     964           30880 :                             VSIFWriteL(&chDiff, 1, 1, fp);
     965                 :                         }
     966               0 :                         else if (nWordSize == 2)
     967                 :                         {
     968               0 :                             CPLAssert(nDiff >= -32768 && nDiff <= 32767);
     969               0 :                             WriteShort(fp, (short)nDiff);
     970                 :                         }
     971                 :                         else
     972                 :                         {
     973               0 :                             WriteInt(fp, nDiff);
     974                 :                         }
     975           30880 :                         nLastVal = nVal;
     976                 :                     }
     977                 :                 }
     978                 :             }
     979                 :             else
     980                 :             {
     981              20 :                 float fMinVal = ((float*)pTileBuffer)[0];
     982              20 :                 float fMaxVal = fMinVal;
     983           82602 :                 for(k=1;k<nReqYSize*nReqXSize;k++)
     984                 :                 {
     985           82582 :                     float fVal = ((float*)pTileBuffer)[k];
     986           82582 :                     if (fVal < fMinVal) fMinVal = fVal;
     987           82582 :                     if (fVal > fMaxVal) fMaxVal = fVal;
     988                 :                 }
     989                 : 
     990              20 :                 float fIntRange = (fMaxVal - fMinVal) / fVertPres;
     991              20 :                 float fScale = (fMinVal == fMaxVal) ? 1 : (fMaxVal - fMinVal) / fIntRange;
     992              20 :                 float fOffset = fMinVal;
     993              20 :                 WriteFloat(fp, fScale); /* scale */
     994              20 :                 WriteFloat(fp, fOffset); /* offset */
     995             602 :                 for(k=0;k<nReqYSize;k++)
     996                 :                 {
     997             582 :                     float fLastVal = ((float*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
     998             582 :                     float fIntLastVal = (fLastVal - fOffset) / fScale;
     999             582 :                     CPLAssert(fIntLastVal >= -2147483648.0f && fIntLastVal <= 2147483647.0f);
    1000             582 :                     int nLastVal = (int)fIntLastVal;
    1001             582 :                     GByte nWordSize = 1;
    1002           82602 :                     for(l=1;l<nReqXSize;l++)
    1003                 :                     {
    1004           82020 :                         float fVal = ((float*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + l];
    1005           82020 :                         float fIntVal = (fVal - fOffset) / fScale;
    1006           82020 :                         CPLAssert(fIntVal >= -2147483648.0f && fIntVal <= 2147483647.0f);
    1007           82020 :                         int nVal = (int)fIntVal;
    1008           82020 :                         int nDiff = nVal - nLastVal;
    1009           82020 :                         CPLAssert((int)((GIntBig)nVal - nLastVal) == nDiff);
    1010           82020 :                         if (nDiff < -32768 || nDiff > 32767)
    1011                 :                         {
    1012               0 :                             nWordSize = 4;
    1013               0 :                             break;
    1014                 :                         }
    1015           82020 :                         if (nDiff < -128 || nDiff > 127)
    1016             402 :                             nWordSize = 2;
    1017           82020 :                         nLastVal = nVal;
    1018                 :                     }
    1019                 : 
    1020             582 :                     VSIFWriteL(&nWordSize, 1, 1, fp);
    1021             582 :                     fLastVal = ((float*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
    1022             582 :                     fIntLastVal = (fLastVal - fOffset) / fScale;
    1023             582 :                     nLastVal = (int)fIntLastVal;
    1024             582 :                     WriteInt(fp, nLastVal);
    1025           82602 :                     for(l=1;l<nReqXSize;l++)
    1026                 :                     {
    1027           82020 :                         float fVal = ((float*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + l];
    1028           82020 :                         float fIntVal = (fVal - fOffset) / fScale;
    1029           82020 :                         int nVal = (int)fIntVal;
    1030           82020 :                         int nDiff = nVal - nLastVal;
    1031           82020 :                         CPLAssert((int)((GIntBig)nVal - nLastVal) == nDiff);
    1032           82020 :                         if (nWordSize == 1)
    1033                 :                         {
    1034           10020 :                             CPLAssert(nDiff >= -128 && nDiff <= 127);
    1035           10020 :                             char chDiff = (char)nDiff;
    1036           10020 :                             VSIFWriteL(&chDiff, 1, 1, fp);
    1037                 :                         }
    1038           72000 :                         else if (nWordSize == 2)
    1039                 :                         {
    1040           72000 :                             CPLAssert(nDiff >= -32768 && nDiff <= 32767);
    1041           72000 :                             WriteShort(fp, (short)nDiff);
    1042                 :                         }
    1043                 :                         else
    1044                 :                         {
    1045               0 :                             WriteInt(fp, nDiff);
    1046                 :                         }
    1047           82020 :                         nLastVal = nVal;
    1048                 :                     }
    1049                 :                 }
    1050                 :             }
    1051                 : 
    1052              36 :             if( pfnProgress && !pfnProgress( (j * nXBlocks + i + 1) * 1.0 / (nXBlocks * nYBlocks), NULL, pProgressData ) )
    1053                 :             {
    1054               0 :                 eErr = CE_Failure;
    1055               0 :                 break;
    1056                 :             }
    1057                 :         }
    1058                 :     }
    1059                 : 
    1060              30 :     CPLFree(pTileBuffer);
    1061                 : 
    1062              30 :     VSIFCloseL(fp);
    1063                 : 
    1064              30 :     return (GDALDataset*) GDALOpen(osFilename.c_str(), GA_ReadOnly);
    1065                 : }
    1066                 : 
    1067                 : /************************************************************************/
    1068                 : /*                         GDALRegister_HF2()                           */
    1069                 : /************************************************************************/
    1070                 : 
    1071            1135 : void GDALRegister_HF2()
    1072                 : 
    1073                 : {
    1074                 :     GDALDriver  *poDriver;
    1075                 : 
    1076            1135 :     if( GDALGetDriverByName( "HF2" ) == NULL )
    1077                 :     {
    1078            1093 :         poDriver = new GDALDriver();
    1079                 :         
    1080            1093 :         poDriver->SetDescription( "HF2" );
    1081                 :         poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
    1082            1093 :                                    "HF2/HFZ heightfield raster" );
    1083                 :         poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, 
    1084            1093 :                                    "frmt_hf2.html" );
    1085            1093 :         poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "hf2" );
    1086                 :         poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
    1087                 : "<CreationOptionList>"
    1088                 : "   <Option name='VERTICAL_PRECISION' type='float' default='0.01' description='Vertical precision.'/>"
    1089                 : "   <Option name='COMPRESS' type='boolean' default='false' description='Set to true to produce a GZip compressed file.'/>"
    1090                 : "   <Option name='BLOCKSIZE' type='int' default='256' description='Tile size.'/>"
    1091            1093 : "</CreationOptionList>");
    1092                 : 
    1093            1093 :         poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
    1094                 : 
    1095            1093 :         poDriver->pfnOpen = HF2Dataset::Open;
    1096            1093 :         poDriver->pfnIdentify = HF2Dataset::Identify;
    1097            1093 :         poDriver->pfnCreateCopy = HF2Dataset::CreateCopy;
    1098                 : 
    1099            1093 :         GetGDALDriverManager()->RegisterDriver( poDriver );
    1100                 :     }
    1101            1135 : }
    1102                 : 

Generated by: LCOV version 1.7