/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2014 Mika Heiskanen <mika.heiskanen@fmi.fi>
 *
 * 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.
 *
 **********************************************************************/

#include <geos/algorithm/PointLocation.h>
#include <geos/algorithm/Orientation.h>
#include <geos/operation/intersection/RectangleIntersection.h>
#include <geos/operation/intersection/Rectangle.h>
#include <geos/operation/intersection/RectangleIntersectionBuilder.h>
#include <geos/operation/overlayng/ElevationModel.h>
#include <geos/operation/predicate/RectangleIntersects.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Polygon.h>
#include <geos/geom/MultiPolygon.h>
#include <geos/geom/Point.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/LineString.h>
#include <geos/geom/LinearRing.h>
#include <geos/geom/Location.h>
#include <geos/util/UnsupportedOperationException.h>
#include <geos/util.h>
#include <list>
#include <stdexcept>

using geos::operation::intersection::Rectangle;
using geos::operation::intersection::RectangleIntersectionBuilder;
using geos::operation::overlayng::ElevationModel;
using namespace geos::geom;
using namespace geos::algorithm;
namespace geos {
namespace operation { // geos::operation
namespace intersection { // geos::operation::intersection

/**
 * \brief Test if two coordinates are different
 */

inline
bool
different(double x1, double y1, double x2, double y2)
{
    return !(x1 == x2 && y1 == y2);
}

/**
 * \brief Calculate a line intersection point
 *
 * Note:
 *   - Calling this with x1,y1 and x2,y2 swapped cuts the other end of the line
 *   - Calling this with x and y swapped cuts in y-direction instead
 *   - Calling with 1<->2 and x<->y swapped works too
 */

inline
void
clip_one_edge(double& x1, double& y1, double& z1, double x2, double y2, double z2, double limit)
{
    if(x2 == limit) {
        y1 = y2;
        x1 = x2;
        z1 = z2;
    }

    if(x1 != x2) {
        double fraction = (limit - x1) / (x2 - x1);
        y1 += (y2 - y1) * fraction;
        z1 += (z2 - z1) * fraction;
        x1 = limit;
    }
}

/**
 * \brief Start point is outside, endpoint is definitely inside
 *
 * Note: Even though one might think using >= etc operators would produce
 *       the same result, that is not the case. We rely on the fact
 *       that nothing is clipped unless the point is truly outside the
 *       rectangle! Without this handling lines ending on the edges of
 *       the rectangle would be very difficult.
 */

void
clip_to_edges(double& x1, double& y1, double& z1,
              double x2, double y2, double z2,
              const Rectangle& rect)
{
    if(x1 < rect.xmin()) {
        clip_one_edge(x1, y1, z1, x2, y2, z2, rect.xmin());
    }
    else if(x1 > rect.xmax()) {
        clip_one_edge(x1, y1, z1, x2, y2, z2, rect.xmax());
    }

    if(y1 < rect.ymin()) {
        clip_one_edge(y1, x1, z1, y2, x2, z2, rect.ymin());
    }
    else if(y1 > rect.ymax()) {
        clip_one_edge(y1, x1, z1, y2, x2, z2, rect.ymax());
    }
}


/**
 * \brief Clip  geometry
 *
 * Here outGeom may also be a MultiPoint
 */

void
RectangleIntersection::clip_point(const geom::Point* g,
                                  RectangleIntersectionBuilder& parts,
                                  const Rectangle& rect)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }

    double x = g->getX();
    double y = g->getY();

    if(rect.position(x, y) == Rectangle::Inside) {
        parts.add(g->clone().release());
    }
}

bool
RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
        RectangleIntersectionBuilder& parts,
        const Rectangle& rect)
{
    auto n = gi->getNumPoints();

    if(gi == nullptr || n < 1) {
        return false;
    }

    // For shorthand code

    std::vector<Coordinate> cs;
    gi->getCoordinatesRO()->toVector(cs);
    //const geom::CoordinateSequence &cs = *(gi->getCoordinatesRO());

    // Keep a record of the point where a line segment entered
    // the rectangle. If the boolean is set, we must insert
    // the point to the beginning of the linestring which then
    // continues inside the rectangle.

    double x0 = 0;
    double y0 = 0;
    double z0 = 0;
    bool add_start = false;

    // Start iterating

    std::size_t i = 0;

    while(i < n) {
        // Establish initial position

        double x = cs[i].x;
        double y = cs[i].y;
        double z = cs[i].z;
        Rectangle::Position pos = rect.position(x, y);

        if(pos == Rectangle::Outside) {
            // Skip points as fast as possible until something has to be checked
            // in more detail.

            ++i;	// we already know it is outside

            if(x < rect.xmin())
                while(i < n && cs[i].x < rect.xmin()) {
                    ++i;
                }

            else if(x > rect.xmax())
                while(i < n && cs[i].x > rect.xmax()) {
                    ++i;
                }

            else if(y < rect.ymin())
                while(i < n && cs[i].y < rect.ymin()) {
                    ++i;
                }

            else if(y > rect.ymax())
                while(i < n && cs[i].y > rect.ymax()) {
                    ++i;
                }

            if(i >= n) {
                return false;
            }

            // Establish new position
            x = cs[i].x;
            y = cs[i].y;
            z = cs[i].z;
            pos = rect.position(x, y);

            // Handle all possible cases
            x0 = cs[i - 1].x;
            y0 = cs[i - 1].y;
            z0 = cs[i - 1].z;
            clip_to_edges(x0, y0, z0, x, y, z, rect);

            if(pos == Rectangle::Inside) {
                add_start = true;	// x0,y0 must have clipped the rectangle
                // Main loop will enter the Inside/Edge section

            }
            else if(pos == Rectangle::Outside) {
                // From Outside to Outside. We need to check whether
                // we created a line segment inside the box. In any
                // case, we will continue the main loop after this
                // which will then enter the Outside section.

                // Clip the other end too
                clip_to_edges(x, y, z, x0, y0, z0, rect);

                Rectangle::Position prev_pos = rect.position(x0, y0);
                pos = rect.position(x, y);

                if(different(x0, y0, x, y) &&		// discard corners etc
                        Rectangle::onEdge(prev_pos) &&		// discard if does not intersect rect
                        Rectangle::onEdge(pos) &&
                        !Rectangle::onSameEdge(prev_pos, pos)	// discard if travels along edge
                  ) {
                    auto coords = detail::make_unique<geom::CoordinateSequence>(2u);
                    coords->setAt(Coordinate(x0, y0, z0), 0);
                    coords->setAt(Coordinate(x, y, z), 1);
                    auto line = _gf->createLineString(std::move(coords));
                    parts.add(line.release());
                }

                // Continue main loop outside the rect

            }
            else {
                // From outside to edge. If the edge we hit first when
                // following the line is not the edge we end at, then
                // clearly we must go through the rectangle and hence
                // a start point must be set.

                Rectangle::Position newpos = rect.position(x0, y0);
                if(!Rectangle::onSameEdge(pos, newpos)) {
                    add_start = true;
                }
                else {
                    // we ignore the travel along the edge and continue
                    // the main loop at the last edge point
                }
            }
        }

        else {
            // The point is now strictly inside or on the edge.
            // Keep iterating until the end or the point goes
            // outside. We may have to output partial linestrings
            // while iterating until we go strictly outside

            auto start_index = i;			// 1st valid original point
            bool go_outside = false;

            while(!go_outside && ++i < n) {
                x = cs[i].x;
                y = cs[i].y;
                z = cs[i].z;

                Rectangle::Position prev_pos = pos;
                pos = rect.position(x, y);

                if(pos == Rectangle::Inside) {
                    // Just keep going
                }
                else if(pos == Rectangle::Outside) {
                    go_outside = true;

                    // Clip the outside point to edges
                    clip_to_edges(x, y, z, cs[i - 1].x, cs[i - 1].y, cs[i - 1].z, rect);
                    pos = rect.position(x, y);

                    // Does the line exit through the inside of the box?

                    bool through_box = (different(x, y, cs[i].x, cs[i].y) &&
                                        !Rectangle::onSameEdge(prev_pos, pos));

                    // Output a LineString if it at least one segment long

                    if(start_index < i - 1 || add_start || through_box) {
                        auto coords = detail::make_unique<CoordinateSequence>();
                        if(add_start) {
                            coords->add(Coordinate(x0, y0, z0));
                            add_start = false;
                        }
                        //line->addSubLineString(&g, start_index, i-1);
                        coords->add(cs.begin() + static_cast<long>(start_index),
                                    cs.begin() + static_cast<long>(i));

                        if(through_box) {
                            coords->add(Coordinate(x, y, z));
                        }

                        auto line = _gf->createLineString(std::move(coords));
                        parts.add(line.release());
                    }
                    // And continue main loop on the outside
                }
                else {
                    // on same edge?
                    if(Rectangle::onSameEdge(prev_pos, pos)) {
                        // Nothing to output if we haven't been elsewhere
                        if(start_index < i - 1 || add_start) {
                        auto coords = detail::make_unique<CoordinateSequence>();
                            //geom::LineString * line = new geom::LineString();
                            if(add_start) {
                                //line->addPoint(x0,y0);
                                coords->add(Coordinate(x0, y0, z0));
                                add_start = false;
                            }
                            //line->addSubLineString(&g, start_index, i-1);
                            coords->add(cs.begin() + static_cast<long>(start_index),
                                        cs.begin() + static_cast<long>(i));

                            auto line = _gf->createLineString(std::move(coords));
                            parts.add(line.release());
                        }
                        start_index = i;
                    }
                    else {
                        // On different edge. Must have gone through the box
                        // then - keep collecting points that generate inside
                        // line segments
                    }
                }
            }

            // Was everything in? If so, generate no output but return true in this case only
            if(start_index == 0 && i >= n) {
                return true;
            }

            // Flush the last line segment if data ended and there is something to flush

            if(!go_outside &&						// meaning data ended
                    (start_index < i - 1 || add_start)) {	// meaning something has to be generated
                auto coords = detail::make_unique<CoordinateSequence>();
                //geom::LineString * line = new geom::LineString();
                if(add_start) {
                    //line->addPoint(x0,y0);
                    coords->add(Coordinate(x0, y0, z0));
                    add_start = false;
                }
                //line->addSubLineString(&g, start_index, i-1);
                coords->add(cs.begin() + static_cast<long>(start_index),
                            cs.begin() + static_cast<long>(i));

                auto line = _gf->createLineString(std::move(coords));
                parts.add(line.release());
            }

        }
    }

    return false;

}

/**
 * \brief Clip polygon, do not close clipped ones
 */

void
RectangleIntersection::clip_polygon_to_linestrings(const geom::Polygon* g,
        RectangleIntersectionBuilder& toParts,
        const Rectangle& rect)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }

    // Clip the exterior first to see what's going on

    RectangleIntersectionBuilder parts(*_gf);

    // If everything was in, just clone the original

    if(clip_linestring_parts(g->getExteriorRing(), parts, rect)) {
        toParts.add(g->clone().release());
        return;
    }

    // Now, if parts is empty, our rectangle may be inside the polygon
    // If not, holes are outside too

    if(parts.empty()) {
        // We could now check whether the rectangle is inside the outer
        // ring to avoid checking the holes. However, if holes are much
        // smaller than the exterior ring just checking the holes
        // separately could be faster.

        if(g->getNumInteriorRing() == 0) {
            return;
        }

    }
    else {
        // The exterior must have been clipped into linestrings.
        // Move them to the actual parts collector, clearing parts.

        parts.reconnect();
        parts.release(toParts);
    }

    // Handle the holes now:
    // - Clipped ones become linestrings
    // - Intact ones become new polygons without holes

    for(std::size_t i = 0, n = g->getNumInteriorRing(); i < n; ++i) {
        if(clip_linestring_parts(g->getInteriorRingN(i), parts, rect)) {
            // clones
            auto hole =g->getInteriorRingN(i)->clone();
            // becomes exterior
            Polygon* poly = _gf->createPolygon(std::move(hole)).release();
            toParts.add(poly);
        }
        else if(!parts.empty()) {
            parts.reconnect();
            parts.release(toParts);
        }
    }
}

/**
 * \brief Clip polygon, close clipped ones
 */

void
RectangleIntersection::clip_polygon_to_polygons(const geom::Polygon* g,
        RectangleIntersectionBuilder& toParts,
        const Rectangle& rect)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }

    // Clip the exterior first to see what's going on

    RectangleIntersectionBuilder parts(*_gf);

    // If everything was in, just clone the original

    const LineString* shell = g->getExteriorRing();
    if(clip_linestring_parts(shell, parts, rect)) {
        toParts.add(g->clone().release());
        return;
    }

    // If there were no intersections, the outer ring might be
    // completely outside.

    using geos::algorithm::Orientation;
    if(parts.empty()) {
        Coordinate rectCenter(rect.xmin(), rect.ymin());
        rectCenter.x += (rect.xmax() - rect.xmin()) / 2;
        rectCenter.y += (rect.ymax() - rect.ymin()) / 2;
        if(PointLocation::locateInRing(rectCenter,
                                       *g->getExteriorRing()->getCoordinatesRO())
                != Location::INTERIOR) {
            return;
        }
    }
    else {
        // TODO: make CCW checking part of clip_linestring_parts ?
        if(Orientation::isCCW(shell->getCoordinatesRO())) {
            parts.reverseLines();
        }
    }

    // Must do this to make sure all end points are on the edges

    parts.reconnect();

    // Handle the holes now:
    // - Clipped ones become part of the exterior
    // - Intact ones become holes in new polygons formed by exterior parts


    for(std::size_t i = 0, n = g->getNumInteriorRing(); i < n; ++i) {
        RectangleIntersectionBuilder holeparts(*_gf);
        const LinearRing* hole = g->getInteriorRingN(i);
        if(clip_linestring_parts(hole, holeparts, rect)) {
            // becomes exterior
            auto cloned = hole->clone();
            Polygon* poly = _gf->createPolygon(std::move(cloned)).release();
            parts.add(poly);
        }
        else {
            if(!holeparts.empty()) {
                // TODO: make CCW checking part of clip_linestring_parts ?
                if(! Orientation::isCCW(hole->getCoordinatesRO())) {
                    holeparts.reverseLines();
                }
                holeparts.reconnect();
                holeparts.release(parts);
            }
            else {

                Coordinate rectCenter(rect.xmin(), rect.ymin());
                rectCenter.x += (rect.xmax() - rect.xmin()) / 2;
                rectCenter.y += (rect.ymax() - rect.ymin()) / 2;
                if(PointLocation::isInRing(rectCenter,
                                           g->getInteriorRingN(i)->getCoordinatesRO())) {
                    // Completely inside the hole
                    return;
                }
            }
        }
    }

    parts.reconnectPolygons(rect);
    parts.release(toParts);
}

/**
 * \brief Clip  geometry
 */

void
RectangleIntersection::clip_polygon(const geom::Polygon* g,
                                    RectangleIntersectionBuilder& parts,
                                    const Rectangle& rect,
                                    bool keep_polygons)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }

    if(keep_polygons) {
        clip_polygon_to_polygons(g, parts, rect);
    }
    else {
        clip_polygon_to_linestrings(g, parts, rect);
    }
}

/**
 * \brief Clip geometry
 */

void
RectangleIntersection::clip_linestring(const geom::LineString* g,
                                       RectangleIntersectionBuilder& parts,
                                       const Rectangle& rect)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }

    // If everything was in, just clone the original

    if(clip_linestring_parts(g, parts, rect)) {
        parts.add(g->clone().release());
    }

}

void
RectangleIntersection::clip_multipoint(const geom::MultiPoint* g,
                                       RectangleIntersectionBuilder& parts,
                                       const Rectangle& rect)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }
    for(std::size_t i = 0, n = g->getNumGeometries(); i < n; ++i) {
        clip_point(g->getGeometryN(i), parts, rect);
    }
}

void
RectangleIntersection::clip_multilinestring(const geom::MultiLineString* g,
        RectangleIntersectionBuilder& parts,
        const Rectangle& rect)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }

    for(std::size_t i = 0, n = g->getNumGeometries(); i < n; ++i) {
        clip_linestring(g->getGeometryN(i), parts, rect);
    }
}

void
RectangleIntersection::clip_multipolygon(const geom::MultiPolygon* g,
        RectangleIntersectionBuilder& parts,
        const Rectangle& rect,
        bool keep_polygons)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }

    for(std::size_t i = 0, n = g->getNumGeometries(); i < n; ++i) {
        clip_polygon(g->getGeometryN(i), parts, rect, keep_polygons);
    }
}

void
RectangleIntersection::clip_geometrycollection(
    const geom::GeometryCollection* g,
    RectangleIntersectionBuilder& parts,
    const Rectangle& rect,
    bool keep_polygons)
{
    if(g == nullptr || g->isEmpty()) {
        return;
    }

    for(std::size_t i = 0, n = g->getNumGeometries(); i < n; ++i) {
        clip_geom(g->getGeometryN(i),
                  parts, rect, keep_polygons);
    }
}

void
RectangleIntersection::clip_geom(const geom::Geometry* g,
                                 RectangleIntersectionBuilder& parts,
                                 const Rectangle& rect,
                                 bool keep_polygons)
{
    if(const Point* p1 = dynamic_cast<const geom::Point*>(g)) {
        return clip_point(p1, parts, rect);
    }
    else if(const MultiPoint* p2 = dynamic_cast<const geom::MultiPoint*>(g)) {
        return clip_multipoint(p2, parts, rect);
    }
    else if(const LineString* p3 = dynamic_cast<const geom::LineString*>(g)) {
        return clip_linestring(p3, parts, rect);
    }
    else if(const MultiLineString* p4 = dynamic_cast<const geom::MultiLineString*>(g)) {
        return clip_multilinestring(p4, parts, rect);
    }
    else if(const Polygon* p5 = dynamic_cast<const geom::Polygon*>(g)) {
        return clip_polygon(p5, parts, rect, keep_polygons);
    }
    else if(const MultiPolygon* p6 = dynamic_cast<const geom::MultiPolygon*>(g)) {
        return clip_multipolygon(p6, parts, rect, keep_polygons);
    }
    else if(const GeometryCollection* p7 = dynamic_cast<const geom::GeometryCollection*>(g)) {
        return clip_geometrycollection(p7, parts, rect, keep_polygons);
    }
    else {
        throw util::UnsupportedOperationException("Encountered an unknown geometry component when clipping polygons");
    }
}

/* public static */
std::unique_ptr<geom::Geometry>
RectangleIntersection::clipBoundary(const geom::Geometry& g, const Rectangle& rect)
{
    RectangleIntersection ri(g, rect);
    return ri.clipBoundary();
}

std::unique_ptr<geom::Geometry>
RectangleIntersection::clipBoundary()
{
    RectangleIntersectionBuilder parts(*_gf);

    bool keep_polygons = false;
    clip_geom(&_geom, parts, _rect, keep_polygons);

    return parts.build();
}

/* public static */
std::unique_ptr<geom::Geometry>
RectangleIntersection::clip(const geom::Geometry& g, const Rectangle& rect)
{
    RectangleIntersection ri(g, rect);
    std::unique_ptr<geom::Geometry> result = ri.clip();

    if (g.hasZ()) {
        std::unique_ptr<ElevationModel> elevModel = ElevationModel::create(g);
        elevModel->populateZ(*result);
    }

    return result;
}

std::unique_ptr<geom::Geometry>
RectangleIntersection::clip()
{
    RectangleIntersectionBuilder parts(*_gf);

    bool keep_polygons = true;
    clip_geom(&_geom, parts, _rect, keep_polygons);

    return parts.build();
}

RectangleIntersection::RectangleIntersection(const geom::Geometry& geom, const Rectangle& rect)
    : _geom(geom), _rect(rect),
      _gf(geom.getFactory())
{
}

} // namespace geos::operation::intersection
} // namespace geos::operation
} // namespace geos