/******************************************************************************
 *
 * Project:  Interlis 1 Translator
 * Purpose:  Implements OGRILI1Layer class.
 * Author:   Pirmin Kalberer, Sourcepole AG
 *
 ******************************************************************************
 * Copyright (c) 2004, Pirmin Kalberer, Sourcepole AG
 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_conv.h"
#include "cpl_string.h"
#include "ogr_geos.h"
#include "ogr_ili1.h"

#include <map>
#include <vector>

/************************************************************************/
/*                           OGRILI1Layer()                              */
/************************************************************************/

OGRILI1Layer::OGRILI1Layer(OGRFeatureDefn *poFeatureDefnIn,
                           const GeomFieldInfos &oGeomFieldInfosIn,
                           OGRILI1DataSource *poDSIn)
    : poFeatureDefn(poFeatureDefnIn), oGeomFieldInfos(oGeomFieldInfosIn),
      nFeatures(0), papoFeatures(nullptr), nFeatureIdx(0), bGeomsJoined(FALSE),
      poDS(poDSIn)
{
    SetDescription(poFeatureDefn->GetName());
    poFeatureDefn->Reference();
}

/************************************************************************/
/*                           ~OGRILI1Layer()                           */
/************************************************************************/

OGRILI1Layer::~OGRILI1Layer()
{
    for (int i = 0; i < nFeatures; i++)
    {
        delete papoFeatures[i];
    }
    CPLFree(papoFeatures);

    if (poFeatureDefn)
        poFeatureDefn->Release();
}

OGRErr OGRILI1Layer::AddFeature(OGRFeature *poFeature)
{
    nFeatures++;

    papoFeatures = static_cast<OGRFeature **>(
        CPLRealloc(papoFeatures, sizeof(void *) * nFeatures));

    papoFeatures[nFeatures - 1] = poFeature;

    return OGRERR_NONE;
}

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

void OGRILI1Layer::ResetReading()
{
    nFeatureIdx = 0;
}

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

OGRFeature *OGRILI1Layer::GetNextFeature()
{
    if (!bGeomsJoined)
        JoinGeomLayers();

    while (nFeatureIdx < nFeatures)
    {
        OGRFeature *poFeature = GetNextFeatureRef();
        if (poFeature)
            return poFeature->Clone();
    }
    return nullptr;
}

OGRFeature *OGRILI1Layer::GetNextFeatureRef()
{
    OGRFeature *poFeature = nullptr;
    if (nFeatureIdx < nFeatures)
    {
        poFeature = papoFeatures[nFeatureIdx++];
        // apply filters
        if ((m_poFilterGeom == nullptr ||
             FilterGeometry(poFeature->GetGeometryRef())) &&
            (m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
            return poFeature;
    }
    return nullptr;
}

/************************************************************************/
/*                             GetFeatureRef()                          */
/************************************************************************/

OGRFeature *OGRILI1Layer::GetFeatureRef(GIntBig nFID)

{
    ResetReading();

    OGRFeature *poFeature = nullptr;
    while ((poFeature = GetNextFeatureRef()) != nullptr)
    {
        if (poFeature->GetFID() == nFID)
            return poFeature;
    }

    return nullptr;
}

OGRFeature *OGRILI1Layer::GetFeatureRef(const char *fid)

{
    ResetReading();

    OGRFeature *poFeature = nullptr;
    while ((poFeature = GetNextFeatureRef()) != nullptr)
    {
        if (!strcmp(poFeature->GetFieldAsString(0), fid))
            return poFeature;
    }

    return nullptr;
}

/************************************************************************/
/*                          GetFeatureCount()                           */
/************************************************************************/

GIntBig OGRILI1Layer::GetFeatureCount(int bForce)
{
    if (m_poFilterGeom == nullptr && m_poAttrQuery == nullptr
        /* && poAreaLineLayer == NULL*/)
    {
        return nFeatures;
    }
    else
    {
        return OGRLayer::GetFeatureCount(bForce);
    }
}

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

int OGRILI1Layer::TestCapability(const char *pszCap) const
{
    if (EQUAL(pszCap, OLCCurveGeometries))
        return TRUE;
    if (EQUAL(pszCap, OLCZGeometries))
        return TRUE;

    return FALSE;
}

/************************************************************************/
/*                         Internal routines                            */
/************************************************************************/

void OGRILI1Layer::JoinGeomLayers()
{
    bGeomsJoined = true;

    CPLConfigOptionSetter oSetter("OGR_ARC_STEPSIZE", "0.96",
                                  /* bSetOnlyIfUndefined = */ true);

    for (GeomFieldInfos::const_iterator it = oGeomFieldInfos.begin();
         it != oGeomFieldInfos.end(); ++it)
    {
        OGRFeatureDefn *geomFeatureDefn = it->second.GetGeomTableDefnRef();
        if (geomFeatureDefn)
        {
            CPLDebug("OGR_ILI", "Join geometry table %s of field '%s'",
                     geomFeatureDefn->GetName(), it->first.c_str());
            OGRILI1Layer *poGeomLayer =
                poDS->GetLayerByName(geomFeatureDefn->GetName());
            const int nGeomFieldIndex =
                GetLayerDefn()->GetGeomFieldIndex(it->first.c_str());
            if (it->second.iliGeomType == "Surface")
            {
                JoinSurfaceLayer(poGeomLayer, nGeomFieldIndex);
            }
            else if (it->second.iliGeomType == "Area")
            {
                CPLString pointField = it->first + "__Point";
                const int nPointFieldIndex =
                    GetLayerDefn()->GetGeomFieldIndex(pointField.c_str());
                PolygonizeAreaLayer(poGeomLayer, nGeomFieldIndex,
                                    nPointFieldIndex);
            }
        }
    }
}

void OGRILI1Layer::JoinSurfaceLayer(OGRILI1Layer *poSurfaceLineLayer,
                                    int nSurfaceFieldIndex)
{
    CPLDebug("OGR_ILI", "Joining surface layer %s with geometries",
             GetLayerDefn()->GetName());
    OGRwkbGeometryType geomType =
        GetLayerDefn()->GetGeomFieldDefn(nSurfaceFieldIndex)->GetType();

    std::map<OGRFeature *, std::vector<OGRCurve *>> oMapFeatureToGeomSet;

    poSurfaceLineLayer->ResetReading();

    // First map: for each target curvepolygon, find all belonging curves
    while (OGRFeature *linefeature = poSurfaceLineLayer->GetNextFeatureRef())
    {
        // OBJE entries with same _RefTID are polygon rings of same feature
        OGRFeature *feature = nullptr;
        if (poFeatureDefn->GetFieldDefn(0)->GetType() == OFTString)
        {
            feature = GetFeatureRef(linefeature->GetFieldAsString(1));
        }
        else
        {
            GIntBig reftid = linefeature->GetFieldAsInteger64(1);
            feature = GetFeatureRef(reftid);
        }
        if (feature)
        {
            OGRGeometry *poGeom = linefeature->GetGeomFieldRef(0);
            OGRMultiCurve *curves = dynamic_cast<OGRMultiCurve *>(poGeom);
            if (curves)
            {
                for (auto &&curve : curves)
                {
                    if (!curve->IsEmpty())
                        oMapFeatureToGeomSet[feature].push_back(curve);
                }
            }
        }
        else
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Couldn't join feature FID " CPL_FRMT_GIB,
                     linefeature->GetFieldAsInteger64(1));
        }
    }

    // Now for each target polygon, assemble the curves together.
    std::map<OGRFeature *, std::vector<OGRCurve *>>::const_iterator oIter =
        oMapFeatureToGeomSet.begin();
    for (; oIter != oMapFeatureToGeomSet.end(); ++oIter)
    {
        OGRFeature *feature = oIter->first;
        std::vector<OGRCurve *> oCurves = oIter->second;

        std::vector<OGRCurve *> oSetDestCurves;
        double dfLargestArea = 0.0;
        OGRCurve *poLargestCurve = nullptr;
        while (true)
        {
            std::vector<OGRCurve *>::iterator oIterCurves = oCurves.begin();
            if (oIterCurves == oCurves.end())
                break;

            OGRPoint endPointCC;
            OGRCompoundCurve *poCC = new OGRCompoundCurve();

            bool bFirst = true;
            while (true)
            {
                bool bNewCurveAdded = false;
                const double dfEps = 1e-14;
                for (oIterCurves = oCurves.begin();
                     oIterCurves != oCurves.end(); ++oIterCurves)
                {
                    OGRCurve *curve = *oIterCurves;
                    OGRPoint startPoint;
                    OGRPoint endPoint;
                    curve->StartPoint(&startPoint);
                    curve->EndPoint(&endPoint);
                    if (bFirst ||
                        (fabs(startPoint.getX() - endPointCC.getX()) < dfEps &&
                         fabs(startPoint.getY() - endPointCC.getY()) < dfEps))
                    {
                        bFirst = false;

                        curve->EndPoint(&endPointCC);

                        const OGRwkbGeometryType eCurveType =
                            wkbFlatten(curve->getGeometryType());
                        if (eCurveType == wkbCompoundCurve)
                        {
                            OGRCompoundCurve *poCCSub =
                                curve->toCompoundCurve();
                            for (auto &&subCurve : poCCSub)
                            {
                                poCC->addCurve(subCurve);
                            }
                        }
                        else
                        {
                            poCC->addCurve(curve);
                        }
                        oCurves.erase(oIterCurves);
                        bNewCurveAdded = true;
                        break;
                    }
                    else if (fabs(endPoint.getX() - endPointCC.getX()) <
                                 dfEps &&
                             fabs(endPoint.getY() - endPointCC.getY()) < dfEps)
                    {
                        curve->StartPoint(&endPointCC);

                        const OGRwkbGeometryType eCurveType =
                            wkbFlatten(curve->getGeometryType());
                        if (eCurveType == wkbLineString ||
                            eCurveType == wkbCircularString)
                        {
                            OGRSimpleCurve *poSC =
                                curve->toSimpleCurve()->clone();
                            poSC->reversePoints();
                            poCC->addCurveDirectly(poSC);
                        }
                        else if (eCurveType == wkbCompoundCurve)
                        {
                            // Reverse the order of the elements of the
                            // compound curve
                            OGRCompoundCurve *poCCSub =
                                curve->toCompoundCurve();
                            for (int i = poCCSub->getNumCurves() - 1; i >= 0;
                                 --i)
                            {
                                OGRSimpleCurve *poSC = poCCSub->getCurve(i)
                                                           ->toSimpleCurve()
                                                           ->clone();
                                poSC->reversePoints();
                                poCC->addCurveDirectly(poSC);
                            }
                        }

                        oCurves.erase(oIterCurves);
                        bNewCurveAdded = true;
                        break;
                    }
                }
                if (!bNewCurveAdded || oCurves.empty() || poCC->get_IsClosed())
                    break;
            }

            if (!poCC->get_IsClosed())
            {
                char *pszJSon = poCC->exportToJson();
                CPLError(CE_Warning, CPLE_AppDefined,
                         "A ring %s for feature " CPL_FRMT_GIB " in layer %s "
                         "was not closed. Dropping it",
                         pszJSon, feature->GetFID(), GetName());
                delete poCC;
                CPLFree(pszJSon);
            }
            else
            {
                double dfArea = poCC->get_Area();
                if (dfArea >= dfLargestArea)
                {
                    dfLargestArea = dfArea;
                    poLargestCurve = poCC;
                }
                oSetDestCurves.push_back(poCC);
            }
        }

        // Now build the final polygon by first inserting the largest ring.
        OGRCurvePolygon *poPoly =
            (geomType == wkbPolygon) ? new OGRPolygon() : new OGRCurvePolygon();
        if (poLargestCurve)
        {
            std::vector<OGRCurve *>::iterator oIterCurves =
                oSetDestCurves.begin();
            for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
            {
                OGRCurve *poCurve = *oIterCurves;
                if (poCurve == poLargestCurve)
                {
                    oSetDestCurves.erase(oIterCurves);
                    break;
                }
            }

            if (geomType == wkbPolygon)
            {
                poLargestCurve = OGRCurve::CastToLinearRing(poLargestCurve);
            }
            OGRErr error = poPoly->addRingDirectly(poLargestCurve);
            if (error != OGRERR_NONE)
            {
                char *pszJSon = poLargestCurve->exportToJson();
                CPLError(CE_Warning, CPLE_AppDefined,
                         "Cannot add ring %s to feature " CPL_FRMT_GIB
                         " in layer %s",
                         pszJSon, feature->GetFID(), GetName());
                CPLFree(pszJSon);
            }

            oIterCurves = oSetDestCurves.begin();
            for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
            {
                OGRCurve *poCurve = *oIterCurves;
                if (geomType == wkbPolygon)
                {
                    poCurve = OGRCurve::CastToLinearRing(poCurve);
                }
                error = poPoly->addRingDirectly(poCurve);
                if (error != OGRERR_NONE)
                {
                    char *pszJSon = poCurve->exportToJson();
                    CPLError(CE_Warning, CPLE_AppDefined,
                             "Cannot add ring %s to feature " CPL_FRMT_GIB
                             " in layer %s",
                             pszJSon, feature->GetFID(), GetName());
                    CPLFree(pszJSon);
                }
            }
        }

        feature->SetGeomFieldDirectly(nSurfaceFieldIndex, poPoly);
    }

    ResetReading();
}

OGRMultiPolygon *OGRILI1Layer::Polygonize(OGRGeometryCollection *poLines,
                                          bool
#if defined(HAVE_GEOS)
                                              fix_crossing_lines
#endif
)
{
    if (poLines->getNumGeometries() == 0)
    {
        return new OGRMultiPolygon();
    }

#if defined(HAVE_GEOS)
    GEOSGeom *ahInGeoms = nullptr;
    OGRGeometryCollection *poNoncrossingLines = poLines;
    GEOSGeom hResultGeom = nullptr;
    OGRGeometry *poMP = nullptr;

    if (fix_crossing_lines && poLines->getNumGeometries() > 0)
    {
        CPLDebug("OGR_ILI", "Fixing crossing lines");
        // A union of the geometry collection with one line fixes
        // invalid geometries.
        OGRGeometry *poUnion = poLines->Union(poLines->getGeometryRef(0));
        if (poUnion != nullptr)
        {
            if (wkbFlatten(poUnion->getGeometryType()) ==
                    wkbGeometryCollection ||
                wkbFlatten(poUnion->getGeometryType()) == wkbMultiLineString)
            {
                poNoncrossingLines =
                    dynamic_cast<OGRGeometryCollection *>(poUnion);
                CPLDebug("OGR_ILI", "Fixed lines: %d",
                         poNoncrossingLines->getNumGeometries() -
                             poLines->getNumGeometries());
            }
            else
            {
                delete poUnion;
            }
        }
    }

    GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();

    ahInGeoms = static_cast<GEOSGeom *>(
        CPLCalloc(sizeof(void *), poNoncrossingLines->getNumGeometries()));
    for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
        ahInGeoms[i] =
            poNoncrossingLines->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);

    hResultGeom = GEOSPolygonize_r(hGEOSCtxt, ahInGeoms,
                                   poNoncrossingLines->getNumGeometries());

    for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
        GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
    CPLFree(ahInGeoms);
    if (poNoncrossingLines != poLines)
        delete poNoncrossingLines;

    if (hResultGeom == nullptr)
    {
        OGRGeometry::freeGEOSContext(hGEOSCtxt);
        return new OGRMultiPolygon();
    }

    poMP = OGRGeometryFactory::createFromGEOS(hGEOSCtxt, hResultGeom);

    GEOSGeom_destroy_r(hGEOSCtxt, hResultGeom);
    OGRGeometry::freeGEOSContext(hGEOSCtxt);

    poMP = OGRGeometryFactory::forceToMultiPolygon(poMP);
    if (poMP && wkbFlatten(poMP->getGeometryType()) == wkbMultiPolygon)
        return dynamic_cast<OGRMultiPolygon *>(poMP);

    delete poMP;
    return new OGRMultiPolygon();

#else
    return new OGRMultiPolygon();
#endif
}

void OGRILI1Layer::PolygonizeAreaLayer(OGRILI1Layer *poAreaLineLayer,
                                       int
#if defined(HAVE_GEOS)
                                           nAreaFieldIndex
#endif
                                       ,
                                       int
#if defined(HAVE_GEOS)
                                           nPointFieldIndex
#endif
)
{
    // add all lines from poAreaLineLayer to collection
    OGRGeometryCollection *gc = new OGRGeometryCollection();
    poAreaLineLayer->ResetReading();
    while (OGRFeature *feature = poAreaLineLayer->GetNextFeatureRef())
        gc->addGeometry(feature->GetGeometryRef());

    // polygonize lines
    CPLDebug("OGR_ILI", "Polygonizing layer %s with %d multilines",
             poAreaLineLayer->GetLayerDefn()->GetName(),
             gc->getNumGeometries());
    poAreaLineLayer = nullptr;
    OGRMultiPolygon *polys = Polygonize(gc, false);
    CPLDebug("OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
    if (polys->getNumGeometries() != GetFeatureCount())
    {
        CPLDebug("OGR_ILI", "Feature count of layer %s: " CPL_FRMT_GIB,
                 GetLayerDefn()->GetName(), GetFeatureCount());
        CPLDebug("OGR_ILI", "Polygonizing again with crossing line fix");
        delete polys;
        polys = Polygonize(gc, true);  // try again with crossing line fix
        CPLDebug("OGR_ILI", "Resulting polygons: %d",
                 polys->getNumGeometries());
    }
    delete gc;

    // associate polygon feature with data row according to centroid
#if defined(HAVE_GEOS)
    OGRPolygon emptyPoly;

    CPLDebug("OGR_ILI", "Associating layer %s with area polygons",
             GetLayerDefn()->GetName());
    GEOSGeom *ahInGeoms = static_cast<GEOSGeom *>(
        CPLCalloc(sizeof(void *), polys->getNumGeometries()));
    GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
    for (int i = 0; i < polys->getNumGeometries(); i++)
    {
        ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
        if (!GEOSisValid_r(hGEOSCtxt, ahInGeoms[i]))
            ahInGeoms[i] = nullptr;
    }
    for (int nFidx = 0; nFidx < nFeatures; nFidx++)
    {
        OGRFeature *feature = papoFeatures[nFidx];
        OGRGeometry *geomRef = feature->GetGeomFieldRef(nPointFieldIndex);
        if (!geomRef)
        {
            continue;
        }
        GEOSGeom point =
            reinterpret_cast<GEOSGeom>(geomRef->exportToGEOS(hGEOSCtxt));

        int i = 0;
        for (; i < polys->getNumGeometries(); i++)
        {
            if (ahInGeoms[i] && GEOSWithin_r(hGEOSCtxt, point, ahInGeoms[i]))
            {
                feature->SetGeomField(nAreaFieldIndex,
                                      polys->getGeometryRef(i));
                break;
            }
        }
        if (i == polys->getNumGeometries())
        {
            CPLDebug("OGR_ILI", "Association between area and point failed.");
            feature->SetGeometry(&emptyPoly);
        }
        GEOSGeom_destroy_r(hGEOSCtxt, point);
    }
    for (int i = 0; i < polys->getNumGeometries(); i++)
        GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
    CPLFree(ahInGeoms);
    OGRGeometry::freeGEOSContext(hGEOSCtxt);
#endif
    delete polys;
}

/************************************************************************/
/*                             GetDataset()                             */
/************************************************************************/

GDALDataset *OGRILI1Layer::GetDataset()
{
    return poDS;
}