#include "PlanimetryCalculus.h"
#include "Circle.h"


namespace PlanimetryCalculus {
    namespace  {
        /*!
            \fn IntersectType intersects(const Circle &circle, const QLineF &line, QList<QPointF> &intersectionPoints, bool calcInstersect)
            Returns the Intersection type of \a circle and \a line.
            Stores the intersection points in \a intersectionPoints if \a calcIntersect is \c true.
            Returns \c Error if either line or circe \c {isNull() == true}.

            \sa QPointF, Circle
        */
        IntersectType intersects(const Circle &circle, const QLineF &line, QList<QPointF> &intersectionPoints, bool calcInstersect)
        {
            if (!circle.isNull() && ! line.isNull()) {
                QPointF translationVector = line.p1();
                double angleWLDegree = line.angle(); // angle between wold and line coordinate system

                QPointF originCircleL = circle.origin() - translationVector;
                rotateReference(originCircleL, -angleWLDegree); // circle origin in line corrdinate system

                double y = originCircleL.y();
                double r = circle.radius();
                if (qAbs(y) > r)
                    return NoIntersection;
                else if ( qFuzzyCompare(qFabs(y), r) ) { // tangent
                    double x_ori = originCircleL.x();

                    if (x_ori >= 0 && x_ori <= line.length()) {
                        if (calcInstersect) {
                            QPointF intersectionPt = QPointF(x_ori, 0);
                            rotateReference(intersectionPt, angleWLDegree);
                            intersectionPoints.append(intersectionPt + translationVector);
                        }

                        return Tangent;
                    }

                    return NoIntersection;
                } else { // sekant
                    double x_ori = originCircleL.x();
                    double y_ori = originCircleL.y();
                    double delta = qSqrt(qPow(r, 2)-qPow(y_ori, 2));
                    double x1    = x_ori + delta; // x coordinate (line system) of fist intersection point
                    double x2    = x_ori - delta;// x coordinate (line system) of second intersection point
                    bool doesIntersect = false; // remember if actual intersection was on the line

                    if (x1 >= 0 && x1 <= line.length()) { // check if intersection point is on the line
                        if (calcInstersect) {
                            QPointF intersectionPt = QPointF(x1, 0); // first intersection point (line system)
                            rotateReference(intersectionPt, angleWLDegree);
                            intersectionPoints.append(intersectionPt + translationVector); // transform (to world system) and append first intersection point
                        }
                        doesIntersect = true;
                    }
                    if (x2 >= 0 && x2 <= line.length()) { // check if intersection point is on the line
                        if (calcInstersect) {
                            QPointF intersectionPt = QPointF(x2, 0); // second intersection point (line system)
                            rotateReference(intersectionPt, angleWLDegree);
                            intersectionPoints.append(intersectionPt + translationVector); // transform (to world system) and append second intersection point
                        }
                        doesIntersect = true;
                    }

                    return doesIntersect ? Secant : NoIntersection;
                }
            }

            return Error;
        }
    } // end anonymous namespace

    /*!
        \fn void rotatePoint(QPointF &point, double alpha)
        Rotates the \a point counter clockwisely by the angle \a alpha (in radiants).
    */
    void rotateReference(QPointF &point, double alpha)
    {
        if (!point.isNull()) {
            double x = point.x();
            double y = point.y();

            point.setX(x*qCos(alpha) - y*qSin(alpha));
            point.setY(x*qSin(alpha) + y*qCos(alpha));
        }
    }

    void rotateReference(QList<QPointF> &points, double alpha)
    {
        for (int i = 0; i < points.size(); i++) {
            rotateReference(points[i], alpha);
        }
    }

    /*!
        \fn void rotatePointDegree(QPointF &point, double alpha)
        Rotates the \a point counter clockwisely by the angle \a alpha (in degrees).
    */
    void rotateReferenceDegree(QPointF &point, double alpha)
    {
        rotateReference(point, alpha/180*M_PI);
    }

    void rotateReferenceDegree(QList<QPointF> &points, double alpha)
    {
        for (int i = 0; i < points.size(); i++) {
            rotateReferenceDegree(points[i], alpha);
        }
    }

    /*!
        \fn IntersectType intersects(const Circle &circle1, const Circle &circle2)
        Returns the intersection type of the two cirles \a circle1 and \a circle2.

        \note Returns Error if circle.isNull() returns true;

        \sa Circle
    */
    IntersectType intersects(const Circle &circle1, const Circle &circle2)
    {
        // r1 == 0 || r2 == 0 results in indefined behavior
        if (!circle1.isNull() && !circle2.isNull()) {
            double r1 = circle1.radius();
            double r2 = circle2.radius();
            double d  = distance(circle1.origin(), circle2.origin());
            double r = 0;
            double R = 0;
            if (r1 > r2) {
                R = r1; // large
                r = r2; // small
            } else {
                // this branch is also choosen if r1 == r2
                R = r2;
                r = r1;
            }

            if        (r + d < R) {
                // this branch is also reached if d < rLarge && rSmall == 0
                return InsideNoIntersection;
            } else if (qFuzzyCompare(r + d, R)) {
                if    (qFuzzyIsNull(d))
                    return CirclesEqual;
                else
                    return InsideTouching;
            } else if (d < R) {
                return InsideIntersection;
            } else if (d - r < R) {
                return OutsideIntersection;
            } else if (qFuzzyCompare(d - r, R)) {
                return OutsideTouching;
            } else {
                return OutsideNoIntersection;
            }
        }
        return Error;
    }

    /*!
        \fn IntersectType intersects(const Circle &circle1, const Circle &circle2, QList<QPointF> intersectionPoints)
        Calculates the intersection points of two circles if present and stores the result in \a intersectionPoints.
        Returns the intersection type of the two cirles \a circle1 and \a circle2.

        The function assumes that the list \a intersectionPoints is empty.

        \note Returns Error if circle.isNull() returns true;

        \sa Circle
    */
    IntersectType intersects(const Circle &circle1, const Circle &circle2, QList<QPointF> &intersectionPoints)
    {
        IntersectType returnValue = intersects(circle1, circle2);

        if (   returnValue == InsideNoIntersection
            || returnValue == OutsideNoIntersection
            || returnValue == CirclesEqual
            || returnValue == Error  ) {
            return returnValue; // No intersection Points, or infinitly many (in case of CirclesEqual).
        } else {
            double r1    = circle1.radius();
            double r2    = circle2.radius();
            double d     = distance(circle1.origin(), circle2.origin());
            double alpha = angle(circle1.origin(), circle2.origin());
            double r = 0;
            double R = 0;
            if (r1 > r2) {
                R = r1;
                r = r2;
            } else {
                // this branch is also choosen if r1 == r2
                R = r2;
                r = r1;
            }

            if        (   returnValue == InsideTouching
                       || returnValue == OutsideTouching) {
                // Intersection point in coordinate system of circle 1.
                // Coordinate system circle1: origin = circle1.origin(), x-axis towars circle2.origin() y-axis such that the
                // coordinate system is dextrorse with z-axis outward faceing with respect to the drawing plane.
                intersectionPoints.append(QPointF(0, r1));
            } else { //triggered if (   returnValue == InsideIntersection
                     //              || returnValue == OutsideIntersection)
                // See fist branch for explanation

                // this equations are obtained by solving x^2+y^2=R^2 and (x - d)^2+y^2=r^2
                double x  = (qPow(d, 2) - qPow(r, 2) + qPow(R, 2))/2/d;
                double y = 1/2/d*qSqrt(4*qPow(d*R, 2) - qPow(qPow(d, 2) - qPow(r, 2) +  qPow(R, 2), 2));

                intersectionPoints.append(QPointF(x, y));
                intersectionPoints.append(QPointF(x, -y));
            }
            // Transform the coordinate to the world coordinate system. Alpha is the angle between world and circle1 coordinate system.
            rotateReference(intersectionPoints, alpha);

            return returnValue;
        }
    }


    /*!
        \fn IntersectType intersects(const Circle &circle, const QLineF &line)
        Returns the Intersection type of \a circle and \a line.
        Returns \c Error if either line or circe \c {isNull() == true}.

        \sa QPointF, Circle
    */
    IntersectType intersects(const Circle &circle, const QLineF &line)
    {
        QList<QPointF> dummyList;
        return intersects(circle, line, dummyList, false /* calculate intersection points*/);
    }

    IntersectType intersects(const Circle &circle, const QLineF &line, QList<QPointF> &intersectionPoints)
    {
        return intersects(circle, line, intersectionPoints, true /* calculate intersection points*/);
    }

    /*!
        \fn double distance(const QPointF &p1, const QPointF p2)
        Calculates the distance (2-norm) between \a p1 and \a p2.
        \sa QPointF
    */
    double distance(const QPointF &p1, const QPointF p2)
    {
        double dx = p2.x()-p1.x();
        double dy = p2.y()-p1.y();

        return qSqrt(dx*dx+dy*dy);
    }

    /*!
        \fn double distance(const QPointF &p1, const QPointF p2)
        Calculates the angle (in radiants) between the line defined by \a p1 and \a p2 and the x-axis according to the following rule.
        Angle = qAtan2(dy, dx), where dx = p2.x()-p1.x() and dy = p2.y()-p1.y().

        \note The order of \a p1 and \a p2 matters. Swapping \a p1 and \a p2 will result in a angle of oposite sign.
        \sa QPointF
    */
    double angle(const QPointF &p1, const QPointF p2)
    {
        double dx = p2.x()-p1.x();
        double dy = p2.y()-p1.y();

        return qAtan2(dy, dx);
    }

    /*!
        \fn double distance(const QPointF &p1, const QPointF p2)
        Calculates the angle (in degrees) between the line defined by \a p1 and \a p2 and the x-axis according to the following rule.
        Angle = qAtan2(dy, dx)*180/pi, where dx = p2.x()-p1.x() and dy = p2.y()-p1.y().

        \note The order of \a p1 and \a p2 matters. Swapping \a p1 and \a p2 will result in a angle of oposite sign.
        \sa QPointF
    */
    double angleDegree(const QPointF &p1, const QPointF p2)
    {
        return angle(p1, p2)*180/M_PI;
    }

    double truncateAngle(double angle)
    {
        while (angle < 0     ) { angle += 2*M_PI;}
        while (angle > 2*M_PI) { angle -= 2*M_PI;}

        return angle;
    }

    double truncateAngleDegree(double angle)
    {
        return truncateAngle(angle/180*M_PI);
    }


    /*!
     * \fn IntersectType intersects(const QLineF &line1, const QLineF &line2, QPointF &intersectionPt)
     * Determines wheter \a line1 and \a line2 intersect and of which type the intersection is.
     * Stores the intersection point in \a intersectionPt
     * Returns the intersection type (\c IntersectType).
     *
     * Intersection Types:
     * \c NoIntersection
     * \c CornerCornerIntersection; A intersection is present such that two of the lines cornes touch.
     * \c EdgeCornerIntersection; A intersection is present such that a corner and a edge touch.
     * \c EdgeEdgeIntersection; A intersection is present such two edges intersect.
     * \c LinesParallel
     * \c LinesEqual
     * \c Error; Returned if at least on line delivers isNULL() == true.
     *
     *
     * \sa QGeoCoordinate
     */
    IntersectType intersects(const QLineF &line1, const QLineF &line2, QPointF &intersectionPt)
    {
        if (line1.isNull() || line2.isNull())
            return Error;

        // line 1 coordinate system: origin line1.p1(), x-axis towards line1.p2()
        QPointF translationVector = line1.p1(); // translation vector between world and line1 system
        double alpha = line1.angle();
        double l1 = line1.length();

        QLineF line2L1 = line1;
        line2L1.translate(-translationVector);
        rotateReference(line2L1, -alpha);

        double x1 = line2L1.x1();
        double x2 = line2L1.x2();
        double y1 = line2L1.y1();
        double y2 = line2L1.y2();
        double dx = (x2-x1);
        double dy = (y2-y1);
        double xNull = 0; // (xNull, 0) intersection point in line1 system

        if (!qFuzzyIsNull(dx)) {
            double k = dy/dx;

            if (qFuzzyIsNull(k)) {
                if (qFuzzyIsNull(x1) && qFuzzyIsNull(y1) && qFuzzyCompare(x2, l1))
                    return LinesEqual;
                else
                    return LinesParallel;
            }


            double d = (y1*x2-y2*x1)/dx;
            xNull = -d/k;
        } else {
            // lines orthogonal
            if (signum(y1) != signum(y2)){
                xNull = x1;
            }
            else
                return NoIntersection;
        }

        IntersectType returnValue;

        if (xNull >= x1 && xNull <= x2){
            // determine intersection type

            if(qFuzzyIsNull(xNull) || qFuzzyCompare(xNull, l1)) {
                if(qFuzzyIsNull(y1) || qFuzzyIsNull(y2))
                    returnValue = CornerCornerIntersection;
                else
                    returnValue = EdgeCornerIntersection;
            } else if (xNull > 0 && xNull < l1) {
                if(qFuzzyIsNull(y1) || qFuzzyIsNull(y2))
                    returnValue = EdgeCornerIntersection;
                else
                    returnValue = EdgeEdgeIntersection;
            } else
                return NoIntersection;
        } else
            return NoIntersection;

        intersectionPt = QPointF(xNull, 0); // intersection point in line1 system
        rotateReference(intersectionPt, alpha);
        intersectionPt += translationVector;

        return returnValue;
    }

    /*!
     * \overload QList<IntersectionType> intersects(const QList<QGeoCoordinate> &polygon, const QList<QGeoCoordinate> &line, QList<QGeoCoordinate> &intersectionList, QList<QPair<int, int> > &neighbourList)
     * Checks if \a polygon intersect with \a line.
     * Stores the intersection points in \a intersectionList.
     *
     * Stores the indices of the closest two \a area vetices for each of coorespoinding intersection points in \a neighbourList. *
     * For example if an intersection point is found between the first and the second vertex of the \a area the intersection point will
     * be stored in \a intersectionList and the indices 1 and 2 will be stored in \a neighbourList.
     * \a neighbourList has entries of type \c {QPair<int, int>}, where \c{pair.first} would contain 1 and \c{pair.second} would contain 2, when
     * relating to the above example.
     *
     * Returns the \c IntersectionType of each intersection point within a QList.
     *
     * \sa QPair, QList
     */
    QList<IntersectionType> intersects(const QList<QGeoCoordinate> &polygon, const Line &line, QList<QGeoCoordinate> &intersectionList, QList<QPair<int, int> > &neighbourList)
    {
        if (polygon.size() >= 3) { // are line a proper line and poly a proper poly?other,
            intersectionList.clear();
            neighbourList.clear();
            QList<IntersectionType> intersectionTypeList;
            // Assemble a line form each tow consecutive polygon vertices and check whether it intersects with line
            for (int i = 0; i < polygon.size(); i++) {

                Line interatorLine;
                QGeoCoordinate currentVertex    = polygon[i];
                QGeoCoordinate nextVertex       = polygon[nextPolygonIndex(polygon.size(), i)];
                interatorLine.first             = currentVertex;
                interatorLine.second            = nextVertex;

                QGeoCoordinate intersectionPoint;
                IntersectionType returnValue =  intersects(line, interatorLine, intersectionPoint);
                if (    returnValue == IntersectionType::EdgeIntersection
                     || returnValue == IntersectionType::InteriorIntersection) {
                    intersectionList.append(intersectionPoint);

                    QPair<int, int>     neighbours;
                    neighbours.first    = i;
                    neighbours.second   = nextPolygonIndex(polygon.size(), i);
                    neighbourList.append(neighbours);

                    intersectionTypeList.append(returnValue);
                }
            }

            if (intersectionList.count() > 0) {
                return intersectionTypeList;
            } else {
                return QList<IntersectionType>();
            }
        } else {
            qWarning("WimaArea::intersects(line, poly): line->count() != 2 || poly->count() < 3");
            return QList<IntersectionType>();
        }
    }

    /*!
     * \overload bool intersects(const QList<QGeoCoordinate> &polygon, const QList<QGeoCoordinate> &line)
     * Returns \c true if any intersection between \a polygon and \a line exists, \c false else.
     *
     * \sa QPair, QList
     */
    bool intersects(const QList<QGeoCoordinate> &polygon, const Line &line)
    {
        QList<QGeoCoordinate> dummyGeo;
        QList<QPair<int, int>> dummyNeighbour;
        intersects(polygon, line, dummyGeo, dummyNeighbour);

        if (dummyGeo.size() > 0)
            return true;

        return false;
    }

    void rotateReference(QLineF &line, double alpha)
    {
        line.setP1(rotateReturn(line.p1(), alpha));
        line.setP2(rotateReturn(line.p2(), alpha));
    }

    QPointF rotateReturn(QPointF point, double alpha)
    {
        rotateReference(point, alpha);
        return point;
    }

    QList<QPointF> rotateReturn(QList<QPointF> points, double alpha)
    {
        rotateReference(points, alpha);
        return points;
    }

    QLineF rotateReturn(QLineF line, double alpha)
    {
        rotateReference(line, alpha);
        return line;
    }

    /*!
     * \fntemplate <typename T> int signum(T val)
     * Returns the signum of a value \a val.
     *
     * \sa QPair, QList
     */
    template <typename T> int signum(T val)
    {
        return (T(0) < val) - (val < T(0));
    }

} // end namespace PlanimetryCalculus



/*!
    \class PlanimetryCalculus
    \inmodule Wima

    \brief The \c PlanimetryCalculus class provides routines handy for planimetrical (2D) calculations.
*/