text:     /******************************************************************************
   *
   * Project:  BSB Reader
   * Purpose:  BSBDataset implementation for BSB format.
   * Author:   Frank Warmerdam, warmerdam@pobox.com
   *
   ******************************************************************************
   * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
   * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
   *
   * Permission is hereby granted, free of charge, to any person obtaining a
   * copy of this software and associated documentation files (the "Software"),
   * to deal in the Software without restriction, including without limitation
   * the rights to use, copy, modify, merge, publish, distribute, sublicense,
   * and/or sell copies of the Software, and to permit persons to whom the
   * Software is furnished to do so, subject to the following conditions:
   *
   * The above copyright notice and this permission notice shall be included
   * in all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
   * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
   * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
   * DEALINGS IN THE SOFTWARE.
   ****************************************************************************/
  
  #include "bsb_read.h"
  #include "cpl_string.h"
  #include "gdal_frmts.h"
  #include "gdal_pam.h"
  #include "ogr_spatialref.h"
  
  #include <cstdlib>
  #include <algorithm>
  
  CPL_CVSID("$Id: bsbdataset.cpp  $")
  
  //Disabled as people may worry about the BSB patent
  //#define BSB_CREATE
  
  /************************************************************************/
  /* ==================================================================== */
  /*                              BSBDataset                              */
  /* ==================================================================== */
  /************************************************************************/
  
  class BSBRasterBand;
  
  class BSBDataset final: public GDALPamDataset
  {
      int         nGCPCount;
      GDAL_GCP    *pasGCPList;
      CPLString   osGCPProjection;
  
      double      adfGeoTransform[6];
      int         bGeoTransformSet;
  
      void        ScanForGCPs( bool isNos, const char *pszFilename );
      void        ScanForGCPsNos( const char *pszFilename );
      void        ScanForGCPsBSB();
  
      void        ScanForCutline();
  
      static int IdentifyInternal( GDALOpenInfo *, bool & isNosOut );
  
    public:
      BSBDataset();
      ~BSBDataset() override;
  
      BSBInfo     *psInfo;
  
      static GDALDataset *Open( GDALOpenInfo * );
      static int Identify( GDALOpenInfo * );
  
      int GetGCPCount() override;
      const char *_GetGCPProjection() override;
      const OGRSpatialReference* GetSpatialRef() const override {
          return GetSpatialRefFromOldGetProjectionRef();
      }
      const GDAL_GCP *GetGCPs() override;
  
      CPLErr GetGeoTransform( double * padfTransform ) override;
      const char *_GetProjectionRef() override;
      const OGRSpatialReference* GetGCPSpatialRef() const override {
          return GetGCPSpatialRefFromOldGetGCPProjection();
      }
  };
  
  /************************************************************************/
  /* ==================================================================== */
  /*                            BSBRasterBand                             */
  /* ==================================================================== */
  /************************************************************************/
  
  class BSBRasterBand final: public GDALPamRasterBand
  {
      GDALColorTable      oCT;
  
    public:
      explicit    BSBRasterBand( BSBDataset * );
  
      CPLErr IReadBlock( int, int, void * ) override;
      GDALColorTable *GetColorTable() override;
      GDALColorInterp GetColorInterpretation() override;
  };
  
  /************************************************************************/
  /*                           BSBRasterBand()                            */
  /************************************************************************/
  
  BSBRasterBand::BSBRasterBand( BSBDataset *poDSIn )
  
  {
      poDS = poDSIn;
      nBand = 1;
  
      eDataType = GDT_Byte;
  
      nBlockXSize = poDS->GetRasterXSize();
      nBlockYSize = 1;
  
      // Note that the first color table entry is dropped, everything is
      // shifted down.
      for( int i = 0; i < poDSIn->psInfo->nPCTSize-1; i++ )
      {
          GDALColorEntry oColor = {
              poDSIn->psInfo->pabyPCT[i*3+0+3],
              poDSIn->psInfo->pabyPCT[i*3+1+3],
              poDSIn->psInfo->pabyPCT[i*3+2+3],
              255
          };
  
          oCT.SetColorEntry( i, &oColor );
      }
  }
  
  /************************************************************************/
  /*                             IReadBlock()                             */
  /************************************************************************/
  
  CPLErr BSBRasterBand::IReadBlock( CPL_UNUSED int nBlockXOff,
                                    int nBlockYOff,
                                    void * pImage )
  {
      BSBDataset *poGDS = (BSBDataset *) poDS;
      GByte *pabyScanline = (GByte*) pImage;
  
      if( BSBReadScanline( poGDS->psInfo, nBlockYOff, pabyScanline ) )
      {
          for( int i = 0; i < nBlockXSize; i++ )
          {
              /* The indices start at 1, except in case of some charts */
              /* where there are missing values, which are filled to 0 */
              /* by BSBReadScanline */
              if (pabyScanline[i] > 0)
                  pabyScanline[i] -= 1;
          }
  
          return CE_None;
      }
  
      return CE_Failure;
  }
  
  /************************************************************************/
  /*                           GetColorTable()                            */
  /************************************************************************/
  
  GDALColorTable *BSBRasterBand::GetColorTable()
  
  {
      return &oCT;
  }
  
  /************************************************************************/
  /*                       GetColorInterpretation()                       */
  /************************************************************************/
  
  GDALColorInterp BSBRasterBand::GetColorInterpretation()
  
  {
      return GCI_PaletteIndex;
  }
  
  /************************************************************************/
  /* ==================================================================== */
  /*                              BSBDataset                              */
  /* ==================================================================== */
  /************************************************************************/
  
  /************************************************************************/
  /*                           BSBDataset()                               */
  /************************************************************************/
  
  BSBDataset::BSBDataset() :
      nGCPCount(0),
      pasGCPList(nullptr),
      osGCPProjection(SRS_WKT_WGS84_LAT_LONG),
      bGeoTransformSet(FALSE),
      psInfo(nullptr)
  {
      adfGeoTransform[0] = 0.0;     /* X Origin (top left corner) */
      adfGeoTransform[1] = 1.0;     /* X Pixel size */
      adfGeoTransform[2] = 0.0;
      adfGeoTransform[3] = 0.0;     /* Y Origin (top left corner) */
      adfGeoTransform[4] = 0.0;
      adfGeoTransform[5] = 1.0;     /* Y Pixel Size */
  }
  
  /************************************************************************/
  /*                            ~BSBDataset()                             */
  /************************************************************************/
  
  BSBDataset::~BSBDataset()
  
  {
      FlushCache(true);
  
      GDALDeinitGCPs( nGCPCount, pasGCPList );
      CPLFree( pasGCPList );
  
      if( psInfo != nullptr )
          BSBClose( psInfo );
  }
  
  /************************************************************************/
  /*                          GetGeoTransform()                           */
  /************************************************************************/
  
  CPLErr BSBDataset::GetGeoTransform( double * padfTransform )
  
  {
      memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
  
      if( bGeoTransformSet )
          return CE_None;
  
      return CE_Failure;
  }
  
  /************************************************************************/
  /*                          GetProjectionRef()                          */
  /************************************************************************/
  
  const char *BSBDataset::_GetProjectionRef()
  
  {
      if( bGeoTransformSet )
          return osGCPProjection;
  
      return "";
  }
  
  /************************************************************************/
  /*                     GDALHeuristicDatelineWrap()                      */
  /************************************************************************/
  
  static void
  GDALHeuristicDatelineWrap( int nPointCount, double *padfX )
  
  {
      if( nPointCount < 2 )
          return;
  
  /* -------------------------------------------------------------------- */
  /*      Work out what the longitude range will be centering on the      */
  /*      prime meridian (-180 to 180) and centering on the dateline      */
  /*      (0 to 360).                                                     */
  /* -------------------------------------------------------------------- */
      /* Following inits are useless but keep GCC happy */
      double dfX_PM_Min = 0.0;
      double dfX_PM_Max = 0.0;
      double dfX_Dateline_Min = 0.0;
      double dfX_Dateline_Max = 0.0;
  
      for( int i = 0; i < nPointCount; i++ )
      {
          double dfX_PM = padfX[i];
          if( dfX_PM > 180 )
              dfX_PM -= 360.0;
  
          double dfX_Dateline = padfX[i];
          if( dfX_Dateline < 0 )
              dfX_Dateline += 360.0;
  
          if( i == 0 )
          {
              dfX_PM_Min = dfX_PM;
              dfX_PM_Max = dfX_PM;
              dfX_Dateline_Min = dfX_Dateline;
              dfX_Dateline_Max = dfX_Dateline;
          }
          else
          {
              dfX_PM_Min = std::min(dfX_PM_Min, dfX_PM);
              dfX_PM_Max = std::max(dfX_PM_Max, dfX_PM);
              dfX_Dateline_Min = std::min(dfX_Dateline_Min, dfX_Dateline);
              dfX_Dateline_Max = std::max(dfX_Dateline_Max, dfX_Dateline);
          }
      }
  
  /* -------------------------------------------------------------------- */
  /*      Do nothing if the range is always fairly small - no apparent    */
  /*      wrapping issues.                                                */
  /* -------------------------------------------------------------------- */
      if( (dfX_PM_Max - dfX_PM_Min) < 270.0
          && (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0 )
          return;
  
  /* -------------------------------------------------------------------- */
  /*      Do nothing if both approach have a wide range - best not to    */
  /*      fiddle if we aren't sure we are improving things.               */
  /* -------------------------------------------------------------------- */
      if( (dfX_PM_Max - dfX_PM_Min) > 270.0
          && (dfX_Dateline_Max - dfX_Dateline_Min) > 270.0 )
          return;
  
  /* -------------------------------------------------------------------- */
  /*      Pick which way to transform things.                             */
  /* -------------------------------------------------------------------- */
      bool bUsePMWrap;
  
      if( (dfX_PM_Max - dfX_PM_Min) > 270.0
          && (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0 )
          bUsePMWrap = false;
      else
          bUsePMWrap = true;
  
  /* -------------------------------------------------------------------- */
  /*      Apply rewrapping.                                               */
  /* -------------------------------------------------------------------- */
      for( int i = 0; i < nPointCount; i++ )
      {
          if( bUsePMWrap )
          {
              if( padfX[i] > 180 )
                  padfX[i] -= 360.0;
          }
          else
          {
              if( padfX[i] < 0 )
                  padfX[i] += 360.0;
          }
      }
  }
  
  /************************************************************************/
  /*                   GDALHeuristicDatelineWrapGCPs()                    */
  /************************************************************************/
  
  static void
  GDALHeuristicDatelineWrapGCPs( int nPointCount, GDAL_GCP *pasGCPList )
  {
      std::vector<double> oadfX;
  
      oadfX.resize( nPointCount );
      for( int i = 0; i < nPointCount; i++ )
          oadfX[i] = pasGCPList[i].dfGCPX;
  
      GDALHeuristicDatelineWrap( nPointCount, &(oadfX[0]) );
  
      for( int i = 0; i < nPointCount; i++ )
          pasGCPList[i].dfGCPX = oadfX[i];
  }
  
  /************************************************************************/
  /*                            ScanForGCPs()                             */
  /************************************************************************/
  
  void BSBDataset::ScanForGCPs( bool isNos, const char *pszFilename )
  
  {
  /* -------------------------------------------------------------------- */
  /*      Collect GCPs as appropriate to source.                          */
  /* -------------------------------------------------------------------- */
      nGCPCount = 0;
  
      if ( isNos )
      {
          ScanForGCPsNos(pszFilename);
      } else {
          ScanForGCPsBSB();
      }
  
  /* -------------------------------------------------------------------- */
  /*      Apply heuristics to re-wrap GCPs to maintain continuity        */
  /*      over the international dateline.                                */
  /* -------------------------------------------------------------------- */
      if( nGCPCount > 1 )
          GDALHeuristicDatelineWrapGCPs( nGCPCount, pasGCPList );
  
  /* -------------------------------------------------------------------- */
  /*      Collect coordinate system related parameters from header.       */
  /* -------------------------------------------------------------------- */
      const char *pszKNP=nullptr;
      const char *pszKNQ=nullptr;
  
      for( int i = 0; psInfo->papszHeader[i] != nullptr; i++ )
      {
          if( STARTS_WITH_CI(psInfo->papszHeader[i], "KNP/") )
          {
              pszKNP = psInfo->papszHeader[i];
              SetMetadataItem( "BSB_KNP", pszKNP + 4 );
          }
          if( STARTS_WITH_CI(psInfo->papszHeader[i], "KNQ/") )
          {
              pszKNQ = psInfo->papszHeader[i];
              SetMetadataItem( "BSB_KNQ", pszKNQ + 4 );
          }
      }
  
  /* -------------------------------------------------------------------- */
  /*      Can we derive a reasonable coordinate system definition for     */
  /*      this file?  For now we keep it simple, just handling            */
  /*      mercator. In the future we should consider others.              */
  /* -------------------------------------------------------------------- */
      CPLString osUnderlyingSRS;
      if( pszKNP != nullptr )
      {
          const char *pszPR = strstr(pszKNP,"PR=");
          const char *pszGD = strstr(pszKNP,"GD=");
          const char *pszGEOGCS = SRS_WKT_WGS84_LAT_LONG;
          CPLString osPP;
  
          // Capture the PP string.
          const char *pszValue = strstr(pszKNP,"PP=");
          const char *pszEnd = pszValue ? strstr(pszValue,",") : nullptr;
          if( pszValue && pszEnd )
              osPP.assign(pszValue+3,pszEnd-pszValue-3);
  
          // Look at the datum
          if( pszGD == nullptr )
          {
              /* no match. We'll default to EPSG:4326 */
          }
          else if( STARTS_WITH_CI(pszGD, "GD=European 1950") )
          {
              pszGEOGCS = "GEOGCS[\"ED50\",DATUM[\"European_Datum_1950\",SPHEROID[\"International 1924\",6378388,297,AUTHORITY[\"EPSG\",\"7022\"]],TOWGS84[-87,-98,-121,0,0,0,0],AUTHORITY[\"EPSG\",\"6230\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4230\"]]";
          }
  
          // Look at the projection
          if( pszPR == nullptr )
          {
              /* no match */
          }
          else if( STARTS_WITH_CI(pszPR, "PR=MERCATOR") && nGCPCount > 0 )
          {
              // We somewhat arbitrarily select our first GCPX as our
              // central meridian.  This is mostly helpful to ensure
              // that regions crossing the dateline will be contiguous
              // in mercator.
              osUnderlyingSRS.Printf( "PROJCS[\"Global Mercator\",%s,PROJECTION[\"Mercator_2SP\"],PARAMETER[\"standard_parallel_1\",0],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%d],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]",
                  pszGEOGCS, (int) pasGCPList[0].dfGCPX );
          }
  
          else if( STARTS_WITH_CI(pszPR, "PR=TRANSVERSE MERCATOR")
                   && !osPP.empty() )
          {
  
              osUnderlyingSRS.Printf(
                  "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%s],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]",
                  pszGEOGCS, osPP.c_str() );
          }
  
          else if( STARTS_WITH_CI(pszPR, "PR=UNIVERSAL TRANSVERSE MERCATOR")
                   && !osPP.empty() )
          {
              // This is not *really* UTM unless the central meridian
              // matches a zone which it does not in some (most?) maps.
              osUnderlyingSRS.Printf(
                  "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%s],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]",
                  pszGEOGCS, osPP.c_str() );
          }
  
          else if( STARTS_WITH_CI(pszPR, "PR=POLYCONIC") && !osPP.empty() )
          {
              osUnderlyingSRS.Printf(
                  "PROJCS[\"unnamed\",%s,PROJECTION[\"Polyconic\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%s],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]",
                  pszGEOGCS, osPP.c_str() );
          }
  
          else if( STARTS_WITH_CI(pszPR, "PR=LAMBERT CONFORMAL CONIC")
                   && !osPP.empty() && pszKNQ != nullptr )
          {
              CPLString osP2, osP3;
  
              // Capture the KNQ/P2 string.
              pszValue = strstr(pszKNQ,"P2=");
              if( pszValue )
                  pszEnd = strstr(pszValue,",");
              if( pszValue && pszEnd )
                  osP2.assign(pszValue+3,pszEnd-pszValue-3);
  
              // Capture the KNQ/P3 string.
              pszValue = strstr(pszKNQ,"P3=");
              if( pszValue )
                  pszEnd = strstr(pszValue,",");
              if( pszValue )
              {
                  if( pszEnd )
                      osP3.assign(pszValue+3,pszEnd-pszValue-3);
                  else
                      osP3.assign(pszValue+3);
              }
  
              if( !osP2.empty() && !osP3.empty() )
                  osUnderlyingSRS.Printf(
                      "PROJCS[\"unnamed\",%s,PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",%s],PARAMETER[\"standard_parallel_2\",%s],PARAMETER[\"latitude_of_origin\",0.0],PARAMETER[\"central_meridian\",%s],PARAMETER[\"false_easting\",0.0],PARAMETER[\"false_northing\",0.0],UNIT[\"Meter\",1]]",
                      pszGEOGCS, osP2.c_str(), osP3.c_str(), osPP.c_str() );
          }
      }
  
  /* -------------------------------------------------------------------- */
  /*      If we got an alternate underlying coordinate system, try        */
  /*      converting the GCPs to that coordinate system.                  */
  /* -------------------------------------------------------------------- */
      if( osUnderlyingSRS.length() > 0 )
      {
          OGRSpatialReference oGeog_SRS, oProjected_SRS;
  
          oProjected_SRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
          oProjected_SRS.SetFromUserInput( osUnderlyingSRS );
          oGeog_SRS.CopyGeogCSFrom( &oProjected_SRS );
          oGeog_SRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
  
          OGRCoordinateTransformation *poCT
              = OGRCreateCoordinateTransformation( &oGeog_SRS,
                                                   &oProjected_SRS );
          if( poCT != nullptr )
          {
              for( int i = 0; i < nGCPCount; i++ )
              {
                  poCT->Transform( 1,
                                   &(pasGCPList[i].dfGCPX),
                                   &(pasGCPList[i].dfGCPY),
                                   &(pasGCPList[i].dfGCPZ) );
              }
  
              osGCPProjection = osUnderlyingSRS;
  
              delete poCT;
          }
          else
              CPLErrorReset();
      }
  
  /* -------------------------------------------------------------------- */
  /*      Attempt to prepare a geotransform from the GCPs.                */
  /* -------------------------------------------------------------------- */
      if( GDALGCPsToGeoTransform( nGCPCount, pasGCPList, adfGeoTransform,
                                  FALSE ) )
      {
          bGeoTransformSet = TRUE;
      }
  }
  
  /************************************************************************/
  /*                           ScanForGCPsNos()                           */
  /*                                                                      */
  /*      Nos files have an accompanying .geo file, that contains some    */
  /*      of the information normally contained in the header section     */
  /*      with BSB files. we try and open a file with the same name,      */
  /*      but a .geo extension, and look for lines like...                */
  /*      PointX=long lat line pixel    (using the same naming system     */
  /*      as BSB) Point1=-22.0000 64.250000 197 744                       */
  /************************************************************************/
  
  void BSBDataset::ScanForGCPsNos( const char *pszFilename )
  {
      const char *extension = CPLGetExtension(pszFilename);
  
      // pseudointelligently try and guess whether we want a .geo or a .GEO
      const char *geofile = nullptr;
      if (extension[1] == 'O')
      {
          geofile = CPLResetExtension( pszFilename, "GEO");
      } else {
          geofile = CPLResetExtension( pszFilename, "geo");
      }
  
      FILE *gfp = VSIFOpen( geofile, "r" );  // Text files
      if( gfp == nullptr )
      {
          CPLError( CE_Failure, CPLE_OpenFailed,
                    "Couldn't find a matching .GEO file: %s", geofile );
          return;
      }
  
      char *thisLine = (char *) CPLMalloc( 80 ); // FIXME
  
      // Count the GCPs (reference points) and seek the file pointer 'gfp' to the starting point
      int fileGCPCount=0;
      while (fgets(thisLine, 80, gfp))
      {
          if( STARTS_WITH_CI(thisLine, "Point") )
              fileGCPCount++;
      }
      VSIRewind( gfp );
  
      // Memory has not been allocated to fileGCPCount yet
      pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),fileGCPCount+1);
  
      while (fgets(thisLine, 80, gfp))
      {
          if( STARTS_WITH_CI(thisLine, "Point") )
          {
              // got a point line, turn it into a gcp
              char **Tokens
                  = CSLTokenizeStringComplex(thisLine, "= ", FALSE, FALSE);
              if (CSLCount(Tokens) >= 5)
              {
                  GDALInitGCPs( 1, pasGCPList + nGCPCount );
                  pasGCPList[nGCPCount].dfGCPX = CPLAtof(Tokens[1]);
                  pasGCPList[nGCPCount].dfGCPY = CPLAtof(Tokens[2]);
                  pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(Tokens[4]);
                  pasGCPList[nGCPCount].dfGCPLine = CPLAtof(Tokens[3]);
  
                  CPLFree( pasGCPList[nGCPCount].pszId );
                  char szName[50];
                  snprintf( szName, sizeof(szName), "GCP_%d", nGCPCount+1 );
                  pasGCPList[nGCPCount].pszId = CPLStrdup( szName );
  
                  nGCPCount++;
              }
              CSLDestroy(Tokens);
          }
      }
  
      CPLFree(thisLine);
      VSIFClose(gfp);
  }
  
  /************************************************************************/
  /*                            ScanForGCPsBSB()                          */
  /************************************************************************/
  
  void BSBDataset::ScanForGCPsBSB()
  {
  /* -------------------------------------------------------------------- */
  /*      Collect standalone GCPs.  They look like:                       */
  /*                                                                      */
  /*      REF/1,115,2727,32.346666666667,-60.881666666667                 */
  /*      REF/n,pixel,line,lat,long                                       */
  /* -------------------------------------------------------------------- */
      int fileGCPCount=0;
  
      // Count the GCPs (reference points) in psInfo->papszHeader
      for( int i = 0; psInfo->papszHeader[i] != nullptr; i++ )
          if( STARTS_WITH_CI(psInfo->papszHeader[i], "REF/") )
              fileGCPCount++;
  
      // Memory has not been allocated to fileGCPCount yet
      pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),fileGCPCount+1);
  
      for( int i = 0; psInfo->papszHeader[i] != nullptr; i++ )
      {
  
          if( !STARTS_WITH_CI(psInfo->papszHeader[i], "REF/") )
              continue;
  
          char **papszTokens
              = CSLTokenizeStringComplex( psInfo->papszHeader[i]+4, ",",
                                          FALSE, FALSE );
  
          if( CSLCount(papszTokens) > 4 )
          {
              GDALInitGCPs( 1, pasGCPList + nGCPCount );
  
              pasGCPList[nGCPCount].dfGCPX = CPLAtof(papszTokens[4]);
              pasGCPList[nGCPCount].dfGCPY = CPLAtof(papszTokens[3]);
              pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(papszTokens[1]);
              pasGCPList[nGCPCount].dfGCPLine = CPLAtof(papszTokens[2]);
  
              CPLFree( pasGCPList[nGCPCount].pszId );
              if( CSLCount(papszTokens) > 5 )
              {
                  pasGCPList[nGCPCount].pszId = CPLStrdup(papszTokens[5]);
              }
              else
              {
                  char szName[50];
                  snprintf( szName, sizeof(szName), "GCP_%d", nGCPCount+1 );
                  pasGCPList[nGCPCount].pszId = CPLStrdup( szName );
              }
  
              nGCPCount++;
          }
          CSLDestroy( papszTokens );
      }
  }
  
  /************************************************************************/
  /*                            ScanForCutline()                          */
  /************************************************************************/
  
  void BSBDataset::ScanForCutline()
  {
      /* PLY: Border Polygon Record - coordinates of the panel within the 
      * raster image, given in chart datum lat/long. Any shape polygon.
      * They look like:
      *      PLY/1,32.346666666667,-60.881666666667 
      *      PLY/n,lat,long 
      *
      * If found then we return it via a BSB_CUTLINE metadata item as a WKT POLYGON.
      */
  
      std::string wkt;
      for( int i = 0; psInfo->papszHeader[i] != nullptr; i++ )
      {
          if( !STARTS_WITH_CI(psInfo->papszHeader[i], "PLY/") )
              continue;
  
          const CPLStringList aosTokens(
              CSLTokenizeString2( psInfo->papszHeader[i]+4, ",", 0 ));
  
          if( aosTokens.size() >= 3 )
          {
              if (wkt.empty())
                  wkt = "POLYGON ((";
              else
                  wkt += ',';
              wkt += aosTokens[2];
              wkt += ' ';
              wkt += aosTokens[1];
          }
      }
  
      if (!wkt.empty())
      {
          wkt += "))";
          SetMetadataItem("BSB_CUTLINE", wkt.c_str());
      }
  }
  
  /************************************************************************/
  /*                          IdentifyInternal()                          */
  /************************************************************************/
  
  int BSBDataset::IdentifyInternal( GDALOpenInfo * poOpenInfo, bool& isNosOut )
  
  {
  /* -------------------------------------------------------------------- */
  /*      Check for BSB/ keyword.                                         */
  /* -------------------------------------------------------------------- */
      isNosOut = false;
  
      if( poOpenInfo->nHeaderBytes < 1000 )
          return FALSE;
  
      int i = 0;
      for( ; i < poOpenInfo->nHeaderBytes - 4; i++ )
      {
          if( poOpenInfo->pabyHeader[i+0] == 'B'
              && poOpenInfo->pabyHeader[i+1] == 'S'
              && poOpenInfo->pabyHeader[i+2] == 'B'
              && poOpenInfo->pabyHeader[i+3] == '/' )
              break;
          if( poOpenInfo->pabyHeader[i+0] == 'N'
              && poOpenInfo->pabyHeader[i+1] == 'O'
              && poOpenInfo->pabyHeader[i+2] == 'S'
              && poOpenInfo->pabyHeader[i+3] == '/' )
          {
              isNosOut = true;
              break;
          }
          if( poOpenInfo->pabyHeader[i+0] == 'W'
              && poOpenInfo->pabyHeader[i+1] == 'X'
              && poOpenInfo->pabyHeader[i+2] == '\\'
              && poOpenInfo->pabyHeader[i+3] == '8' )
              break;
      }
  
      if( i == poOpenInfo->nHeaderBytes - 4 )
          return FALSE;
  
      /* Additional test to avoid false positive. See #2881 */
      const char* pszRA = strstr((const char*)poOpenInfo->pabyHeader + i, "RA=");
      if (pszRA == nullptr) /* This may be a NO1 file */
          pszRA = strstr((const char*)poOpenInfo->pabyHeader + i, "[JF");
      if (pszRA == nullptr || pszRA - ((const char*)poOpenInfo->pabyHeader + i) > 100 )
          return FALSE;
  
      return TRUE;
  }
  
  /************************************************************************/
  /*                              Identify()                              */
  /************************************************************************/
  
  int BSBDataset::Identify( GDALOpenInfo * poOpenInfo )
  
  {
      bool isNos;
      return IdentifyInternal(poOpenInfo, isNos);
  }
  
  /************************************************************************/
  /*                                Open()                                */
  /************************************************************************/
  
  GDALDataset *BSBDataset::Open( GDALOpenInfo * poOpenInfo )
  
  {
      bool isNos = false;
      if (!IdentifyInternal(poOpenInfo, isNos))
          return nullptr;
  
      if( poOpenInfo->eAccess == GA_Update )
      {
          CPLError( CE_Failure, CPLE_NotSupported,
                    "The BSB driver does not support update access to existing"
                    " datasets.\n" );
          return nullptr;
      }
  
  /* -------------------------------------------------------------------- */
  /*      Create a corresponding GDALDataset.                             */
  /* -------------------------------------------------------------------- */
      BSBDataset *poDS = new BSBDataset();
  
  /* -------------------------------------------------------------------- */
  /*      Open the file.                                                  */
  /* -------------------------------------------------------------------- */
      poDS->psInfo = BSBOpen( poOpenInfo->pszFilename );
      if( poDS->psInfo == nullptr )
      {
          delete poDS;
          return nullptr;
      }
  
      poDS->nRasterXSize = poDS->psInfo->nXSize;
      poDS->nRasterYSize = poDS->psInfo->nYSize;
  
  /* -------------------------------------------------------------------- */
  /*      Create band information objects.                                */
  /* -------------------------------------------------------------------- */
      poDS->SetBand( 1, new BSBRasterBand( poDS ));
  
      poDS->ScanForGCPs( isNos, poOpenInfo->pszFilename );
      
  /* -------------------------------------------------------------------- */
  /*      Set CUTLINE metadata if a bounding polygon is available         */
  /* -------------------------------------------------------------------- */
      poDS->ScanForCutline();
  
  /* -------------------------------------------------------------------- */
  /*      Initialize any PAM information.                                 */
  /* -------------------------------------------------------------------- */
      poDS->SetDescription( poOpenInfo->pszFilename );
      poDS->TryLoadXML();
  
      poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
  
      return poDS;
  }
  
  /************************************************************************/
  /*                            GetGCPCount()                             */
  /************************************************************************/
  
  int BSBDataset::GetGCPCount()
  
  {
      return nGCPCount;
  }
  
  /************************************************************************/
  /*                          GetGCPProjection()                          */
  /************************************************************************/
  
  const char *BSBDataset::_GetGCPProjection()
  
  {
      return osGCPProjection;
  }
  
  /************************************************************************/
  /*                               GetGCP()                               */
  /************************************************************************/
  
  const GDAL_GCP *BSBDataset::GetGCPs()
  
  {
      return pasGCPList;
  }
  
  #ifdef BSB_CREATE
  
  /************************************************************************/
  /*                             BSBIsSRSOK()                             */
  /************************************************************************/
  
  static int BSBIsSRSOK(const char *pszWKT)
  {
      bool bOK = false;
      OGRSpatialReference oSRS, oSRS_WGS84, oSRS_NAD83;
  
      if( pszWKT != NULL && pszWKT[0] != '\0' )
      {
          oSRS.importFromWkt( pszWKT );
  
          oSRS_WGS84.SetWellKnownGeogCS( "WGS84" );
          oSRS_NAD83.SetWellKnownGeogCS( "NAD83" );
          if ( (oSRS.IsSameGeogCS(&oSRS_WGS84) || oSRS.IsSameGeogCS(&oSRS_NAD83)) &&
                oSRS.IsGeographic() && oSRS.GetPrimeMeridian() == 0.0 )
          {
              bOK = true;
          }
      }
  
      if (!bOK)
      {
          CPLError(CE_Warning, CPLE_NotSupported,
                  "BSB only supports WGS84 or NAD83 geographic projections.\n");
      }
  
      return bOK;
  }
  
  /************************************************************************/
  /*                           BSBCreateCopy()                            */
  /************************************************************************/
  
  static GDALDataset *
  BSBCreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
                 int bStrict, char ** /*papszOptions*/,
                 GDALProgressFunc /*pfnProgress*/, void * /*pProgressData*/ )
  
  {
  /* -------------------------------------------------------------------- */
  /*      Some some rudimentary checks                                    */
  /* -------------------------------------------------------------------- */
      const int nBands = poSrcDS->GetRasterCount();
      if( nBands != 1 )
      {
          CPLError( CE_Failure, CPLE_NotSupported,
                    "BSB driver only supports one band images.\n" );
  
          return NULL;
      }
  
      if( poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte
          && bStrict )
      {
          CPLError( CE_Failure, CPLE_NotSupported,
                    "BSB driver doesn't support data type %s. "
                    "Only eight bit bands supported.\n",
                    GDALGetDataTypeName(
                        poSrcDS->GetRasterBand(1)->GetRasterDataType()) );
  
          return NULL;
      }
  
      const int nXSize = poSrcDS->GetRasterXSize();
      const int nYSize = poSrcDS->GetRasterYSize();
  
  /* -------------------------------------------------------------------- */
  /*      Open the output file.                                           */
  /* -------------------------------------------------------------------- */
      BSBInfo *psBSB = BSBCreate( pszFilename, 0, 200, nXSize, nYSize );
      if( psBSB == NULL )
          return NULL;
  
  /* -------------------------------------------------------------------- */
  /*      Prepare initial color table.colortable.                         */
  /* -------------------------------------------------------------------- */
      GDALRasterBand      *poBand = poSrcDS->GetRasterBand(1);
      unsigned char       abyPCT[771];
      int                 nPCTSize;
      int                 anRemap[256];
  
      abyPCT[0] = 0;
      abyPCT[1] = 0;
      abyPCT[2] = 0;
  
      if( poBand->GetColorTable() == NULL )
      {
          /* map greyscale down to 63 grey levels. */
          for( int iColor = 0; iColor < 256; iColor++ )
          {
              int nOutValue = (int) (iColor / 4.1) + 1;
  
              anRemap[iColor] = nOutValue;
              abyPCT[nOutValue*3 + 0] = (unsigned char) iColor;
              abyPCT[nOutValue*3 + 1] = (unsigned char) iColor;
              abyPCT[nOutValue*3 + 2] = (unsigned char) iColor;
          }
          nPCTSize = 64;
      }
      else
      {
          GDALColorTable *poCT = poBand->GetColorTable();
          int nColorTableSize = poCT->GetColorEntryCount();
          if (nColorTableSize > 255)
              nColorTableSize = 255;
  
          for( int iColor = 0; iColor < nColorTableSize; iColor++ )
          {
              GDALColorEntry sEntry;
  
              poCT->GetColorEntryAsRGB( iColor, &sEntry );
  
              anRemap[iColor] = iColor + 1;
              abyPCT[(iColor+1)*3 + 0] = (unsigned char) sEntry.c1;
              abyPCT[(iColor+1)*3 + 1] = (unsigned char) sEntry.c2;
              abyPCT[(iColor+1)*3 + 2] = (unsigned char) sEntry.c3;
          }
  
          nPCTSize = nColorTableSize + 1;
  
          // Add entries for pixel values which apparently will not occur.
          for( int iColor = nPCTSize; iColor < 256; iColor++ )
              anRemap[iColor] = 1;
      }
  
  /* -------------------------------------------------------------------- */
  /*      Boil out all duplicate entries.                                 */
  /* -------------------------------------------------------------------- */
      for( int i = 1; i < nPCTSize-1; i++ )
      {
          for( int j = i+1; j < nPCTSize; j++ )
          {
              if( abyPCT[i*3+0] == abyPCT[j*3+0]
                  && abyPCT[i*3+1] == abyPCT[j*3+1]
                  && abyPCT[i*3+2] == abyPCT[j*3+2] )
              {
                  nPCTSize--;
                  abyPCT[j*3+0] = abyPCT[nPCTSize*3+0];
                  abyPCT[j*3+1] = abyPCT[nPCTSize*3+1];
                  abyPCT[j*3+2] = abyPCT[nPCTSize*3+2];
  
                  for( int k = 0; k < 256; k++ )
                  {
                      // merge matching entries.
                      if( anRemap[k] == j )
                          anRemap[k] = i;
  
                      // shift the last PCT entry into the new hole.
                      if( anRemap[k] == nPCTSize )
                          anRemap[k] = j;
                  }
              }
          }
      }
  
  /* -------------------------------------------------------------------- */
  /*      Boil out all duplicate entries.                                 */
  /* -------------------------------------------------------------------- */
      if( nPCTSize > 128 )
      {
          CPLError( CE_Warning, CPLE_AppDefined,
                    "Having to merge color table entries to reduce %d real\n"
                    "color table entries down to 127 values.",
                    nPCTSize );
      }
  
      while( nPCTSize > 128 )
      {
          int nBestRange = 768;
          int iBestMatch1 = -1;
          int iBestMatch2 = -1;
  
          // Find the closest pair of color table entries.
  
          for( int i = 1; i < nPCTSize-1; i++ )
          {
              for( int j = i+1; j < nPCTSize; j++ )
              {
                  int nRange = std::abs(abyPCT[i*3+0] - abyPCT[j*3+0])
                      + std::abs(abyPCT[i*3+1] - abyPCT[j*3+1])
                      + std::abs(abyPCT[i*3+2] - abyPCT[j*3+2]);
  
                  if( nRange < nBestRange )
                  {
                      iBestMatch1 = i;
                      iBestMatch2 = j;
                      nBestRange = nRange;
                  }
              }
          }
  
          // Merge the second entry into the first.
          nPCTSize--;
          abyPCT[iBestMatch2*3+0] = abyPCT[nPCTSize*3+0];
          abyPCT[iBestMatch2*3+1] = abyPCT[nPCTSize*3+1];
          abyPCT[iBestMatch2*3+2] = abyPCT[nPCTSize*3+2];
  
          for( int i = 0; i < 256; i++ )
          {
              // merge matching entries.
              if( anRemap[i] == iBestMatch2 )
                  anRemap[i] = iBestMatch1;
  
              // shift the last PCT entry into the new hole.
              if( anRemap[i] == nPCTSize )
                  anRemap[i] = iBestMatch2;
          }
      }
  
  /* -------------------------------------------------------------------- */
  /*      Write the PCT.                                                  */
  /* -------------------------------------------------------------------- */
      if( !BSBWritePCT( psBSB, nPCTSize, abyPCT ) )
      {
          BSBClose( psBSB );
          return NULL;
      }
  
  /* -------------------------------------------------------------------- */
  /*      Write the GCPs.                                                 */
  /* -------------------------------------------------------------------- */
      double adfGeoTransform[6];
      int nGCPCount = poSrcDS->GetGCPCount();
      if (nGCPCount)
      {
          const char* pszGCPProjection = poSrcDS->GetGCPProjection();
          if ( BSBIsSRSOK(pszGCPProjection) )
          {
              const GDAL_GCP * pasGCPList = poSrcDS->GetGCPs();
              for( int i = 0; i < nGCPCount; i++ )
              {
                  VSIFPrintfL( psBSB->fp,
                              "REF/%d,%f,%f,%f,%f\n",
                              i+1,
                              pasGCPList[i].dfGCPPixel, pasGCPList[i].dfGCPLine,
                              pasGCPList[i].dfGCPY, pasGCPList[i].dfGCPX);
              }
          }
      }
      else if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None)
      {
          const char* pszProjection = poSrcDS->GetProjectionRef();
          if ( BSBIsSRSOK(pszProjection) )
          {
              VSIFPrintfL( psBSB->fp,
                          "REF/%d,%d,%d,%f,%f\n",
                          1,
                          0, 0,
                          adfGeoTransform[3] + 0 * adfGeoTransform[4] + 0 * adfGeoTransform[5],
                          adfGeoTransform[0] + 0 * adfGeoTransform[1] + 0 * adfGeoTransform[2]);
              VSIFPrintfL( psBSB->fp,
                          "REF/%d,%d,%d,%f,%f\n",
                          2,
                          nXSize, 0,
                          adfGeoTransform[3] + nXSize * adfGeoTransform[4] + 0 * adfGeoTransform[5],
                          adfGeoTransform[0] + nXSize * adfGeoTransform[1] + 0 * adfGeoTransform[2]);
              VSIFPrintfL( psBSB->fp,
                          "REF/%d,%d,%d,%f,%f\n",
                          3,
                          nXSize, nYSize,
                          adfGeoTransform[3] + nXSize * adfGeoTransform[4] + nYSize * adfGeoTransform[5],
                          adfGeoTransform[0] + nXSize * adfGeoTransform[1] + nYSize * adfGeoTransform[2]);
              VSIFPrintfL( psBSB->fp,
                          "REF/%d,%d,%d,%f,%f\n",
                          4,
                          0, nYSize,
                          adfGeoTransform[3] + 0 * adfGeoTransform[4] + nYSize * adfGeoTransform[5],
                          adfGeoTransform[0] + 0 * adfGeoTransform[1] + nYSize * adfGeoTransform[2]);
          }
      }
  
  /* -------------------------------------------------------------------- */
  /*      Loop over image, copying image data.                            */
  /* -------------------------------------------------------------------- */
      CPLErr      eErr = CE_None;
  
      GByte *pabyScanline = (GByte *) CPLMalloc( nXSize );
  
      for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ )
      {
          eErr = poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
                                   pabyScanline, nXSize, 1, GDT_Byte,
                                   nBands, nBands * nXSize, NULL );
          if( eErr == CE_None )
          {
              for( int i = 0; i < nXSize; i++ )
                  pabyScanline[i] = (GByte) anRemap[pabyScanline[i]];
  
              if( !BSBWriteScanline( psBSB, pabyScanline ) )
                  eErr = CE_Failure;
          }
      }
  
      CPLFree( pabyScanline );
  
  /* -------------------------------------------------------------------- */
  /*      cleanup                                                         */
  /* -------------------------------------------------------------------- */
      BSBClose( psBSB );
  
      if( eErr != CE_None )
      {
          VSIUnlink( pszFilename );
          return NULL;
      }
  
      return (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
  }
  #endif
  
  /************************************************************************/
  /*                        GDALRegister_BSB()                            */
  /************************************************************************/
  
  void GDALRegister_BSB()
  
  {
      if( GDALGetDriverByName( "BSB" ) != nullptr )
          return;
  
      GDALDriver *poDriver = new GDALDriver();
  
      poDriver->SetDescription( "BSB" );
      poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
      poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
                                 "Maptech BSB Nautical Charts" );
      poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/bsb.html" );
      poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
      poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "kap" );
  #ifdef BSB_CREATE
      poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Byte" );
  #endif
      poDriver->pfnOpen = BSBDataset::Open;
      poDriver->pfnIdentify = BSBDataset::Identify;
  #ifdef BSB_CREATE
      poDriver->pfnCreateCopy = BSBCreateCopy;
  #endif
  
      GetGDALDriverManager()->RegisterDriver( poDriver );
  }
  