LCOV - code coverage report
Current view: top level - frmts/bsb - bsbdataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 253 329 76.9 %
Date: 2025-02-18 14:19:29 Functions: 19 22 86.4 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  BSB Reader
       4             :  * Purpose:  BSBDataset implementation for BSB format.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
       9             :  * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "bsb_read.h"
      15             : #include "cpl_string.h"
      16             : #include "gdal_frmts.h"
      17             : #include "gdal_pam.h"
      18             : #include "ogr_spatialref.h"
      19             : 
      20             : #include <cstdlib>
      21             : #include <algorithm>
      22             : 
      23             : // Disabled as people may worry about the BSB patent
      24             : // #define BSB_CREATE
      25             : 
      26             : /************************************************************************/
      27             : /* ==================================================================== */
      28             : /*                              BSBDataset                              */
      29             : /* ==================================================================== */
      30             : /************************************************************************/
      31             : 
      32             : class BSBRasterBand;
      33             : 
      34             : class BSBDataset final : public GDALPamDataset
      35             : {
      36             :     int nGCPCount;
      37             :     GDAL_GCP *pasGCPList;
      38             :     OGRSpatialReference m_oGCPSRS{};
      39             : 
      40             :     double adfGeoTransform[6];
      41             :     int bGeoTransformSet;
      42             : 
      43             :     void ScanForGCPs(bool isNos, const char *pszFilename);
      44             :     void ScanForGCPsNos(const char *pszFilename);
      45             :     void ScanForGCPsBSB();
      46             : 
      47             :     void ScanForCutline();
      48             : 
      49             :     static int IdentifyInternal(GDALOpenInfo *, bool &isNosOut);
      50             : 
      51             :   public:
      52             :     BSBDataset();
      53             :     ~BSBDataset() override;
      54             : 
      55             :     BSBInfo *psInfo;
      56             : 
      57             :     static GDALDataset *Open(GDALOpenInfo *);
      58             :     static int Identify(GDALOpenInfo *);
      59             : 
      60             :     int GetGCPCount() override;
      61             :     const OGRSpatialReference *GetSpatialRef() const override;
      62             :     const GDAL_GCP *GetGCPs() override;
      63             : 
      64             :     CPLErr GetGeoTransform(double *padfTransform) override;
      65             :     const OGRSpatialReference *GetGCPSpatialRef() const override;
      66             : };
      67             : 
      68             : /************************************************************************/
      69             : /* ==================================================================== */
      70             : /*                            BSBRasterBand                             */
      71             : /* ==================================================================== */
      72             : /************************************************************************/
      73             : 
      74             : class BSBRasterBand final : public GDALPamRasterBand
      75             : {
      76             :     GDALColorTable oCT;
      77             : 
      78             :   public:
      79             :     explicit BSBRasterBand(BSBDataset *);
      80             : 
      81             :     CPLErr IReadBlock(int, int, void *) override;
      82             :     GDALColorTable *GetColorTable() override;
      83             :     GDALColorInterp GetColorInterpretation() override;
      84             : };
      85             : 
      86             : /************************************************************************/
      87             : /*                           BSBRasterBand()                            */
      88             : /************************************************************************/
      89             : 
      90          13 : BSBRasterBand::BSBRasterBand(BSBDataset *poDSIn)
      91             : 
      92             : {
      93          13 :     poDS = poDSIn;
      94          13 :     nBand = 1;
      95             : 
      96          13 :     eDataType = GDT_Byte;
      97             : 
      98          13 :     nBlockXSize = poDS->GetRasterXSize();
      99          13 :     nBlockYSize = 1;
     100             : 
     101             :     // Note that the first color table entry is dropped, everything is
     102             :     // shifted down.
     103        1430 :     for (int i = 0; i < poDSIn->psInfo->nPCTSize - 1; i++)
     104             :     {
     105        1417 :         GDALColorEntry oColor = {poDSIn->psInfo->pabyPCT[i * 3 + 0 + 3],
     106        1417 :                                  poDSIn->psInfo->pabyPCT[i * 3 + 1 + 3],
     107        1417 :                                  poDSIn->psInfo->pabyPCT[i * 3 + 2 + 3], 255};
     108             : 
     109        1417 :         oCT.SetColorEntry(i, &oColor);
     110             :     }
     111          13 : }
     112             : 
     113             : /************************************************************************/
     114             : /*                             IReadBlock()                             */
     115             : /************************************************************************/
     116             : 
     117         250 : CPLErr BSBRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
     118             :                                  void *pImage)
     119             : {
     120         250 :     BSBDataset *poGDS = (BSBDataset *)poDS;
     121         250 :     GByte *pabyScanline = (GByte *)pImage;
     122             : 
     123         250 :     if (BSBReadScanline(poGDS->psInfo, nBlockYOff, pabyScanline))
     124             :     {
     125       12648 :         for (int i = 0; i < nBlockXSize; i++)
     126             :         {
     127             :             /* The indices start at 1, except in case of some charts */
     128             :             /* where there are missing values, which are filled to 0 */
     129             :             /* by BSBReadScanline */
     130       12400 :             if (pabyScanline[i] > 0)
     131       12400 :                 pabyScanline[i] -= 1;
     132             :         }
     133             : 
     134         248 :         return CE_None;
     135             :     }
     136             : 
     137           2 :     return CE_Failure;
     138             : }
     139             : 
     140             : /************************************************************************/
     141             : /*                           GetColorTable()                            */
     142             : /************************************************************************/
     143             : 
     144           0 : GDALColorTable *BSBRasterBand::GetColorTable()
     145             : 
     146             : {
     147           0 :     return &oCT;
     148             : }
     149             : 
     150             : /************************************************************************/
     151             : /*                       GetColorInterpretation()                       */
     152             : /************************************************************************/
     153             : 
     154           0 : GDALColorInterp BSBRasterBand::GetColorInterpretation()
     155             : 
     156             : {
     157           0 :     return GCI_PaletteIndex;
     158             : }
     159             : 
     160             : /************************************************************************/
     161             : /* ==================================================================== */
     162             : /*                              BSBDataset                              */
     163             : /* ==================================================================== */
     164             : /************************************************************************/
     165             : 
     166             : /************************************************************************/
     167             : /*                           BSBDataset()                               */
     168             : /************************************************************************/
     169             : 
     170          13 : BSBDataset::BSBDataset()
     171             :     : nGCPCount(0), pasGCPList(nullptr), bGeoTransformSet(FALSE),
     172          13 :       psInfo(nullptr)
     173             : {
     174          13 :     m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     175          13 :     m_oGCPSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
     176             : 
     177          13 :     adfGeoTransform[0] = 0.0; /* X Origin (top left corner) */
     178          13 :     adfGeoTransform[1] = 1.0; /* X Pixel size */
     179          13 :     adfGeoTransform[2] = 0.0;
     180          13 :     adfGeoTransform[3] = 0.0; /* Y Origin (top left corner) */
     181          13 :     adfGeoTransform[4] = 0.0;
     182          13 :     adfGeoTransform[5] = 1.0; /* Y Pixel Size */
     183          13 : }
     184             : 
     185             : /************************************************************************/
     186             : /*                            ~BSBDataset()                             */
     187             : /************************************************************************/
     188             : 
     189          26 : BSBDataset::~BSBDataset()
     190             : 
     191             : {
     192          13 :     FlushCache(true);
     193             : 
     194          13 :     GDALDeinitGCPs(nGCPCount, pasGCPList);
     195          13 :     CPLFree(pasGCPList);
     196             : 
     197          13 :     if (psInfo != nullptr)
     198          13 :         BSBClose(psInfo);
     199          26 : }
     200             : 
     201             : /************************************************************************/
     202             : /*                          GetGeoTransform()                           */
     203             : /************************************************************************/
     204             : 
     205           1 : CPLErr BSBDataset::GetGeoTransform(double *padfTransform)
     206             : 
     207             : {
     208           1 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     209             : 
     210           1 :     if (bGeoTransformSet)
     211           1 :         return CE_None;
     212             : 
     213           0 :     return CE_Failure;
     214             : }
     215             : 
     216             : /************************************************************************/
     217             : /*                          GetSpatialRef()                             */
     218             : /************************************************************************/
     219             : 
     220           1 : const OGRSpatialReference *BSBDataset::GetSpatialRef() const
     221             : 
     222             : {
     223           1 :     if (bGeoTransformSet)
     224           1 :         return &m_oGCPSRS;
     225             : 
     226           0 :     return nullptr;
     227             : }
     228             : 
     229             : /************************************************************************/
     230             : /*                     GDALHeuristicDatelineWrap()                      */
     231             : /************************************************************************/
     232             : 
     233           3 : static void GDALHeuristicDatelineWrap(int nPointCount, double *padfX)
     234             : 
     235             : {
     236           3 :     if (nPointCount < 2)
     237           3 :         return;
     238             : 
     239             :     /* -------------------------------------------------------------------- */
     240             :     /*      Work out what the longitude range will be centering on the      */
     241             :     /*      prime meridian (-180 to 180) and centering on the dateline      */
     242             :     /*      (0 to 360).                                                     */
     243             :     /* -------------------------------------------------------------------- */
     244             :     /* Following inits are useless but keep GCC happy */
     245           3 :     double dfX_PM_Min = 0.0;
     246           3 :     double dfX_PM_Max = 0.0;
     247           3 :     double dfX_Dateline_Min = 0.0;
     248           3 :     double dfX_Dateline_Max = 0.0;
     249             : 
     250         110 :     for (int i = 0; i < nPointCount; i++)
     251             :     {
     252         107 :         double dfX_PM = padfX[i];
     253         107 :         if (dfX_PM > 180)
     254           0 :             dfX_PM -= 360.0;
     255             : 
     256         107 :         double dfX_Dateline = padfX[i];
     257         107 :         if (dfX_Dateline < 0)
     258           0 :             dfX_Dateline += 360.0;
     259             : 
     260         107 :         if (i == 0)
     261             :         {
     262           3 :             dfX_PM_Min = dfX_PM;
     263           3 :             dfX_PM_Max = dfX_PM;
     264           3 :             dfX_Dateline_Min = dfX_Dateline;
     265           3 :             dfX_Dateline_Max = dfX_Dateline;
     266             :         }
     267             :         else
     268             :         {
     269         104 :             dfX_PM_Min = std::min(dfX_PM_Min, dfX_PM);
     270         104 :             dfX_PM_Max = std::max(dfX_PM_Max, dfX_PM);
     271         104 :             dfX_Dateline_Min = std::min(dfX_Dateline_Min, dfX_Dateline);
     272         104 :             dfX_Dateline_Max = std::max(dfX_Dateline_Max, dfX_Dateline);
     273             :         }
     274             :     }
     275             : 
     276             :     /* -------------------------------------------------------------------- */
     277             :     /*      Do nothing if the range is always fairly small - no apparent    */
     278             :     /*      wrapping issues.                                                */
     279             :     /* -------------------------------------------------------------------- */
     280           3 :     if ((dfX_PM_Max - dfX_PM_Min) < 270.0 &&
     281           3 :         (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0)
     282           3 :         return;
     283             : 
     284             :     /* -------------------------------------------------------------------- */
     285             :     /*      Do nothing if both approach have a wide range - best not to    */
     286             :     /*      fiddle if we aren't sure we are improving things.               */
     287             :     /* -------------------------------------------------------------------- */
     288           0 :     if ((dfX_PM_Max - dfX_PM_Min) > 270.0 &&
     289           0 :         (dfX_Dateline_Max - dfX_Dateline_Min) > 270.0)
     290           0 :         return;
     291             : 
     292             :     /* -------------------------------------------------------------------- */
     293             :     /*      Pick which way to transform things.                             */
     294             :     /* -------------------------------------------------------------------- */
     295             :     bool bUsePMWrap;
     296             : 
     297           0 :     if ((dfX_PM_Max - dfX_PM_Min) > 270.0 &&
     298           0 :         (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0)
     299           0 :         bUsePMWrap = false;
     300             :     else
     301           0 :         bUsePMWrap = true;
     302             : 
     303             :     /* -------------------------------------------------------------------- */
     304             :     /*      Apply rewrapping.                                               */
     305             :     /* -------------------------------------------------------------------- */
     306           0 :     for (int i = 0; i < nPointCount; i++)
     307             :     {
     308           0 :         if (bUsePMWrap)
     309             :         {
     310           0 :             if (padfX[i] > 180)
     311           0 :                 padfX[i] -= 360.0;
     312             :         }
     313             :         else
     314             :         {
     315           0 :             if (padfX[i] < 0)
     316           0 :                 padfX[i] += 360.0;
     317             :         }
     318             :     }
     319             : }
     320             : 
     321             : /************************************************************************/
     322             : /*                   GDALHeuristicDatelineWrapGCPs()                    */
     323             : /************************************************************************/
     324             : 
     325           3 : static void GDALHeuristicDatelineWrapGCPs(int nPointCount, GDAL_GCP *pasGCPList)
     326             : {
     327           6 :     std::vector<double> oadfX;
     328             : 
     329           3 :     oadfX.resize(nPointCount);
     330         110 :     for (int i = 0; i < nPointCount; i++)
     331         107 :         oadfX[i] = pasGCPList[i].dfGCPX;
     332             : 
     333           3 :     GDALHeuristicDatelineWrap(nPointCount, &(oadfX[0]));
     334             : 
     335         110 :     for (int i = 0; i < nPointCount; i++)
     336         107 :         pasGCPList[i].dfGCPX = oadfX[i];
     337           3 : }
     338             : 
     339             : /************************************************************************/
     340             : /*                            ScanForGCPs()                             */
     341             : /************************************************************************/
     342             : 
     343          13 : void BSBDataset::ScanForGCPs(bool isNos, const char *pszFilename)
     344             : 
     345             : {
     346             :     /* -------------------------------------------------------------------- */
     347             :     /*      Collect GCPs as appropriate to source.                          */
     348             :     /* -------------------------------------------------------------------- */
     349          13 :     nGCPCount = 0;
     350             : 
     351          13 :     if (isNos)
     352             :     {
     353           0 :         ScanForGCPsNos(pszFilename);
     354             :     }
     355             :     else
     356             :     {
     357          13 :         ScanForGCPsBSB();
     358             :     }
     359             : 
     360             :     /* -------------------------------------------------------------------- */
     361             :     /*      Apply heuristics to re-wrap GCPs to maintain continuity        */
     362             :     /*      over the international dateline.                                */
     363             :     /* -------------------------------------------------------------------- */
     364          13 :     if (nGCPCount > 1)
     365           3 :         GDALHeuristicDatelineWrapGCPs(nGCPCount, pasGCPList);
     366             : 
     367             :     /* -------------------------------------------------------------------- */
     368             :     /*      Collect coordinate system related parameters from header.       */
     369             :     /* -------------------------------------------------------------------- */
     370          13 :     const char *pszKNP = nullptr;
     371          13 :     const char *pszKNQ = nullptr;
     372             : 
     373        1765 :     for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
     374             :     {
     375        1752 :         if (STARTS_WITH_CI(psInfo->papszHeader[i], "KNP/"))
     376             :         {
     377          13 :             pszKNP = psInfo->papszHeader[i];
     378          13 :             SetMetadataItem("BSB_KNP", pszKNP + 4);
     379             :         }
     380        1752 :         if (STARTS_WITH_CI(psInfo->papszHeader[i], "KNQ/"))
     381             :         {
     382           3 :             pszKNQ = psInfo->papszHeader[i];
     383           3 :             SetMetadataItem("BSB_KNQ", pszKNQ + 4);
     384             :         }
     385             :     }
     386             : 
     387             :     /* -------------------------------------------------------------------- */
     388             :     /*      Can we derive a reasonable coordinate system definition for     */
     389             :     /*      this file?  For now we keep it simple, just handling            */
     390             :     /*      mercator. In the future we should consider others.              */
     391             :     /* -------------------------------------------------------------------- */
     392          26 :     CPLString osUnderlyingSRS;
     393          13 :     if (pszKNP != nullptr)
     394             :     {
     395          13 :         const char *pszPR = strstr(pszKNP, "PR=");
     396          13 :         const char *pszGD = strstr(pszKNP, "GD=");
     397          13 :         const char *pszGEOGCS = SRS_WKT_WGS84_LAT_LONG;
     398          26 :         CPLString osPP;
     399             : 
     400             :         // Capture the PP string.
     401          13 :         const char *pszValue = strstr(pszKNP, "PP=");
     402          13 :         const char *pszEnd = pszValue ? strstr(pszValue, ",") : nullptr;
     403          13 :         if (pszValue && pszEnd)
     404          13 :             osPP.assign(pszValue + 3, pszEnd - pszValue - 3);
     405             : 
     406             :         // Look at the datum
     407          13 :         if (pszGD == nullptr)
     408             :         {
     409             :             /* no match. We'll default to EPSG:4326 */
     410             :         }
     411          13 :         else if (STARTS_WITH_CI(pszGD, "GD=European 1950"))
     412             :         {
     413           0 :             pszGEOGCS =
     414             :                 "GEOGCS[\"ED50\",DATUM[\"European_Datum_1950\",SPHEROID["
     415             :                 "\"International "
     416             :                 "1924\",6378388,297,AUTHORITY[\"EPSG\",\"7022\"]],TOWGS84[-87,-"
     417             :                 "98,-121,0,0,0,0],AUTHORITY[\"EPSG\",\"6230\"]],PRIMEM["
     418             :                 "\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\","
     419             :                 "0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY["
     420             :                 "\"EPSG\",\"4230\"]]";
     421             :         }
     422             : 
     423             :         // Look at the projection
     424          13 :         if (pszPR == nullptr)
     425             :         {
     426             :             /* no match */
     427             :         }
     428          13 :         else if (STARTS_WITH_CI(pszPR, "PR=MERCATOR") && nGCPCount > 0)
     429             :         {
     430             :             // We somewhat arbitrarily select our first GCPX as our
     431             :             // central meridian.  This is mostly helpful to ensure
     432             :             // that regions crossing the dateline will be contiguous
     433             :             // in mercator.
     434           1 :             osUnderlyingSRS.Printf(
     435             :                 "PROJCS[\"Global "
     436             :                 "Mercator\",%s,PROJECTION[\"Mercator_2SP\"],PARAMETER["
     437             :                 "\"standard_parallel_1\",0],PARAMETER[\"latitude_of_origin\",0]"
     438             :                 ",PARAMETER[\"central_meridian\",%d],PARAMETER[\"false_"
     439             :                 "easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]"
     440             :                 "]",
     441           1 :                 pszGEOGCS, (int)pasGCPList[0].dfGCPX);
     442             :         }
     443             : 
     444          13 :         else if (STARTS_WITH_CI(pszPR, "PR=TRANSVERSE MERCATOR") &&
     445           1 :                  !osPP.empty())
     446             :         {
     447             : 
     448             :             osUnderlyingSRS.Printf(
     449             :                 "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],"
     450             :                 "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_"
     451             :                 "meridian\",%s],PARAMETER[\"scale_factor\",1],PARAMETER["
     452             :                 "\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT["
     453             :                 "\"Meter\",1]]",
     454           1 :                 pszGEOGCS, osPP.c_str());
     455             :         }
     456             : 
     457          11 :         else if (STARTS_WITH_CI(pszPR, "PR=UNIVERSAL TRANSVERSE MERCATOR") &&
     458           0 :                  !osPP.empty())
     459             :         {
     460             :             // This is not *really* UTM unless the central meridian
     461             :             // matches a zone which it does not in some (most?) maps.
     462             :             osUnderlyingSRS.Printf(
     463             :                 "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],"
     464             :                 "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_"
     465             :                 "meridian\",%s],PARAMETER[\"scale_factor\",0.9996],PARAMETER["
     466             :                 "\"false_easting\",500000],PARAMETER[\"false_northing\",0],"
     467             :                 "UNIT[\"Meter\",1]]",
     468           0 :                 pszGEOGCS, osPP.c_str());
     469             :         }
     470             : 
     471          11 :         else if (STARTS_WITH_CI(pszPR, "PR=POLYCONIC") && !osPP.empty())
     472             :         {
     473             :             osUnderlyingSRS.Printf(
     474             :                 "PROJCS[\"unnamed\",%s,PROJECTION[\"Polyconic\"],PARAMETER["
     475             :                 "\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%s],"
     476             :                 "PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0]"
     477             :                 ",UNIT[\"Meter\",1]]",
     478           0 :                 pszGEOGCS, osPP.c_str());
     479             :         }
     480             : 
     481          23 :         else if (STARTS_WITH_CI(pszPR, "PR=LAMBERT CONFORMAL CONIC") &&
     482          11 :                  !osPP.empty() && pszKNQ != nullptr)
     483             :         {
     484           2 :             CPLString osP2, osP3;
     485             : 
     486             :             // Capture the KNQ/P2 string.
     487           1 :             pszValue = strstr(pszKNQ, "P2=");
     488           1 :             if (pszValue)
     489           1 :                 pszEnd = strstr(pszValue, ",");
     490           1 :             if (pszValue && pszEnd)
     491           1 :                 osP2.assign(pszValue + 3, pszEnd - pszValue - 3);
     492             : 
     493             :             // Capture the KNQ/P3 string.
     494           1 :             pszValue = strstr(pszKNQ, "P3=");
     495           1 :             if (pszValue)
     496           1 :                 pszEnd = strstr(pszValue, ",");
     497           1 :             if (pszValue)
     498             :             {
     499           1 :                 if (pszEnd)
     500           1 :                     osP3.assign(pszValue + 3, pszEnd - pszValue - 3);
     501             :                 else
     502           0 :                     osP3.assign(pszValue + 3);
     503             :             }
     504             : 
     505           1 :             if (!osP2.empty() && !osP3.empty())
     506             :                 osUnderlyingSRS.Printf(
     507             :                     "PROJCS[\"unnamed\",%s,PROJECTION[\"Lambert_Conformal_"
     508             :                     "Conic_2SP\"],PARAMETER[\"standard_parallel_1\",%s],"
     509             :                     "PARAMETER[\"standard_parallel_2\",%s],PARAMETER["
     510             :                     "\"latitude_of_origin\",0.0],PARAMETER[\"central_"
     511             :                     "meridian\",%s],PARAMETER[\"false_easting\",0.0],PARAMETER["
     512             :                     "\"false_northing\",0.0],UNIT[\"Meter\",1]]",
     513           1 :                     pszGEOGCS, osP2.c_str(), osP3.c_str(), osPP.c_str());
     514             :         }
     515             :     }
     516             : 
     517             :     /* -------------------------------------------------------------------- */
     518             :     /*      If we got an alternate underlying coordinate system, try        */
     519             :     /*      converting the GCPs to that coordinate system.                  */
     520             :     /* -------------------------------------------------------------------- */
     521          13 :     if (osUnderlyingSRS.length() > 0)
     522             :     {
     523           6 :         OGRSpatialReference oGeog_SRS, oProjected_SRS;
     524             : 
     525           3 :         oProjected_SRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     526           3 :         oProjected_SRS.SetFromUserInput(osUnderlyingSRS);
     527           3 :         oGeog_SRS.CopyGeogCSFrom(&oProjected_SRS);
     528           3 :         oGeog_SRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     529             : 
     530             :         OGRCoordinateTransformation *poCT =
     531           3 :             OGRCreateCoordinateTransformation(&oGeog_SRS, &oProjected_SRS);
     532           3 :         if (poCT != nullptr)
     533             :         {
     534         110 :             for (int i = 0; i < nGCPCount; i++)
     535             :             {
     536         107 :                 poCT->Transform(1, &(pasGCPList[i].dfGCPX),
     537         107 :                                 &(pasGCPList[i].dfGCPY),
     538         107 :                                 &(pasGCPList[i].dfGCPZ));
     539             :             }
     540             : 
     541           3 :             m_oGCPSRS.importFromWkt(osUnderlyingSRS.c_str());
     542             : 
     543           3 :             delete poCT;
     544             :         }
     545             :         else
     546           0 :             CPLErrorReset();
     547             :     }
     548             : 
     549             :     /* -------------------------------------------------------------------- */
     550             :     /*      Attempt to prepare a geotransform from the GCPs.                */
     551             :     /* -------------------------------------------------------------------- */
     552          13 :     if (GDALGCPsToGeoTransform(nGCPCount, pasGCPList, adfGeoTransform, FALSE))
     553             :     {
     554           2 :         bGeoTransformSet = TRUE;
     555             :     }
     556          13 : }
     557             : 
     558             : /************************************************************************/
     559             : /*                           ScanForGCPsNos()                           */
     560             : /*                                                                      */
     561             : /*      Nos files have an accompanying .geo file, that contains some    */
     562             : /*      of the information normally contained in the header section     */
     563             : /*      with BSB files. we try and open a file with the same name,      */
     564             : /*      but a .geo extension, and look for lines like...                */
     565             : /*      PointX=long lat line pixel    (using the same naming system     */
     566             : /*      as BSB) Point1=-22.0000 64.250000 197 744                       */
     567             : /************************************************************************/
     568             : 
     569           0 : void BSBDataset::ScanForGCPsNos(const char *pszFilename)
     570             : {
     571           0 :     const std::string extension = CPLGetExtensionSafe(pszFilename);
     572             : 
     573             :     // pseudointelligently try and guess whether we want a .geo or a .GEO
     574           0 :     std::string geofile;
     575           0 :     if (extension.size() >= 2 && extension[1] == 'O')
     576             :     {
     577           0 :         geofile = CPLResetExtensionSafe(pszFilename, "GEO");
     578             :     }
     579             :     else
     580             :     {
     581           0 :         geofile = CPLResetExtensionSafe(pszFilename, "geo");
     582             :     }
     583             : 
     584           0 :     FILE *gfp = VSIFOpen(geofile.c_str(), "r");  // Text files
     585           0 :     if (gfp == nullptr)
     586             :     {
     587           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     588             :                  "Couldn't find a matching .GEO file: %s", geofile.c_str());
     589           0 :         return;
     590             :     }
     591             : 
     592           0 :     char *thisLine = (char *)CPLMalloc(80);  // FIXME
     593             : 
     594             :     // Count the GCPs (reference points) and seek the file pointer 'gfp' to the
     595             :     // starting point
     596           0 :     int fileGCPCount = 0;
     597           0 :     while (fgets(thisLine, 80, gfp))
     598             :     {
     599           0 :         if (STARTS_WITH_CI(thisLine, "Point"))
     600           0 :             fileGCPCount++;
     601             :     }
     602           0 :     VSIRewind(gfp);
     603             : 
     604             :     // Memory has not been allocated to fileGCPCount yet
     605           0 :     pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), fileGCPCount + 1);
     606             : 
     607           0 :     while (fgets(thisLine, 80, gfp))
     608             :     {
     609           0 :         if (STARTS_WITH_CI(thisLine, "Point"))
     610             :         {
     611             :             // got a point line, turn it into a gcp
     612             :             char **Tokens =
     613           0 :                 CSLTokenizeStringComplex(thisLine, "= ", FALSE, FALSE);
     614           0 :             if (CSLCount(Tokens) >= 5)
     615             :             {
     616           0 :                 GDALInitGCPs(1, pasGCPList + nGCPCount);
     617           0 :                 pasGCPList[nGCPCount].dfGCPX = CPLAtof(Tokens[1]);
     618           0 :                 pasGCPList[nGCPCount].dfGCPY = CPLAtof(Tokens[2]);
     619           0 :                 pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(Tokens[4]);
     620           0 :                 pasGCPList[nGCPCount].dfGCPLine = CPLAtof(Tokens[3]);
     621             : 
     622           0 :                 CPLFree(pasGCPList[nGCPCount].pszId);
     623             :                 char szName[50];
     624           0 :                 snprintf(szName, sizeof(szName), "GCP_%d", nGCPCount + 1);
     625           0 :                 pasGCPList[nGCPCount].pszId = CPLStrdup(szName);
     626             : 
     627           0 :                 nGCPCount++;
     628             :             }
     629           0 :             CSLDestroy(Tokens);
     630             :         }
     631             :     }
     632             : 
     633           0 :     CPLFree(thisLine);
     634           0 :     VSIFClose(gfp);
     635             : }
     636             : 
     637             : /************************************************************************/
     638             : /*                            ScanForGCPsBSB()                          */
     639             : /************************************************************************/
     640             : 
     641          13 : void BSBDataset::ScanForGCPsBSB()
     642             : {
     643             :     /* -------------------------------------------------------------------- */
     644             :     /*      Collect standalone GCPs.  They look like:                       */
     645             :     /*                                                                      */
     646             :     /*      REF/1,115,2727,32.346666666667,-60.881666666667                 */
     647             :     /*      REF/n,pixel,line,lat,long                                       */
     648             :     /* -------------------------------------------------------------------- */
     649          13 :     int fileGCPCount = 0;
     650             : 
     651             :     // Count the GCPs (reference points) in psInfo->papszHeader
     652        1765 :     for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
     653        1752 :         if (STARTS_WITH_CI(psInfo->papszHeader[i], "REF/"))
     654         107 :             fileGCPCount++;
     655             : 
     656             :     // Memory has not been allocated to fileGCPCount yet
     657          13 :     pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), fileGCPCount + 1);
     658             : 
     659        1765 :     for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
     660             :     {
     661             : 
     662        1752 :         if (!STARTS_WITH_CI(psInfo->papszHeader[i], "REF/"))
     663        1645 :             continue;
     664             : 
     665         214 :         char **papszTokens = CSLTokenizeStringComplex(
     666         107 :             psInfo->papszHeader[i] + 4, ",", FALSE, FALSE);
     667             : 
     668         107 :         if (CSLCount(papszTokens) > 4)
     669             :         {
     670         107 :             GDALInitGCPs(1, pasGCPList + nGCPCount);
     671             : 
     672         107 :             pasGCPList[nGCPCount].dfGCPX = CPLAtof(papszTokens[4]);
     673         107 :             pasGCPList[nGCPCount].dfGCPY = CPLAtof(papszTokens[3]);
     674         107 :             pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(papszTokens[1]);
     675         107 :             pasGCPList[nGCPCount].dfGCPLine = CPLAtof(papszTokens[2]);
     676             : 
     677         107 :             CPLFree(pasGCPList[nGCPCount].pszId);
     678         107 :             if (CSLCount(papszTokens) > 5)
     679             :             {
     680           0 :                 pasGCPList[nGCPCount].pszId = CPLStrdup(papszTokens[5]);
     681             :             }
     682             :             else
     683             :             {
     684             :                 char szName[50];
     685         107 :                 snprintf(szName, sizeof(szName), "GCP_%d", nGCPCount + 1);
     686         107 :                 pasGCPList[nGCPCount].pszId = CPLStrdup(szName);
     687             :             }
     688             : 
     689         107 :             nGCPCount++;
     690             :         }
     691         107 :         CSLDestroy(papszTokens);
     692             :     }
     693          13 : }
     694             : 
     695             : /************************************************************************/
     696             : /*                            ScanForCutline()                          */
     697             : /************************************************************************/
     698             : 
     699          13 : void BSBDataset::ScanForCutline()
     700             : {
     701             :     /* PLY: Border Polygon Record - coordinates of the panel within the
     702             :      * raster image, given in chart datum lat/long. Any shape polygon.
     703             :      * They look like:
     704             :      *      PLY/1,32.346666666667,-60.881666666667
     705             :      *      PLY/n,lat,long
     706             :      *
     707             :      * If found then we return it via a BSB_CUTLINE metadata item as a WKT
     708             :      * POLYGON.
     709             :      */
     710             : 
     711          26 :     std::string wkt;
     712        1765 :     for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
     713             :     {
     714        1752 :         if (!STARTS_WITH_CI(psInfo->papszHeader[i], "PLY/"))
     715        1731 :             continue;
     716             : 
     717             :         const CPLStringList aosTokens(
     718          42 :             CSLTokenizeString2(psInfo->papszHeader[i] + 4, ",", 0));
     719             : 
     720          21 :         if (aosTokens.size() >= 3)
     721             :         {
     722          21 :             if (wkt.empty())
     723           2 :                 wkt = "POLYGON ((";
     724             :             else
     725          19 :                 wkt += ',';
     726          21 :             wkt += aosTokens[2];
     727          21 :             wkt += ' ';
     728          21 :             wkt += aosTokens[1];
     729             :         }
     730             :     }
     731             : 
     732          13 :     if (!wkt.empty())
     733             :     {
     734           2 :         wkt += "))";
     735           2 :         SetMetadataItem("BSB_CUTLINE", wkt.c_str());
     736             :     }
     737          13 : }
     738             : 
     739             : /************************************************************************/
     740             : /*                          IdentifyInternal()                          */
     741             : /************************************************************************/
     742             : 
     743       56449 : int BSBDataset::IdentifyInternal(GDALOpenInfo *poOpenInfo, bool &isNosOut)
     744             : 
     745             : {
     746             :     /* -------------------------------------------------------------------- */
     747             :     /*      Check for BSB/ keyword.                                         */
     748             :     /* -------------------------------------------------------------------- */
     749       56449 :     isNosOut = false;
     750             : 
     751       56449 :     if (poOpenInfo->nHeaderBytes < 1000)
     752       52141 :         return FALSE;
     753             : 
     754        4308 :     int i = 0;
     755     5653250 :     for (; i < poOpenInfo->nHeaderBytes - 4; i++)
     756             :     {
     757     5648960 :         if (poOpenInfo->pabyHeader[i + 0] == 'B' &&
     758       16793 :             poOpenInfo->pabyHeader[i + 1] == 'S' &&
     759         102 :             poOpenInfo->pabyHeader[i + 2] == 'B' &&
     760          28 :             poOpenInfo->pabyHeader[i + 3] == '/')
     761          26 :             break;
     762     5648940 :         if (poOpenInfo->pabyHeader[i + 0] == 'N' &&
     763       17538 :             poOpenInfo->pabyHeader[i + 1] == 'O' &&
     764         622 :             poOpenInfo->pabyHeader[i + 2] == 'S' &&
     765          19 :             poOpenInfo->pabyHeader[i + 3] == '/')
     766             :         {
     767           0 :             isNosOut = true;
     768           0 :             break;
     769             :         }
     770     5648940 :         if (poOpenInfo->pabyHeader[i + 0] == 'W' &&
     771        5160 :             poOpenInfo->pabyHeader[i + 1] == 'X' &&
     772         142 :             poOpenInfo->pabyHeader[i + 2] == '\\' &&
     773           0 :             poOpenInfo->pabyHeader[i + 3] == '8')
     774           0 :             break;
     775             :     }
     776             : 
     777        4308 :     if (i == poOpenInfo->nHeaderBytes - 4)
     778        4283 :         return FALSE;
     779             : 
     780             :     /* Additional test to avoid false positive. See #2881 */
     781          25 :     const char *pszHeader =
     782             :         reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
     783          25 :     const char *pszShiftedHeader = pszHeader + i;
     784          25 :     const char *pszRA = strstr(pszShiftedHeader, "RA=");
     785          25 :     if (pszRA == nullptr) /* This may be a NO1 file */
     786           0 :         pszRA = strstr(pszShiftedHeader, "[JF");
     787          25 :     if (pszRA == nullptr)
     788           0 :         return FALSE;
     789          25 :     if (pszRA - pszShiftedHeader > 100 && !strstr(pszHeader, "VER/") &&
     790           0 :         !strstr(pszHeader, "KNP/") && !strstr(pszHeader, "KNQ/") &&
     791           0 :         !strstr(pszHeader, "RGB/"))
     792           0 :         return FALSE;
     793             : 
     794          25 :     return TRUE;
     795             : }
     796             : 
     797             : /************************************************************************/
     798             : /*                              Identify()                              */
     799             : /************************************************************************/
     800             : 
     801       56436 : int BSBDataset::Identify(GDALOpenInfo *poOpenInfo)
     802             : 
     803             : {
     804             :     bool isNos;
     805       56436 :     return IdentifyInternal(poOpenInfo, isNos);
     806             : }
     807             : 
     808             : /************************************************************************/
     809             : /*                                Open()                                */
     810             : /************************************************************************/
     811             : 
     812          13 : GDALDataset *BSBDataset::Open(GDALOpenInfo *poOpenInfo)
     813             : 
     814             : {
     815          13 :     bool isNos = false;
     816          13 :     if (!IdentifyInternal(poOpenInfo, isNos))
     817           0 :         return nullptr;
     818             : 
     819          13 :     if (poOpenInfo->eAccess == GA_Update)
     820             :     {
     821           0 :         ReportUpdateNotSupportedByDriver("BSB");
     822           0 :         return nullptr;
     823             :     }
     824             : 
     825             :     /* -------------------------------------------------------------------- */
     826             :     /*      Create a corresponding GDALDataset.                             */
     827             :     /* -------------------------------------------------------------------- */
     828          13 :     BSBDataset *poDS = new BSBDataset();
     829             : 
     830             :     /* -------------------------------------------------------------------- */
     831             :     /*      Open the file.                                                  */
     832             :     /* -------------------------------------------------------------------- */
     833          13 :     poDS->psInfo = BSBOpen(poOpenInfo->pszFilename);
     834          13 :     if (poDS->psInfo == nullptr)
     835             :     {
     836           0 :         delete poDS;
     837           0 :         return nullptr;
     838             :     }
     839             : 
     840          13 :     poDS->nRasterXSize = poDS->psInfo->nXSize;
     841          13 :     poDS->nRasterYSize = poDS->psInfo->nYSize;
     842             : 
     843             :     /* -------------------------------------------------------------------- */
     844             :     /*      Create band information objects.                                */
     845             :     /* -------------------------------------------------------------------- */
     846          13 :     poDS->SetBand(1, new BSBRasterBand(poDS));
     847             : 
     848          13 :     poDS->ScanForGCPs(isNos, poOpenInfo->pszFilename);
     849             : 
     850             :     /* -------------------------------------------------------------------- */
     851             :     /*      Set CUTLINE metadata if a bounding polygon is available         */
     852             :     /* -------------------------------------------------------------------- */
     853          13 :     poDS->ScanForCutline();
     854             : 
     855             :     /* -------------------------------------------------------------------- */
     856             :     /*      Initialize any PAM information.                                 */
     857             :     /* -------------------------------------------------------------------- */
     858          13 :     poDS->SetDescription(poOpenInfo->pszFilename);
     859          13 :     poDS->TryLoadXML();
     860             : 
     861          13 :     poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
     862             : 
     863          13 :     return poDS;
     864             : }
     865             : 
     866             : /************************************************************************/
     867             : /*                            GetGCPCount()                             */
     868             : /************************************************************************/
     869             : 
     870           2 : int BSBDataset::GetGCPCount()
     871             : 
     872             : {
     873           2 :     return nGCPCount;
     874             : }
     875             : 
     876             : /************************************************************************/
     877             : /*                          GetGCPSpatialRef()                          */
     878             : /************************************************************************/
     879             : 
     880           1 : const OGRSpatialReference *BSBDataset::GetGCPSpatialRef() const
     881             : 
     882             : {
     883           1 :     return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
     884             : }
     885             : 
     886             : /************************************************************************/
     887             : /*                               GetGCP()                               */
     888             : /************************************************************************/
     889             : 
     890           1 : const GDAL_GCP *BSBDataset::GetGCPs()
     891             : 
     892             : {
     893           1 :     return pasGCPList;
     894             : }
     895             : 
     896             : #ifdef BSB_CREATE
     897             : 
     898             : /************************************************************************/
     899             : /*                             BSBIsSRSOK()                             */
     900             : /************************************************************************/
     901             : 
     902             : static int BSBIsSRSOK(const char *pszWKT)
     903             : {
     904             :     bool bOK = false;
     905             :     OGRSpatialReference oSRS, oSRS_WGS84, oSRS_NAD83;
     906             : 
     907             :     if (pszWKT != NULL && pszWKT[0] != '\0')
     908             :     {
     909             :         oSRS.importFromWkt(pszWKT);
     910             : 
     911             :         oSRS_WGS84.SetWellKnownGeogCS("WGS84");
     912             :         oSRS_NAD83.SetWellKnownGeogCS("NAD83");
     913             :         if ((oSRS.IsSameGeogCS(&oSRS_WGS84) ||
     914             :              oSRS.IsSameGeogCS(&oSRS_NAD83)) &&
     915             :             oSRS.IsGeographic() && oSRS.GetPrimeMeridian() == 0.0)
     916             :         {
     917             :             bOK = true;
     918             :         }
     919             :     }
     920             : 
     921             :     if (!bOK)
     922             :     {
     923             :         CPLError(CE_Warning, CPLE_NotSupported,
     924             :                  "BSB only supports WGS84 or NAD83 geographic projections.\n");
     925             :     }
     926             : 
     927             :     return bOK;
     928             : }
     929             : 
     930             : /************************************************************************/
     931             : /*                           BSBCreateCopy()                            */
     932             : /************************************************************************/
     933             : 
     934             : static GDALDataset *BSBCreateCopy(const char *pszFilename, GDALDataset *poSrcDS,
     935             :                                   int bStrict, char ** /*papszOptions*/,
     936             :                                   GDALProgressFunc /*pfnProgress*/,
     937             :                                   void * /*pProgressData*/)
     938             : 
     939             : {
     940             :     /* -------------------------------------------------------------------- */
     941             :     /*      Some some rudimentary checks                                    */
     942             :     /* -------------------------------------------------------------------- */
     943             :     const int nBands = poSrcDS->GetRasterCount();
     944             :     if (nBands != 1)
     945             :     {
     946             :         CPLError(CE_Failure, CPLE_NotSupported,
     947             :                  "BSB driver only supports one band images.\n");
     948             : 
     949             :         return NULL;
     950             :     }
     951             : 
     952             :     if (poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte && bStrict)
     953             :     {
     954             :         CPLError(CE_Failure, CPLE_NotSupported,
     955             :                  "BSB driver doesn't support data type %s. "
     956             :                  "Only eight bit bands supported.\n",
     957             :                  GDALGetDataTypeName(
     958             :                      poSrcDS->GetRasterBand(1)->GetRasterDataType()));
     959             : 
     960             :         return NULL;
     961             :     }
     962             : 
     963             :     const int nXSize = poSrcDS->GetRasterXSize();
     964             :     const int nYSize = poSrcDS->GetRasterYSize();
     965             : 
     966             :     /* -------------------------------------------------------------------- */
     967             :     /*      Open the output file.                                           */
     968             :     /* -------------------------------------------------------------------- */
     969             :     BSBInfo *psBSB = BSBCreate(pszFilename, 0, 200, nXSize, nYSize);
     970             :     if (psBSB == NULL)
     971             :         return NULL;
     972             : 
     973             :     /* -------------------------------------------------------------------- */
     974             :     /*      Prepare initial color table.colortable.                         */
     975             :     /* -------------------------------------------------------------------- */
     976             :     GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);
     977             :     unsigned char abyPCT[771];
     978             :     int nPCTSize;
     979             :     int anRemap[256];
     980             : 
     981             :     abyPCT[0] = 0;
     982             :     abyPCT[1] = 0;
     983             :     abyPCT[2] = 0;
     984             : 
     985             :     if (poBand->GetColorTable() == NULL)
     986             :     {
     987             :         /* map greyscale down to 63 grey levels. */
     988             :         for (int iColor = 0; iColor < 256; iColor++)
     989             :         {
     990             :             int nOutValue = (int)(iColor / 4.1) + 1;
     991             : 
     992             :             anRemap[iColor] = nOutValue;
     993             :             abyPCT[nOutValue * 3 + 0] = (unsigned char)iColor;
     994             :             abyPCT[nOutValue * 3 + 1] = (unsigned char)iColor;
     995             :             abyPCT[nOutValue * 3 + 2] = (unsigned char)iColor;
     996             :         }
     997             :         nPCTSize = 64;
     998             :     }
     999             :     else
    1000             :     {
    1001             :         GDALColorTable *poCT = poBand->GetColorTable();
    1002             :         int nColorTableSize = poCT->GetColorEntryCount();
    1003             :         if (nColorTableSize > 255)
    1004             :             nColorTableSize = 255;
    1005             : 
    1006             :         for (int iColor = 0; iColor < nColorTableSize; iColor++)
    1007             :         {
    1008             :             GDALColorEntry sEntry;
    1009             : 
    1010             :             poCT->GetColorEntryAsRGB(iColor, &sEntry);
    1011             : 
    1012             :             anRemap[iColor] = iColor + 1;
    1013             :             abyPCT[(iColor + 1) * 3 + 0] = (unsigned char)sEntry.c1;
    1014             :             abyPCT[(iColor + 1) * 3 + 1] = (unsigned char)sEntry.c2;
    1015             :             abyPCT[(iColor + 1) * 3 + 2] = (unsigned char)sEntry.c3;
    1016             :         }
    1017             : 
    1018             :         nPCTSize = nColorTableSize + 1;
    1019             : 
    1020             :         // Add entries for pixel values which apparently will not occur.
    1021             :         for (int iColor = nPCTSize; iColor < 256; iColor++)
    1022             :             anRemap[iColor] = 1;
    1023             :     }
    1024             : 
    1025             :     /* -------------------------------------------------------------------- */
    1026             :     /*      Boil out all duplicate entries.                                 */
    1027             :     /* -------------------------------------------------------------------- */
    1028             :     for (int i = 1; i < nPCTSize - 1; i++)
    1029             :     {
    1030             :         for (int j = i + 1; j < nPCTSize; j++)
    1031             :         {
    1032             :             if (abyPCT[i * 3 + 0] == abyPCT[j * 3 + 0] &&
    1033             :                 abyPCT[i * 3 + 1] == abyPCT[j * 3 + 1] &&
    1034             :                 abyPCT[i * 3 + 2] == abyPCT[j * 3 + 2])
    1035             :             {
    1036             :                 nPCTSize--;
    1037             :                 abyPCT[j * 3 + 0] = abyPCT[nPCTSize * 3 + 0];
    1038             :                 abyPCT[j * 3 + 1] = abyPCT[nPCTSize * 3 + 1];
    1039             :                 abyPCT[j * 3 + 2] = abyPCT[nPCTSize * 3 + 2];
    1040             : 
    1041             :                 for (int k = 0; k < 256; k++)
    1042             :                 {
    1043             :                     // merge matching entries.
    1044             :                     if (anRemap[k] == j)
    1045             :                         anRemap[k] = i;
    1046             : 
    1047             :                     // shift the last PCT entry into the new hole.
    1048             :                     if (anRemap[k] == nPCTSize)
    1049             :                         anRemap[k] = j;
    1050             :                 }
    1051             :             }
    1052             :         }
    1053             :     }
    1054             : 
    1055             :     /* -------------------------------------------------------------------- */
    1056             :     /*      Boil out all duplicate entries.                                 */
    1057             :     /* -------------------------------------------------------------------- */
    1058             :     if (nPCTSize > 128)
    1059             :     {
    1060             :         CPLError(CE_Warning, CPLE_AppDefined,
    1061             :                  "Having to merge color table entries to reduce %d real\n"
    1062             :                  "color table entries down to 127 values.",
    1063             :                  nPCTSize);
    1064             :     }
    1065             : 
    1066             :     while (nPCTSize > 128)
    1067             :     {
    1068             :         int nBestRange = 768;
    1069             :         int iBestMatch1 = -1;
    1070             :         int iBestMatch2 = -1;
    1071             : 
    1072             :         // Find the closest pair of color table entries.
    1073             : 
    1074             :         for (int i = 1; i < nPCTSize - 1; i++)
    1075             :         {
    1076             :             for (int j = i + 1; j < nPCTSize; j++)
    1077             :             {
    1078             :                 int nRange = std::abs(abyPCT[i * 3 + 0] - abyPCT[j * 3 + 0]) +
    1079             :                              std::abs(abyPCT[i * 3 + 1] - abyPCT[j * 3 + 1]) +
    1080             :                              std::abs(abyPCT[i * 3 + 2] - abyPCT[j * 3 + 2]);
    1081             : 
    1082             :                 if (nRange < nBestRange)
    1083             :                 {
    1084             :                     iBestMatch1 = i;
    1085             :                     iBestMatch2 = j;
    1086             :                     nBestRange = nRange;
    1087             :                 }
    1088             :             }
    1089             :         }
    1090             : 
    1091             :         // Merge the second entry into the first.
    1092             :         nPCTSize--;
    1093             :         abyPCT[iBestMatch2 * 3 + 0] = abyPCT[nPCTSize * 3 + 0];
    1094             :         abyPCT[iBestMatch2 * 3 + 1] = abyPCT[nPCTSize * 3 + 1];
    1095             :         abyPCT[iBestMatch2 * 3 + 2] = abyPCT[nPCTSize * 3 + 2];
    1096             : 
    1097             :         for (int i = 0; i < 256; i++)
    1098             :         {
    1099             :             // merge matching entries.
    1100             :             if (anRemap[i] == iBestMatch2)
    1101             :                 anRemap[i] = iBestMatch1;
    1102             : 
    1103             :             // shift the last PCT entry into the new hole.
    1104             :             if (anRemap[i] == nPCTSize)
    1105             :                 anRemap[i] = iBestMatch2;
    1106             :         }
    1107             :     }
    1108             : 
    1109             :     /* -------------------------------------------------------------------- */
    1110             :     /*      Write the PCT.                                                  */
    1111             :     /* -------------------------------------------------------------------- */
    1112             :     if (!BSBWritePCT(psBSB, nPCTSize, abyPCT))
    1113             :     {
    1114             :         BSBClose(psBSB);
    1115             :         return NULL;
    1116             :     }
    1117             : 
    1118             :     /* -------------------------------------------------------------------- */
    1119             :     /*      Write the GCPs.                                                 */
    1120             :     /* -------------------------------------------------------------------- */
    1121             :     double adfGeoTransform[6];
    1122             :     int nGCPCount = poSrcDS->GetGCPCount();
    1123             :     if (nGCPCount)
    1124             :     {
    1125             :         const char *pszGCPProjection = poSrcDS->GetGCPProjection();
    1126             :         if (BSBIsSRSOK(pszGCPProjection))
    1127             :         {
    1128             :             const GDAL_GCP *pasGCPList = poSrcDS->GetGCPs();
    1129             :             for (int i = 0; i < nGCPCount; i++)
    1130             :             {
    1131             :                 VSIFPrintfL(psBSB->fp, "REF/%d,%f,%f,%f,%f\n", i + 1,
    1132             :                             pasGCPList[i].dfGCPPixel, pasGCPList[i].dfGCPLine,
    1133             :                             pasGCPList[i].dfGCPY, pasGCPList[i].dfGCPX);
    1134             :             }
    1135             :         }
    1136             :     }
    1137             :     else if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None)
    1138             :     {
    1139             :         const char *pszProjection = poSrcDS->GetProjectionRef();
    1140             :         if (BSBIsSRSOK(pszProjection))
    1141             :         {
    1142             :             VSIFPrintfL(psBSB->fp, "REF/%d,%d,%d,%f,%f\n", 1, 0, 0,
    1143             :                         adfGeoTransform[3] + 0 * adfGeoTransform[4] +
    1144             :                             0 * adfGeoTransform[5],
    1145             :                         adfGeoTransform[0] + 0 * adfGeoTransform[1] +
    1146             :                             0 * adfGeoTransform[2]);
    1147             :             VSIFPrintfL(psBSB->fp, "REF/%d,%d,%d,%f,%f\n", 2, nXSize, 0,
    1148             :                         adfGeoTransform[3] + nXSize * adfGeoTransform[4] +
    1149             :                             0 * adfGeoTransform[5],
    1150             :                         adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
    1151             :                             0 * adfGeoTransform[2]);
    1152             :             VSIFPrintfL(psBSB->fp, "REF/%d,%d,%d,%f,%f\n", 3, nXSize, nYSize,
    1153             :                         adfGeoTransform[3] + nXSize * adfGeoTransform[4] +
    1154             :                             nYSize * adfGeoTransform[5],
    1155             :                         adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
    1156             :                             nYSize * adfGeoTransform[2]);
    1157             :             VSIFPrintfL(psBSB->fp, "REF/%d,%d,%d,%f,%f\n", 4, 0, nYSize,
    1158             :                         adfGeoTransform[3] + 0 * adfGeoTransform[4] +
    1159             :                             nYSize * adfGeoTransform[5],
    1160             :                         adfGeoTransform[0] + 0 * adfGeoTransform[1] +
    1161             :                             nYSize * adfGeoTransform[2]);
    1162             :         }
    1163             :     }
    1164             : 
    1165             :     /* -------------------------------------------------------------------- */
    1166             :     /*      Loop over image, copying image data.                            */
    1167             :     /* -------------------------------------------------------------------- */
    1168             :     CPLErr eErr = CE_None;
    1169             : 
    1170             :     GByte *pabyScanline = (GByte *)CPLMalloc(nXSize);
    1171             : 
    1172             :     for (int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++)
    1173             :     {
    1174             :         eErr =
    1175             :             poBand->RasterIO(GF_Read, 0, iLine, nXSize, 1, pabyScanline, nXSize,
    1176             :                              1, GDT_Byte, nBands, nBands * nXSize, NULL);
    1177             :         if (eErr == CE_None)
    1178             :         {
    1179             :             for (int i = 0; i < nXSize; i++)
    1180             :                 pabyScanline[i] = (GByte)anRemap[pabyScanline[i]];
    1181             : 
    1182             :             if (!BSBWriteScanline(psBSB, pabyScanline))
    1183             :                 eErr = CE_Failure;
    1184             :         }
    1185             :     }
    1186             : 
    1187             :     CPLFree(pabyScanline);
    1188             : 
    1189             :     /* -------------------------------------------------------------------- */
    1190             :     /*      cleanup                                                         */
    1191             :     /* -------------------------------------------------------------------- */
    1192             :     BSBClose(psBSB);
    1193             : 
    1194             :     if (eErr != CE_None)
    1195             :     {
    1196             :         VSIUnlink(pszFilename);
    1197             :         return NULL;
    1198             :     }
    1199             : 
    1200             :     return (GDALDataset *)GDALOpen(pszFilename, GA_ReadOnly);
    1201             : }
    1202             : #endif
    1203             : 
    1204             : /************************************************************************/
    1205             : /*                        GDALRegister_BSB()                            */
    1206             : /************************************************************************/
    1207             : 
    1208        1686 : void GDALRegister_BSB()
    1209             : 
    1210             : {
    1211        1686 :     if (GDALGetDriverByName("BSB") != nullptr)
    1212         302 :         return;
    1213             : 
    1214        1384 :     GDALDriver *poDriver = new GDALDriver();
    1215             : 
    1216        1384 :     poDriver->SetDescription("BSB");
    1217        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1218        1384 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Maptech BSB Nautical Charts");
    1219        1384 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/bsb.html");
    1220        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1221        1384 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "kap");
    1222             : #ifdef BSB_CREATE
    1223             :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte");
    1224             : #endif
    1225        1384 :     poDriver->pfnOpen = BSBDataset::Open;
    1226        1384 :     poDriver->pfnIdentify = BSBDataset::Identify;
    1227             : #ifdef BSB_CREATE
    1228             :     poDriver->pfnCreateCopy = BSBCreateCopy;
    1229             : #endif
    1230             : 
    1231        1384 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1232             : }

Generated by: LCOV version 1.14