Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  Implementation for ELAS DIPEx format variant.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2006, Frank Warmerdam
       9             :  * Copyright (c) 2008-2011, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_string.h"
      15             : #include "gdal_frmts.h"
      16             : #include "ogr_spatialref.h"
      17             : #include "rawdataset.h"
      18             : 
      19             : #include <cmath>
      20             : #include <algorithm>
      21             : 
      22             : using std::fill;
      23             : 
      24             : /************************************************************************/
      25             : /* ==================================================================== */
      26             : /*                              DIPExDataset                            */
      27             : /* ==================================================================== */
      28             : /************************************************************************/
      29             : 
      30             : class DIPExDataset final : public GDALPamDataset
      31             : {
      32             :     struct DIPExHeader
      33             :     {
      34             :         GInt32 NBIH = {0};  /* bytes in header, normally 1024 */
      35             :         GInt32 NBPR = {0};  /* bytes per data record (all bands of scanline) */
      36             :         GInt32 IL = {0};    /* initial line - normally 1 */
      37             :         GInt32 LL = {0};    /* last line */
      38             :         GInt32 IE = {0};    /* initial element (pixel), normally 1 */
      39             :         GInt32 LE = {0};    /* last element (pixel) */
      40             :         GInt32 NC = {0};    /* number of channels (bands) */
      41             :         GInt32 H4322 = {0}; /* header record identifier - always 4322. */
      42             :         char unused1[40] = {0};
      43             :         GByte IH19[4] = {0}; /* data type, and size flags */
      44             :         GInt32 IH20 = {0};   /* number of secondary headers */
      45             :         GInt32 SRID = {0};
      46             :         char unused2[12] = {0};
      47             :         double YOffset = {0};
      48             :         double XOffset = {0};
      49             :         double YPixSize = {0};
      50             :         double XPixSize = {0};
      51             :         double Matrix[4] = {0};
      52             :         char unused3[344] = {0};
      53             :         GUInt16 ColorTable[256] = {0}; /* RGB packed with 4 bits each */
      54             :         char unused4[32] = {0};
      55             :     };
      56             : 
      57             :     VSILFILE *fp;
      58             :     OGRSpatialReference m_oSRS{};
      59             : 
      60             :     DIPExHeader sHeader{};
      61             : 
      62             :     GDALDataType eRasterDataType;
      63             : 
      64             :     double adfGeoTransform[6];
      65             : 
      66             :     CPL_DISALLOW_COPY_ASSIGN(DIPExDataset)
      67             : 
      68             :   public:
      69             :     DIPExDataset();
      70             :     ~DIPExDataset() override;
      71             : 
      72             :     CPLErr GetGeoTransform(double *) override;
      73             : 
      74           0 :     const OGRSpatialReference *GetSpatialRef() const override
      75             :     {
      76           0 :         return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
      77             :     }
      78             : 
      79             :     static GDALDataset *Open(GDALOpenInfo *);
      80             : };
      81             : 
      82             : /************************************************************************/
      83             : /* ==================================================================== */
      84             : /*                             DIPExDataset                              */
      85             : /* ==================================================================== */
      86             : /************************************************************************/
      87             : 
      88             : /************************************************************************/
      89             : /*                            DIPExDataset()                             */
      90             : /************************************************************************/
      91             : 
      92           2 : DIPExDataset::DIPExDataset() : fp(nullptr), eRasterDataType(GDT_Unknown)
      93             : {
      94           2 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
      95           2 :     adfGeoTransform[0] = 0.0;
      96           2 :     adfGeoTransform[1] = 1.0;
      97           2 :     adfGeoTransform[2] = 0.0;
      98           2 :     adfGeoTransform[3] = 0.0;
      99           2 :     adfGeoTransform[4] = 0.0;
     100           2 :     adfGeoTransform[5] = 1.0;
     101           2 : }
     102             : 
     103             : /************************************************************************/
     104             : /*                            ~DIPExDataset()                            */
     105             : /************************************************************************/
     106             : 
     107           4 : DIPExDataset::~DIPExDataset()
     108             : 
     109             : {
     110           2 :     if (fp)
     111           2 :         CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
     112           2 :     fp = nullptr;
     113           4 : }
     114             : 
     115             : /************************************************************************/
     116             : /*                                Open()                                */
     117             : /************************************************************************/
     118             : 
     119       30562 : GDALDataset *DIPExDataset::Open(GDALOpenInfo *poOpenInfo)
     120             : 
     121             : {
     122             :     /* -------------------------------------------------------------------- */
     123             :     /*      First we check to see if the file has the expected header       */
     124             :     /*      bytes.                                                          */
     125             :     /* -------------------------------------------------------------------- */
     126       30562 :     if (poOpenInfo->nHeaderBytes < 256 || poOpenInfo->fpL == nullptr)
     127       27838 :         return nullptr;
     128             : 
     129        2724 :     if (CPL_LSBWORD32(
     130             :             *(reinterpret_cast<GInt32 *>(poOpenInfo->pabyHeader + 0))) != 1024)
     131        2722 :         return nullptr;
     132             : 
     133           2 :     if (CPL_LSBWORD32(
     134             :             *(reinterpret_cast<GInt32 *>(poOpenInfo->pabyHeader + 28))) != 4322)
     135           0 :         return nullptr;
     136             : 
     137             :     /* -------------------------------------------------------------------- */
     138             :     /*      Create a corresponding GDALDataset.                             */
     139             :     /* -------------------------------------------------------------------- */
     140           4 :     auto poDS = std::make_unique<DIPExDataset>();
     141             : 
     142           2 :     poDS->eAccess = poOpenInfo->eAccess;
     143           2 :     poDS->fp = poOpenInfo->fpL;
     144           2 :     poOpenInfo->fpL = nullptr;
     145             : 
     146             :     /* -------------------------------------------------------------------- */
     147             :     /*      Read the header information.                                    */
     148             :     /* -------------------------------------------------------------------- */
     149           2 :     if (VSIFReadL(&(poDS->sHeader), 1024, 1, poDS->fp) != 1)
     150             :     {
     151           0 :         CPLError(CE_Failure, CPLE_FileIO,
     152             :                  "Attempt to read 1024 byte header filed on file %s\n",
     153             :                  poOpenInfo->pszFilename);
     154           0 :         return nullptr;
     155             :     }
     156             : 
     157             :     // To avoid cppcheck warnings about unused members
     158           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.NBIH);
     159           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.H4322);
     160           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.unused1);
     161           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.IH20);
     162           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.unused2);
     163           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.Matrix);
     164           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.unused3);
     165           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.ColorTable);
     166           2 :     CPL_IGNORE_RET_VAL(poDS->sHeader.unused4);
     167             : 
     168             :     /* -------------------------------------------------------------------- */
     169             :     /*      Extract information of interest from the header.                */
     170             :     /* -------------------------------------------------------------------- */
     171           2 :     const int nLineOffset = CPL_LSBWORD32(poDS->sHeader.NBPR);
     172             : 
     173           2 :     int nStart = CPL_LSBWORD32(poDS->sHeader.IL);
     174           2 :     int nEnd = CPL_LSBWORD32(poDS->sHeader.LL);
     175           2 :     GIntBig nDiff = static_cast<GIntBig>(nEnd) - nStart + 1;
     176           2 :     if (nDiff <= 0 || nDiff > INT_MAX)
     177             :     {
     178           0 :         return nullptr;
     179             :     }
     180           2 :     poDS->nRasterYSize = static_cast<int>(nDiff);
     181             : 
     182           2 :     nStart = CPL_LSBWORD32(poDS->sHeader.IE);
     183           2 :     nEnd = CPL_LSBWORD32(poDS->sHeader.LE);
     184           2 :     nDiff = static_cast<GIntBig>(nEnd) - nStart + 1;
     185           2 :     if (nDiff <= 0 || nDiff > INT_MAX)
     186             :     {
     187           0 :         return nullptr;
     188             :     }
     189           2 :     poDS->nRasterXSize = static_cast<int>(nDiff);
     190             : 
     191           2 :     const int nBands = CPL_LSBWORD32(poDS->sHeader.NC);
     192             : 
     193           4 :     if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
     194           2 :         !GDALCheckBandCount(nBands, FALSE))
     195             :     {
     196           0 :         return nullptr;
     197             :     }
     198             : 
     199           2 :     const int nDIPExDataType = (poDS->sHeader.IH19[1] & 0x7e) >> 2;
     200           2 :     const int nBytesPerSample = poDS->sHeader.IH19[0];
     201             : 
     202           2 :     if (nDIPExDataType == 0 && nBytesPerSample == 1)
     203           2 :         poDS->eRasterDataType = GDT_Byte;
     204           0 :     else if (nDIPExDataType == 1 && nBytesPerSample == 1)
     205           0 :         poDS->eRasterDataType = GDT_Byte;
     206           0 :     else if (nDIPExDataType == 16 && nBytesPerSample == 4)
     207           0 :         poDS->eRasterDataType = GDT_Float32;
     208           0 :     else if (nDIPExDataType == 17 && nBytesPerSample == 8)
     209           0 :         poDS->eRasterDataType = GDT_Float64;
     210             :     else
     211             :     {
     212           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     213             :                  "Unrecognized image data type %d, with BytesPerSample=%d.",
     214             :                  nDIPExDataType, nBytesPerSample);
     215           0 :         return nullptr;
     216             :     }
     217             : 
     218           2 :     if (nLineOffset <= 0 || nLineOffset > INT_MAX / nBands)
     219             :     {
     220           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     221             :                  "Invalid values: nLineOffset = %d, nBands = %d.", nLineOffset,
     222             :                  nBands);
     223           0 :         return nullptr;
     224             :     }
     225             : 
     226             :     /* -------------------------------------------------------------------- */
     227             :     /*      Create band information objects.                                */
     228             :     /* -------------------------------------------------------------------- */
     229           2 :     CPLErrorReset();
     230           4 :     for (int iBand = 0; iBand < nBands; iBand++)
     231             :     {
     232             :         auto poBand = RawRasterBand::Create(
     233           4 :             poDS.get(), iBand + 1, poDS->fp, 1024 + iBand * nLineOffset,
     234           2 :             nBytesPerSample, nLineOffset * nBands, poDS->eRasterDataType,
     235             :             RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
     236           2 :             RawRasterBand::OwnFP::NO);
     237           2 :         if (!poBand)
     238           0 :             return nullptr;
     239           2 :         poDS->SetBand(iBand + 1, std::move(poBand));
     240             :     }
     241             : 
     242             :     /* -------------------------------------------------------------------- */
     243             :     /*      Extract the projection coordinates, if present.                 */
     244             :     /* -------------------------------------------------------------------- */
     245           2 :     CPL_LSBPTR64(&(poDS->sHeader.XPixSize));
     246           2 :     CPL_LSBPTR64(&(poDS->sHeader.YPixSize));
     247           2 :     CPL_LSBPTR64(&(poDS->sHeader.XOffset));
     248           2 :     CPL_LSBPTR64(&(poDS->sHeader.YOffset));
     249             : 
     250           2 :     if (poDS->sHeader.XOffset != 0)
     251             :     {
     252           0 :         poDS->adfGeoTransform[0] = poDS->sHeader.XOffset;
     253           0 :         poDS->adfGeoTransform[1] = poDS->sHeader.XPixSize;
     254           0 :         poDS->adfGeoTransform[2] = 0.0;
     255           0 :         poDS->adfGeoTransform[3] = poDS->sHeader.YOffset;
     256           0 :         poDS->adfGeoTransform[4] = 0.0;
     257           0 :         poDS->adfGeoTransform[5] = -1.0 * std::abs(poDS->sHeader.YPixSize);
     258             : 
     259           0 :         poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[1] * 0.5;
     260           0 :         poDS->adfGeoTransform[3] -= poDS->adfGeoTransform[5] * 0.5;
     261             :     }
     262             :     else
     263             :     {
     264           2 :         poDS->adfGeoTransform[0] = 0.0;
     265           2 :         poDS->adfGeoTransform[1] = 1.0;
     266           2 :         poDS->adfGeoTransform[2] = 0.0;
     267           2 :         poDS->adfGeoTransform[3] = 0.0;
     268           2 :         poDS->adfGeoTransform[4] = 0.0;
     269           2 :         poDS->adfGeoTransform[5] = 1.0;
     270             :     }
     271             : 
     272             :     /* -------------------------------------------------------------------- */
     273             :     /*      Look for SRID.                                                  */
     274             :     /* -------------------------------------------------------------------- */
     275           2 :     CPL_LSBPTR32(&(poDS->sHeader.SRID));
     276             : 
     277           2 :     if (poDS->sHeader.SRID > 0 && poDS->sHeader.SRID < 33000)
     278             :     {
     279           4 :         OGRSpatialReference oSR;
     280           2 :         oSR.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     281           2 :         if (oSR.importFromEPSG(poDS->sHeader.SRID) == OGRERR_NONE)
     282             :         {
     283           2 :             poDS->m_oSRS = std::move(oSR);
     284             :         }
     285             :     }
     286             : 
     287             :     /* -------------------------------------------------------------------- */
     288             :     /*      Initialize any PAM information.                                 */
     289             :     /* -------------------------------------------------------------------- */
     290           2 :     poDS->SetDescription(poOpenInfo->pszFilename);
     291           2 :     poDS->TryLoadXML();
     292             : 
     293             :     /* -------------------------------------------------------------------- */
     294             :     /*      Check for external overviews.                                   */
     295             :     /* -------------------------------------------------------------------- */
     296           4 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename,
     297           2 :                                 poOpenInfo->GetSiblingFiles());
     298             : 
     299           2 :     return poDS.release();
     300             : }
     301             : 
     302             : /************************************************************************/
     303             : /*                          GetGeoTransform()                           */
     304             : /************************************************************************/
     305             : 
     306           0 : CPLErr DIPExDataset::GetGeoTransform(double *padfTransform)
     307             : 
     308             : {
     309           0 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     310             : 
     311           0 :     return CE_None;
     312             : }
     313             : 
     314             : /************************************************************************/
     315             : /*                          GDALRegister_DIPEx()                        */
     316             : /************************************************************************/
     317             : 
     318        1595 : void GDALRegister_DIPEx()
     319             : 
     320             : {
     321        1595 :     if (GDALGetDriverByName("DIPEx") != nullptr)
     322         302 :         return;
     323             : 
     324        1293 :     GDALDriver *poDriver = new GDALDriver();
     325             : 
     326        1293 :     poDriver->SetDescription("DIPEx");
     327        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     328        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "DIPEx");
     329        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     330             : 
     331        1293 :     poDriver->pfnOpen = DIPExDataset::Open;
     332             : 
     333        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     334             : }