/******************************************************************************
 *
 * Project:  BNA Translator
 * Purpose:  Implements OGRBNALayer class.
 * Author:   Even Rouault, even dot rouault at mines dash paris dot org
 *
 ******************************************************************************
 * Copyright (c) 2007-2010, Even Rouault <even dot rouault at mines-paris dot org>
 *
 * 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 "ogr_bna.h"
#include "cpl_conv.h"
#include "cpl_string.h"
#include "cpl_csv.h"
#include "ogr_p.h"

CPL_CVSID("$Id: ogrbnalayer.cpp 7a3b8689d54765d306ea9cc4c773b4c967db980f 2018-10-14 19:50:43 +0200 Even Rouault $")

/************************************************************************/
/*                            OGRBNALayer()                             */
/*                                                                      */
/*      Note that the OGRBNALayer assumes ownership of the passed       */
/*      file pointer.                                                   */
/************************************************************************/

OGRBNALayer::OGRBNALayer( const char *pszFilename,
                          const char* layerName,
                          BNAFeatureType bnaFeatureTypeIn,
                          OGRwkbGeometryType eLayerGeomType,
                          int bWriterIn,
                          OGRBNADataSource* poDSIn,
                          int nIDsIn) :
    poDS(poDSIn),
    bWriter(bWriterIn),
    nIDs(nIDsIn),
    eof(false),
    failed(false),
    curLine(0),
    nNextFID(0),
    nFeatures(0),
    partialIndexTable(TRUE),
    offsetAndLineFeaturesTable(nullptr)
{
    static const char* const iKnowHowToCount[]
        = { "Primary", "Secondary", "Third", "Fourth", "Fifth" };
    char tmp[32];

    poFeatureDefn = new OGRFeatureDefn(
        CPLSPrintf( "%s_%s",
                    CPLGetBasename( pszFilename ),
                    layerName ));
    poFeatureDefn->Reference();
    poFeatureDefn->SetGeomType( eLayerGeomType );
    SetDescription( poFeatureDefn->GetName() );
    bnaFeatureType = bnaFeatureTypeIn;

    if (! bWriter )
    {
        for(int i=0;i<nIDs;i++)
        {
            if (i < static_cast<int>(
                        sizeof(iKnowHowToCount)/sizeof(iKnowHowToCount[0]) ) )
            {
                snprintf(tmp, sizeof(tmp), "%s ID", iKnowHowToCount[i]);
                OGRFieldDefn oFieldID(tmp, OFTString );
                poFeatureDefn->AddFieldDefn( &oFieldID );
            }
            else
            {
                snprintf(tmp, sizeof(tmp), "%dth ID", i+1);
                OGRFieldDefn oFieldID(tmp, OFTString );
                poFeatureDefn->AddFieldDefn( &oFieldID );
            }
        }

        if (bnaFeatureType == BNA_ELLIPSE)
        {
            OGRFieldDefn oFieldMajorRadius( "Major radius", OFTReal );
            poFeatureDefn->AddFieldDefn( &oFieldMajorRadius );

            OGRFieldDefn oFieldMinorRadius( "Minor radius", OFTReal );
            poFeatureDefn->AddFieldDefn( &oFieldMinorRadius );
        }

        fpBNA = VSIFOpenL( pszFilename, "rb" );
        if( fpBNA == nullptr )
            return;
    }
    else
    {
        fpBNA = nullptr;
    }
}

/************************************************************************/
/*                            ~OGRBNALayer()                            */
/************************************************************************/

OGRBNALayer::~OGRBNALayer()

{
    poFeatureDefn->Release();

    CPLFree(offsetAndLineFeaturesTable);

    if (fpBNA)
        VSIFCloseL( fpBNA );
}

/************************************************************************/
/*                         SetFeatureIndexTable()                       */
/************************************************************************/
void OGRBNALayer::SetFeatureIndexTable(
    int nFeaturesIn,
    OffsetAndLine* offsetAndLineFeaturesTableIn,
    int partialIndexTableIn )
{
    nFeatures = nFeaturesIn;
    offsetAndLineFeaturesTable = offsetAndLineFeaturesTableIn;
    partialIndexTable = CPL_TO_BOOL(partialIndexTableIn);
}

/************************************************************************/
/*                            ResetReading()                            */
/************************************************************************/

void OGRBNALayer::ResetReading()

{
    if( fpBNA == nullptr )
        return;
    eof = false;
    failed = false;
    curLine = 0;
    nNextFID = 0;
    CPL_IGNORE_RET_VAL(VSIFSeekL( fpBNA, 0, SEEK_SET ));
}

/************************************************************************/
/*                           GetNextFeature()                           */
/************************************************************************/

OGRFeature *OGRBNALayer::GetNextFeature()
{
    if( failed || eof || fpBNA == nullptr )
        return nullptr;

    while( true )
    {
        int ok = FALSE;
        const int offset = static_cast<int>( VSIFTellL(fpBNA) );
        const int line = curLine;
        if (nNextFID < nFeatures)
        {
            if( VSIFSeekL( fpBNA, offsetAndLineFeaturesTable[nNextFID].offset,
                           SEEK_SET ) < 0 )
                return nullptr;
            curLine = offsetAndLineFeaturesTable[nNextFID].line;
        }
        BNARecord* record
            = BNA_GetNextRecord(fpBNA, &ok, &curLine, TRUE, bnaFeatureType);
        if (ok == FALSE)
        {
            BNA_FreeRecord(record);
            failed = true;
            return nullptr;
        }
        if (record == nullptr)
        {
            /* end of file */
            eof = true;

            /* and we have finally build the whole index table */
            partialIndexTable = false;
            return nullptr;
        }

        if (record->featureType == bnaFeatureType)
        {
            if (nNextFID >= nFeatures)
            {
                nFeatures++;
                offsetAndLineFeaturesTable =
                    static_cast<OffsetAndLine *>(
                        CPLRealloc( offsetAndLineFeaturesTable,
                                    nFeatures * sizeof(OffsetAndLine) ) );
                offsetAndLineFeaturesTable[nFeatures-1].offset = offset;
                offsetAndLineFeaturesTable[nFeatures-1].line = line;
            }

            OGRFeature *poFeature
                = BuildFeatureFromBNARecord(record, nNextFID++);

            BNA_FreeRecord(record);

            if(   (m_poFilterGeom == nullptr
                || FilterGeometry( poFeature->GetGeometryRef() ) )
            && (m_poAttrQuery == nullptr
                || m_poAttrQuery->Evaluate( poFeature )) )
            {
                 return poFeature;
            }

            delete poFeature;
        }
        else
        {
            BNA_FreeRecord(record);
        }
    }
}

/************************************************************************/
/*                      WriteFeatureAttributes()                        */
/************************************************************************/

void OGRBNALayer::WriteFeatureAttributes( VSILFILE* fp, OGRFeature *poFeature )
{
    int nbOutID = poDS->GetNbOutId();
    if (nbOutID < 0)
        nbOutID = poFeatureDefn->GetFieldCount();
    for(int i = 0; i < nbOutID; i++ )
    {
        if (i < poFeatureDefn->GetFieldCount())
        {
            OGRFieldDefn *poFieldDefn = poFeatureDefn->GetFieldDefn( i );
            if( poFeature->IsFieldSetAndNotNull( i ) )
            {
                if (poFieldDefn->GetType() == OFTReal)
                {
                    char szBuffer[64];
                    OGRFormatDouble(szBuffer, sizeof(szBuffer),
                                    poFeature->GetFieldAsDouble(i), '.');
                    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "\"%s\",", szBuffer));
                }
                else
                {
                    const char *pszRaw = poFeature->GetFieldAsString( i );
                    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "\"%s\",", pszRaw));
                }
            }
            else
            {
                CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "\"\","));
            }
        }
        else
        {
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "\"\","));
        }
    }
}

/************************************************************************/
/*                             WriteCoord()                             */
/************************************************************************/

void OGRBNALayer::WriteCoord(VSILFILE* fp, double dfX, double dfY)
{
    char szBuffer[64];
    OGRFormatDouble(szBuffer, sizeof(szBuffer), dfX, '.', poDS->GetCoordinatePrecision());
    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", szBuffer));
    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", poDS->GetCoordinateSeparator()));
    OGRFormatDouble(szBuffer, sizeof(szBuffer), dfY, '.', poDS->GetCoordinatePrecision());
    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", szBuffer));
}

/************************************************************************/
/*                           ICreateFeature()                            */
/************************************************************************/

OGRErr OGRBNALayer::ICreateFeature( OGRFeature *poFeature )

{
    OGRGeometry     *poGeom = poFeature->GetGeometryRef();
    char eol[3];
    const char* partialEol = (poDS->GetMultiLine()) ? eol : poDS->GetCoordinateSeparator();

    if (poGeom == nullptr || poGeom->IsEmpty() )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "OGR BNA driver cannot write features with empty geometries.");
        return OGRERR_FAILURE;
    }

    if( poDS->GetUseCRLF() )
    {
        eol[0] = 13;
        eol[1] = 10;
        eol[2] = 0;
    }
    else
    {
        eol[0] = 10;
        eol[1] = 0;
    }

    VSILFILE* fp = poDS->GetOutputFP();
    if ( ! bWriter || fp == nullptr )
    {
        return OGRERR_FAILURE;
    }

    if( poFeature->GetFID() == OGRNullFID )
        poFeature->SetFID( nFeatures++ );

    int nbPairPerLine = poDS->GetNbPairPerLine();

    switch( poGeom->getGeometryType() )
    {
        case wkbPoint:
        case wkbPoint25D:
        {
            OGRPoint* point = poGeom->toPoint();
            WriteFeatureAttributes(fp, poFeature);
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "1"));
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", partialEol));
            WriteCoord(fp, point->getX(), point->getY());
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", eol));
            break;
        }

        case wkbPolygon:
        case wkbPolygon25D:
        {
            OGRPolygon* polygon = poGeom->toPolygon();
            OGRLinearRing* ring = polygon->getExteriorRing();
            if (ring == nullptr)
            {
                return OGRERR_FAILURE;
            }

            double firstX = ring->getX(0);
            double firstY = ring->getY(0);
            int nBNAPoints = ring->getNumPoints();
            bool is_ellipse = false;

            /* This code tries to detect an ellipse in a polygon geometry */
            /* This will only work presumably on ellipses already read from a BNA file */
            /* Mostly a BNA to BNA feature... */
            if (poDS->GetEllipsesAsEllipses() &&
                polygon->getNumInteriorRings() == 0 &&
                nBNAPoints == 361)
            {
                double oppositeX = ring->getX(180);
                double oppositeY = ring->getY(180);
                double quarterX = ring->getX(90);
                double quarterY = ring->getY(90);
                double antiquarterX = ring->getX(270);
                double antiquarterY = ring->getY(270);
                double center1X = 0.5*(firstX + oppositeX);
                double center1Y = 0.5*(firstY + oppositeY);
                double center2X = 0.5*(quarterX + antiquarterX);
                double center2Y = 0.5*(quarterY + antiquarterY);
                if (fabs(center1X - center2X) < 1e-5 && fabs(center1Y - center2Y) < 1e-5 &&
                    fabs(oppositeY - firstY) < 1e-5 &&
                    fabs(quarterX - antiquarterX) < 1e-5)
                {
                    double major_radius = fabs(firstX - center1X);
                    double minor_radius = fabs(quarterY - center1Y);
                    is_ellipse = true;
                    for(int i=0;i<360;i++)
                    {
                        if (!(fabs(center1X + major_radius * cos(i * (M_PI / 180)) - ring->getX(i)) < 1e-5 &&
                              fabs(center1Y + minor_radius * sin(i * (M_PI / 180)) - ring->getY(i)) < 1e-5))
                        {
                            is_ellipse = false;
                            break;
                        }
                    }
                    if( is_ellipse )
                    {
                        WriteFeatureAttributes(fp, poFeature);
                        CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "2"));
                        CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", partialEol));
                        WriteCoord(fp, center1X, center1Y);
                        CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", partialEol));
                        WriteCoord(fp, major_radius, minor_radius);
                        CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", eol));
                    }
                }
            }

            if( !is_ellipse )
            {
                int nInteriorRings = polygon->getNumInteriorRings();
                for(int i=0;i<nInteriorRings;i++)
                {
                    nBNAPoints += polygon->getInteriorRing(i)->getNumPoints() + 1;
                }
                if (nBNAPoints <= 3)
                {
                    CPLError( CE_Failure, CPLE_AppDefined, "Invalid geometry" );
                    return OGRERR_FAILURE;
                }
                WriteFeatureAttributes(fp, poFeature);
                CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%d", nBNAPoints));
                int n = ring->getNumPoints();
                int nbPair = 0;
                for(int i=0;i<n;i++)
                {
                    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " "));
                    WriteCoord(fp, ring->getX(i), ring->getY(i));
                    nbPair++;
                }
                for(int i=0;i<nInteriorRings;i++)
                {
                    ring = polygon->getInteriorRing(i);
                    n = ring->getNumPoints();
                    for(int j=0;j<n;j++)
                    {
                        CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " "));
                        WriteCoord(fp, ring->getX(j), ring->getY(j));
                        nbPair++;
                    }
                    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " "));
                    WriteCoord(fp, firstX, firstY);
                    nbPair++;
                }
                CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", eol));
            }
            break;
        }

        case wkbMultiPolygon:
        case wkbMultiPolygon25D:
        {
            OGRMultiPolygon* multipolygon = poGeom->toMultiPolygon();
            int nBNAPoints = 0;
            double firstX = 0.0;
            double firstY = 0.0;
            for( auto&& polygon: *multipolygon )
            {
                OGRLinearRing* ring = polygon->getExteriorRing();
                if (ring == nullptr)
                    continue;

                if (nBNAPoints)
                    nBNAPoints ++;
                else
                {
                    firstX = ring->getX(0);
                    firstY = ring->getY(0);
                }
                nBNAPoints += ring->getNumPoints();
                int nInteriorRings = polygon->getNumInteriorRings();
                for(int j=0;j<nInteriorRings;j++)
                {
                    nBNAPoints += polygon->getInteriorRing(j)->getNumPoints() + 1;
                }
            }
            if (nBNAPoints <= 3)
            {
                CPLError( CE_Failure, CPLE_AppDefined, "Invalid geometry" );
                return OGRERR_FAILURE;
            }
            WriteFeatureAttributes(fp, poFeature);
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%d", nBNAPoints));
            int nbPair = 0;
            bool bFirst = true;
            for( auto&& polygon: *multipolygon )
            {
                OGRLinearRing* ring = polygon->getExteriorRing();
                if (ring == nullptr)
                    continue;

                int n = ring->getNumPoints();
                int nInteriorRings = polygon->getNumInteriorRings();
                for(int j=0;j<n;j++)
                {
                    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " "));
                    WriteCoord(fp, ring->getX(j), ring->getY(j));
                    nbPair++;
                }
                if (!bFirst)
                {
                    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " "));
                    WriteCoord(fp, firstX, firstY);
                    nbPair++;
                }
                bFirst = false;
                for(int j=0;j<nInteriorRings;j++)
                {
                    ring = polygon->getInteriorRing(j);
                    n = ring->getNumPoints();
                    for(int k=0;k<n;k++)
                    {
                        CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " "));
                        WriteCoord(fp, ring->getX(k), ring->getY(k));
                        nbPair++;
                    }
                    CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " "));
                    WriteCoord(fp, firstX, firstY);
                    nbPair++;
                }
            }
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", eol));
            break;
        }

        case wkbLineString:
        case wkbLineString25D:
        {
            OGRLineString* line = poGeom->toLineString();
            int n = line->getNumPoints();
            if (n < 2)
            {
                CPLError( CE_Failure, CPLE_AppDefined, "Invalid geometry" );
                return OGRERR_FAILURE;
            }
            WriteFeatureAttributes(fp, poFeature);
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "-%d", n));
            int nbPair = 0;
            for(int i=0;i<n;i++)
            {
                CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", partialEol));
                WriteCoord(fp, line->getX(i), line->getY(i));
                nbPair++;
            }
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", eol));
            break;
        }

        default:
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Unsupported geometry type : %s.",
                      poGeom->getGeometryName() );

            return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
        }
    }

    return OGRERR_NONE;
}

/************************************************************************/
/*                            CreateField()                             */
/************************************************************************/

OGRErr OGRBNALayer::CreateField( OGRFieldDefn *poField, int /* bApproxOK */ )
{
    if( !bWriter || nFeatures != 0)
        return OGRERR_FAILURE;

    poFeatureDefn->AddFieldDefn( poField );

    return OGRERR_NONE;
}

/************************************************************************/
/*                           BuildFeatureFromBNARecord()                */
/************************************************************************/
OGRFeature *OGRBNALayer::BuildFeatureFromBNARecord (BNARecord* record, long fid)
{
    OGRFeature  *poFeature = new OGRFeature( poFeatureDefn );
    for( int i = 0; i < nIDs; i++ )
    {
        poFeature->SetField( i, record->ids[i] ? record->ids[i] : "");
    }
    poFeature->SetFID( fid );
    if (bnaFeatureType == BNA_POINT)
    {
        poFeature->SetGeometryDirectly(
            new OGRPoint( record->tabCoords[0][0], record->tabCoords[0][1] ) );
    }
    else if (bnaFeatureType == BNA_POLYLINE)
    {
        OGRLineString* lineString = new OGRLineString ();
        lineString->setCoordinateDimension(2);
        lineString->setNumPoints(record->nCoords);
        for( int i = 0; i < record->nCoords; i++)
        {
            lineString->setPoint( i, record->tabCoords[i][0],
                                  record->tabCoords[i][1] );
        }
        poFeature->SetGeometryDirectly(lineString);
    }
    else if (bnaFeatureType == BNA_POLYGON)
    {
        double firstX = record->tabCoords[0][0];
        double firstY = record->tabCoords[0][1];
        int isFirstPolygon = 1;
        double secondaryFirstX = 0.0;
        double secondaryFirstY = 0.0;
        int iStartIndex = 0;

        /* record->nCoords is really a safe upper bound */
        int nbPolygons = 0;
        OGRPolygon** tabPolygons =
            static_cast<OGRPolygon**>(
                CPLMalloc( record->nCoords * sizeof(OGRPolygon*) ) );

        int i = 1;
        for( ; i < record->nCoords; i++ )
        {
            if (isFirstPolygon == 1 &&
                record->tabCoords[i][0] == firstX &&
                record->tabCoords[i][1] == firstY)
            {
                OGRPolygon* polygon = new OGRPolygon ();
                OGRLinearRing* ring = new OGRLinearRing ();
                ring->setCoordinateDimension(2);
                ring->setNumPoints(i - iStartIndex + 1, false);
                for( int j = iStartIndex; j <= i ; ++j )
                {
                    ring->setPoint( j - iStartIndex,
                                    record->tabCoords[j][0],
                                    record->tabCoords[j][1] );
                }
                polygon->addRingDirectly(ring);
                tabPolygons[nbPolygons] = polygon;
                nbPolygons++;

                if (i == record->nCoords - 1)
                {
                    break;
                }

                isFirstPolygon = 0;

                i++;
                secondaryFirstX = record->tabCoords[i][0];
                secondaryFirstY = record->tabCoords[i][1];
                iStartIndex = i;
            }
            else if (isFirstPolygon == 0 &&
                    record->tabCoords[i][0] == secondaryFirstX &&
                    record->tabCoords[i][1] == secondaryFirstY)
            {

                OGRPolygon* polygon = new OGRPolygon ();
                OGRLinearRing* ring = new OGRLinearRing ();
                ring->setCoordinateDimension(2);
                ring->setNumPoints(i - iStartIndex + 1, false);
                for( int j = iStartIndex; j <= i ; ++j )
                {
                    ring->setPoint( j - iStartIndex,
                                    record->tabCoords[j][0],
                                    record->tabCoords[j][1] );
                }
                polygon->addRingDirectly(ring);

                for( int j = 0; j < nbPolygons; j++ )
                {
                    if( polygon->Equals(tabPolygons[j]) )
                    {
                        delete polygon;
                        polygon = nullptr;
                        break;
                    }
                }
                if( polygon )
                {
                    tabPolygons[nbPolygons] = polygon;
                    nbPolygons++;
                }

                if (i < record->nCoords - 1)
                {
                    // After the closing of a subpolygon, the first
                    // coordinates of the first polygon should be
                    // recalled... in theory
                    if (record->tabCoords[i+1][0] == firstX
                        && record->tabCoords[i+1][1] == firstY)
                    {
                        if (i + 1 == record->nCoords - 1)
                            break;
                        i ++;
                    }
                    else
                    {
#if 0
                        CPLError(CE_Warning, CPLE_AppDefined,
                                 "Geometry of polygon of fid %d starting at "
                                 "line %d is not strictly conformant. "
                                 "Trying to go on...\n",
                                 fid,
                                 offsetAndLineFeaturesTable[fid].line + 1);
#endif
                    }

                    i ++;
                    secondaryFirstX = record->tabCoords[i][0];
                    secondaryFirstY = record->tabCoords[i][1];
                    iStartIndex = i;
                }
                else
                {
#if 0
                    CPLError( CE_Warning, CPLE_AppDefined,
                              "Geometry of polygon of fid %d starting at "
                              "line %d is not strictly conformant. "
                              "Trying to go on...\n",
                              fid,
                              offsetAndLineFeaturesTable[fid].line + 1);
#endif
                }
            }
        }
        if (i == record->nCoords)
        {
            /* Let's be a bit tolerant about non-closing polygons. */
            if (isFirstPolygon)
            {
                CPLAssert( iStartIndex == 0 );
                OGRLinearRing* ring = new OGRLinearRing ();
                ring->setCoordinateDimension(2);
                ring->setNumPoints(record->nCoords + 1, false);
                for( int j = 0; j < record->nCoords ; ++j )
                {
                    ring->setPoint( j,
                                    record->tabCoords[j][0],
                                    record->tabCoords[j][1] );
                }

                ring->setPoint( record->nCoords,
                                record->tabCoords[0][0],
                                record->tabCoords[0][1] );

                OGRPolygon* polygon = new OGRPolygon ();
                polygon->addRingDirectly(ring);

                tabPolygons[nbPolygons] = polygon;
                nbPolygons++;
            }
        }

        if (nbPolygons == 1)
        {
            // Special optimization here : we directly put the polygon into
            // the multipolygon.  This should save quite a few useless copies
            OGRMultiPolygon* multipolygon = new OGRMultiPolygon();
            multipolygon->addGeometryDirectly(tabPolygons[0]);
            poFeature->SetGeometryDirectly(multipolygon);
        }
        else
        {
            int isValidGeometry = FALSE;
            poFeature->SetGeometryDirectly(
                OGRGeometryFactory::organizePolygons(
                    reinterpret_cast<OGRGeometry**>( tabPolygons ),
                    nbPolygons, &isValidGeometry, nullptr ) );

            if (!isValidGeometry)
            {
                CPLError( CE_Warning, CPLE_AppDefined,
                          "Geometry of polygon of fid %ld starting at line %d "
                          "cannot be translated to Simple Geometry. "
                          "All polygons will be contained in a multipolygon.\n",
                          fid,
                          offsetAndLineFeaturesTable[fid].line + 1);
            }
        }

        CPLFree(tabPolygons);
    }
    else
    {
        // Circle or ellipses are not part of the OGR Simple Geometry, so we
        // discretize them into polygons by 1 degree step.
        OGRPolygon* polygon = new OGRPolygon ();
        OGRLinearRing* ring = new OGRLinearRing ();
        ring->setCoordinateDimension(2);
        double center_x = record->tabCoords[0][0];
        double center_y = record->tabCoords[0][1];
        double major_radius = record->tabCoords[1][0];
        double minor_radius = record->tabCoords[1][1];
        if (minor_radius == 0)
            minor_radius = major_radius;
        for( int i = 0; i < 360; i++ )
        {
            ring->addPoint(center_x + major_radius * cos(i * (M_PI / 180)),
                           center_y + minor_radius * sin(i * (M_PI / 180)) );
        }
        ring->addPoint(center_x + major_radius, center_y);
        polygon->addRingDirectly  (  ring );
        poFeature->SetGeometryDirectly(polygon);

        poFeature->SetField( nIDs, major_radius);
        poFeature->SetField( nIDs+1, minor_radius);
    }

    return poFeature;
}

/************************************************************************/
/*                           FastParseUntil()                           */
/************************************************************************/
void OGRBNALayer::FastParseUntil ( int interestFID)
{
    if( partialIndexTable )
    {
        ResetReading();

        BNARecord *record = nullptr;

        if (nFeatures > 0)
        {
            if( VSIFSeekL( fpBNA,
                           offsetAndLineFeaturesTable[nFeatures-1].offset,
                           SEEK_SET ) < 0 )
                return;
            curLine = offsetAndLineFeaturesTable[nFeatures-1].line;

            /* Just skip the last read one */
            int ok = FALSE;
            record = BNA_GetNextRecord(fpBNA, &ok, &curLine, TRUE,
                                       BNA_READ_NONE);
            BNA_FreeRecord(record);
        }

        while(1)
        {
            int ok = FALSE;
            int offset = static_cast<int>( VSIFTellL(fpBNA) );
            int line = curLine;
            record = BNA_GetNextRecord(fpBNA, &ok, &curLine, TRUE,
                                       BNA_READ_NONE);
            if (ok == FALSE)
            {
                failed = true;
                return;
            }
            if (record == nullptr)
            {
                /* end of file */
                eof = true;

                /* and we have finally build the whole index table */
                partialIndexTable = false;
                return;
            }

            if (record->featureType == bnaFeatureType)
            {
                nFeatures++;
                offsetAndLineFeaturesTable =
                    static_cast<OffsetAndLine *>(
                        CPLRealloc( offsetAndLineFeaturesTable,
                                    nFeatures * sizeof(OffsetAndLine) ) );
                offsetAndLineFeaturesTable[nFeatures-1].offset = offset;
                offsetAndLineFeaturesTable[nFeatures-1].line = line;

                BNA_FreeRecord(record);

                if (nFeatures - 1 == interestFID)
                  return;
            }
            else
            {
                BNA_FreeRecord(record);
            }
        }
    }
}

/************************************************************************/
/*                           GetFeature()                               */
/************************************************************************/

OGRFeature *  OGRBNALayer::GetFeature( GIntBig nFID )
{
    if (nFID < 0 || !CPL_INT64_FITS_ON_INT32(nFID))
        return nullptr;

    FastParseUntil( static_cast<int>( nFID ) );

    if (nFID >= nFeatures)
        return nullptr;

    if( VSIFSeekL( fpBNA, offsetAndLineFeaturesTable[nFID].offset, SEEK_SET ) < 0 )
        return nullptr;

    curLine = offsetAndLineFeaturesTable[nFID].line;
    int ok = FALSE;
    BNARecord* record
        = BNA_GetNextRecord(fpBNA, &ok, &curLine, TRUE, bnaFeatureType);

    OGRFeature *poFeature = BuildFeatureFromBNARecord(record, (int)nFID);

    BNA_FreeRecord(record);

    return poFeature;
}

/************************************************************************/
/*                           TestCapability()                           */
/************************************************************************/

int OGRBNALayer::TestCapability( const char * pszCap )

{
    if( EQUAL(pszCap,OLCSequentialWrite) )
        return bWriter;
    else if( EQUAL(pszCap,OLCCreateField) )
        return bWriter && nFeatures == 0;
    else
        return FALSE;
}