/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2016 Daniel Baston
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Public Licence as published
 * by the Free Software Foundation.
 * See the COPYING file for more information.
 *
 **********************************************************************
 *
 * Last port: precision/MinimumClearance.java (f6187ee2 JTS-1.14)
 *
 **********************************************************************/

#include <geos/algorithm/Distance.h>
#include <geos/precision/MinimumClearance.h>
#include <geos/index/strtree/STRtree.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/operation/distance/FacetSequenceTreeBuilder.h>
#include <geos/geom/LineSegment.h>

using namespace geos::geom;
using namespace geos::operation::distance;
using namespace geos::index::strtree;

namespace geos {
namespace precision {

MinimumClearance::MinimumClearance(const Geometry* g) : inputGeom(g) {}

double
MinimumClearance::getDistance()
{
    compute();
    return minClearance;
}

std::unique_ptr<LineString>
MinimumClearance::getLine()
{
    compute();

    // return empty line string if no min pts were found
    if (minClearance == DoubleInfinity) {
        return inputGeom->getFactory()->createLineString();
    }

    return inputGeom->getFactory()->createLineString(minClearancePts->clone());
}

void
MinimumClearance::compute()
{
    class MinClearanceDistance {
    private:
        double minDist;
        std::vector<Coordinate> minPts;

        void
        updatePts(const Coordinate& p, const Coordinate& seg0, const Coordinate& seg1)
        {
            LineSegment seg(seg0, seg1);

            minPts[0] = p;
            seg.closestPoint(p, minPts[1]);
        }

    public:
        MinClearanceDistance() :
            minDist(DoubleInfinity),
            minPts(std::vector<Coordinate>(2))
        {}

        const std::vector<Coordinate>*
        getCoordinates()
        {
            return &minPts;
        }

        double operator()(const FacetSequence* fs1, const FacetSequence* fs2)
        {
            minDist = DoubleInfinity;
            return distance(fs1, fs2);
        }

        double
        distance(const FacetSequence* fs1, const FacetSequence* fs2)
        {
            // Compute MinClearance distance metric

            vertexDistance(fs1, fs2);
            if(fs1->size() == 1 && fs2->size() == 1) {
                return minDist;
            }
            if(minDist <= 0.0) {
                return minDist;
            }

            segmentDistance(fs1, fs2);
            if(minDist <= 0.0) {
                return minDist;
            }

            segmentDistance(fs2, fs1);
            return minDist;
        }

        double
        vertexDistance(const FacetSequence* fs1, const FacetSequence* fs2)
        {
            for(std::size_t i1 = 0; i1 < fs1->size(); i1++) {
                for(std::size_t i2 = 0; i2 < fs2->size(); i2++) {
                    const Coordinate* p1 = fs1->getCoordinate(i1);
                    const Coordinate* p2 = fs2->getCoordinate(i2);
                    if(!p1->equals2D(*p2)) {
                        double d = p1->distance(*p2);
                        if(d < minDist) {
                            minDist = d;
                            minPts[0] = *p1;
                            minPts[1] = *p2;
                            if(d == 0.0) {
                                return d;
                            }
                        }
                    }
                }
            }
            return minDist;
        }

        double
        segmentDistance(const FacetSequence* fs1, const FacetSequence* fs2)
        {
            for(std::size_t i1 = 0; i1 < fs1->size(); i1++) {
                for(std::size_t i2 = 1; i2 < fs2->size(); i2++) {
                    const Coordinate* p = fs1->getCoordinate(i1);

                    const Coordinate* seg0 = fs2->getCoordinate(i2 - 1);
                    const Coordinate* seg1 = fs2->getCoordinate(i2);

                    if(!(p->equals2D(*seg0) || p->equals2D(*seg1))) {
                        double d = geos::algorithm::Distance::pointToSegment(*p, *seg0, *seg1);
                        if(d < minDist) {
                            minDist = d;
                            updatePts(*p, *seg0, *seg1);
                            if(d == 0.0) {
                                return d;
                            }
                        }
                    }
                }
            }
            return minDist;
        }
    };

    // already computed
    if(minClearancePts != nullptr) {
        return;
    }

    // initialize to "No Distance Exists" state
    minClearancePts = detail::make_unique<CoordinateSequence>(2u, 2u);
    minClearance = DoubleInfinity;

    // handle empty geometries
    if(inputGeom->isEmpty()) {
        return;
    }

    auto tree = FacetSequenceTreeBuilder::build(inputGeom);
    MinClearanceDistance mcd;
    auto nearest = tree->nearestNeighbour(mcd);

    if (nearest.first == nullptr || nearest.second == nullptr) {
        throw util::GEOSException("Failed to find nearest items");
    }

    minClearance = mcd.distance(nearest.first, nearest.second);

    const std::vector<Coordinate>* minClearancePtsVec = mcd.getCoordinates();
    minClearancePts->setAt((*minClearancePtsVec)[0], 0);
    minClearancePts->setAt((*minClearancePtsVec)[1], 1);
}


}
}