/**********************************************************************

 *

 * GEOS - Geometry Engine Open Source

 * http://geos.osgeo.org

 *

 * Copyright (C) 2009 2011 Sandro Santilli <strk@kbt.io>

 * Copyright (C) 2005-2006 Refractions Research Inc.

 * Copyright (C) 2001-2002 Vivid Solutions Inc.

 *

 * 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: geom/LineSegment.java r18 (JTS-1.11)

 *

 **********************************************************************/

#include <geos/constants.h>

#include <geos/geom/LineSegment.h>

#include <geos/geom/LineString.h> // for toGeometry

#include <geos/geom/CoordinateSequence.h>

#include <geos/geom/GeometryFactory.h>

#include <geos/algorithm/LineIntersector.h>

#include <geos/algorithm/Intersection.h>

#include <geos/util/IllegalStateException.h>

#include <geos/profiler.h>

#include <geos/util.h>

#include <algorithm> // for max

#include <sstream>

#include <cmath>

namespace geos {

namespace geom { // geos::geom

/*public*/

void

LineSegment::reverse()

{

    std::swap(p0, p1);

}

double

LineSegment::projectionFactor(const CoordinateXY& p) const

{

    return projectionFactor(p0, p1, p);

}

/*public*/

double

LineSegment::projectionFactor(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q)

{

    if(q == p0) {

        return 0.0;

    }

    if(q == p1) {

        return 1.0;

    }

    if(p0 == p1) {

        return 0.0;

    }

    // Otherwise, use comp.graphics.algorithms Frequently Asked Questions method

    /*(1)             AC dot AB

                   r = ---------

                         ||AB||^2

                r has the following meaning:

                r=0 P = A

                r=1 P = B

                r<0 P is on the backward extension of AB

                r>1 P is on the forward extension of AB

                0<r<1 P is interior to AB

        */

    double dx = p1.x - p0.x;

    double dy = p1.y - p0.y;

    double len2 = dx * dx + dy * dy;

    double r = ((q.x - p0.x) * dx + (q.y - p0.y) * dy) / len2;

    return r;

}

/*public*/

double

LineSegment::segmentFraction(const CoordinateXY& inputPt) const

{

    double segFrac = projectionFactor(inputPt);

    if(segFrac < 0.0) {

        segFrac = 0.0;

    }

    else if(segFrac > 1.0) {

        segFrac = 1.0;

    }

    return segFrac;

}

/*public*/

void

LineSegment::project(const Coordinate& p, Coordinate& ret) const

{

    if(p == p0 || p == p1) {

        ret = p;

        return;

    }

    double r = projectionFactor(p);

    ret = Coordinate(p0.x + r * (p1.x - p0.x), p0.y + r * (p1.y - p0.y));

}

CoordinateXY

LineSegment::project(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q)

{

    if(q == p0 || q == p1) {

        return q;

    }

    double r = projectionFactor(p0, p1, q);

    double x = p0.x + r * (p1.x - p0.x);

    double y = p0.y + r * (p1.y - p0.y);

    return CoordinateXY(x, y);

}

CoordinateXY

LineSegment::project(const CoordinateXY& p) const

{

    return project(p0, p1, p);

}

/*private*/

CoordinateXY

LineSegment::project(const CoordinateXY& p0, const CoordinateXY& p1, double factor)

{

    if(factor == 1.0) {

        return p1;

    }

    return CoordinateXY(p0.x + factor * (p1.x - p0.x), p0.y + factor * (p1.y - p0.y));

}

void

LineSegment::project(double factor, CoordinateXY& ret) const

{

    ret = project(p0, p1, factor);

}

bool

LineSegment::project(const LineSegment& seg, LineSegment& ret) const

{

    double pf0 = projectionFactor(seg.p0);

    double pf1 = projectionFactor(seg.p1);

    // check if segment projects at all

    if(pf0 >= 1.0 && pf1 >= 1.0) {

        return false;

    }

    if(pf0 <= 0.0 && pf1 <= 0.0) {

        return false;

    }

    Coordinate newp0;

    project(pf0, newp0);

    if (pf0 < 0.0) newp0 = p0;

    if (pf0 > 1.0) newp0 = p1;

    Coordinate newp1;

    project(pf1, newp1);

    if (pf1 < 0.0) newp1 = p0;

    if (pf1 > 1.0) newp1 = p1;

    ret.setCoordinates(newp0, newp1);

    return true;

}

CoordinateXY

LineSegment::closestPoint(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q)

{

    double factor = projectionFactor(p0, p1, q);

    if(factor > 0 && factor < 1) {

        return project(p0, p1, factor);

    }

    double dist0 = p0.distance(q);

    double dist1 = p1.distance(q);

    if(dist0 < dist1) {

        return p0;

    }

    return p1;

}

//Coordinate*

void

LineSegment::closestPoint(const CoordinateXY& p, CoordinateXY& ret) const

{

    double factor = projectionFactor(p);

    if(factor > 0 && factor < 1) {

        project(factor, ret);

        return;

    }

    double dist0 = p0.distance(p);

    double dist1 = p1.distance(p);

    if(dist0 < dist1) {

        ret = p0;

        return;

    }

    ret = p1;

}

/*public*/

bool

LineSegment::equalsTopo(const LineSegment& other) const

{

    return (p0 == other.p0 && p1 == other.p1) || (p0 == other.p1 && p1 == other.p0);

}

/*public*/

int

LineSegment::orientationIndex(const LineSegment& seg) const

{

    int orient0 = algorithm::Orientation::index(p0, p1, seg.p0);

    int orient1 = algorithm::Orientation::index(p0, p1, seg.p1);

    // this handles the case where the points are L or collinear

    if(orient0 >= 0 && orient1 >= 0) {

        return std::max(orient0, orient1);

    }

    // this handles the case where the points are R or collinear

    if(orient0 <= 0 && orient1 <= 0) {

        return std::min(orient0, orient1);

    }

    // points lie on opposite sides ==> indeterminate orientation

    return 0;

}

std::array<Coordinate, 2>

LineSegment::closestPoints(const LineSegment& line)

{

    // test for intersection

    Coordinate intPt = intersection(line);

    if(!intPt.isNull()) {

        return {{ intPt, intPt }};

    }

    /*

     * if no intersection closest pair contains at least one endpoint.

     * Test each endpoint in turn.

     */

    std::array<Coordinate, 2> closestPt;

    double dist;

    Coordinate close00;

    closestPoint(line.p0, close00);

    double minDistance = close00.distance(line.p0);

    closestPt[0] = close00;

    closestPt[1] = line.p0;

    Coordinate close01;

    closestPoint(line.p1, close01);

    dist = close01.distance(line.p1);

    if(dist < minDistance) {

        minDistance = dist;

        closestPt[0] = close01;

        closestPt[1] = line.p1;

    }

    Coordinate close10;

    line.closestPoint(p0, close10);

    dist = close10.distance(p0);

    if(dist < minDistance) {

        minDistance = dist;

        closestPt[0] = p0;

        closestPt[1] = close10;

    }

    Coordinate close11;

    line.closestPoint(p1, close11);

    dist = close11.distance(p1);

    if(dist < minDistance) {

        closestPt[0] = p1;

        closestPt[1] = close11;

    }

    return closestPt;

}

std::array<CoordinateXY, 2>

LineSegment::closestPoints(const CoordinateXY& p0, const CoordinateXY& p1,

                           const CoordinateXY& q0, const CoordinateXY& q1)

{

    // test for intersection

    algorithm::LineIntersector li;

    li.computeIntersection(p0, p1, q0, q1);

    if(li.hasIntersection()) {

        return {CoordinateXY{li.getIntersection(0)},

                   CoordinateXY{li.getIntersection(0)}};

    }

    /*

     * if no intersection closest pair contains at least one endpoint.

     * Test each endpoint in turn.

     */

    std::array<CoordinateXY, 2> closestPt;

    CoordinateXY close00 = closestPoint(p0, p1, q0);

    double minDistance = close00.distance(q0);

    closestPt[0] = close00;

    closestPt[1] = q0;

    CoordinateXY close01 = closestPoint(p0, p1, q1);

    double dist = close01.distance(q1);

    if(dist < minDistance) {

        minDistance = dist;

        closestPt[0] = close01;

        closestPt[1] = q1;

    }

    CoordinateXY close10 = closestPoint(q0, q1, p0);

    dist = close10.distance(p0);

    if(dist < minDistance) {

        minDistance = dist;

        closestPt[0] = p0;

        closestPt[1] = close10;

    }

    CoordinateXY close11 = closestPoint(q0, q1, p1);

    dist = close11.distance(p1);

    if(dist < minDistance) {

        closestPt[0] = p1;

        closestPt[1] = close11;

    }

    return closestPt;

}

Coordinate

LineSegment::intersection(const LineSegment& line) const

{

    algorithm::LineIntersector li;

    li.computeIntersection(p0, p1, line.p0, line.p1);

    if(li.hasIntersection()) {

        return li.getIntersection(0);

    }

    Coordinate rv;

    rv.setNull();

    return rv;

}

Coordinate

LineSegment::lineIntersection(const LineSegment& line) const

{

    // TODO return a CoordinateXY here.

    return Coordinate(algorithm::Intersection::intersection(p0, p1, line.p0, line.p1));

}

/* public */

void

LineSegment::pointAlongOffset(double segmentLengthFraction,

                              double offsetDistance,

                              Coordinate& ret) const

{

    // the point on the segment line

    double segx = p0.x + segmentLengthFraction * (p1.x - p0.x);

    double segy = p0.y + segmentLengthFraction * (p1.y - p0.y);

    double dx = p1.x - p0.x;

    double dy = p1.y - p0.y;

    double len = std::sqrt(dx * dx + dy * dy);

    double ux = 0.0;

    double uy = 0.0;

    if(offsetDistance != 0.0) {

        if(len <= 0.0) {

            throw util::IllegalStateException("Cannot compute offset from zero-length line segment");

        }

        // u is the vector that is the length of the offset,

        // in the direction of the segment

        ux = offsetDistance * dx / len;

        uy = offsetDistance * dy / len;

    }

    // the offset point is the seg point plus the offset

    // vector rotated 90 degrees CCW

    double offsetx = segx - uy;

    double offsety = segy + ux;

    ret = Coordinate(offsetx, offsety);

}

/* public */

LineSegment

LineSegment::offset(double offsetDistance)

{

    Coordinate offset0, offset1;

    pointAlongOffset(0, offsetDistance, offset0);

    pointAlongOffset(1, offsetDistance, offset1);

    LineSegment ls(offset0, offset1);

    return ls;

}

/* public */

double

LineSegment::distancePerpendicularOriented(const CoordinateXY& p) const

{

    if (p0.equals2D(p1))

        return p0.distance(p);

    double dist = distancePerpendicular(p);

    if (orientationIndex(p) < 0)

        return -dist;

    return dist;

}

/* public */

std::unique_ptr<LineString>

LineSegment::toGeometry(const GeometryFactory& gf) const

{

    auto cl = detail::make_unique<CoordinateSequence>(2u);

    cl->setAt(p0, 0);

    cl->setAt(p1, 1);

    return gf.createLineString(std::move(cl));

}

/* public */

std::ostream&

LineSegment::operator<< (std::ostream& os)

{

    return os << "LINESEGMENT("

        << p0.x << " " << p0.y << ","

        << p1.x << " " << p1.y << ")";

}

} // namespace geos::geom

} // namespace geos