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

 Copyright (C) 2010-2015 the Office National des Forêts (ONF), France
                     All rights reserved.

 Contact : alexandre.piboule@onf.fr

 Developers : Alexandre Piboule (ONF)

 This file is part of PluginShared library 2.0.

 PluginShared is free library: you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 PluginShared is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU Lesser General Public License
 along with PluginShared.  If not, see <http://www.gnu.org/licenses/lgpl.html>.

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

#include "ct_delaunaytriangle.h"

#include <stdlib.h>     /* srand, rand */

CT_DelaunayTriangle::CT_DelaunayTriangle(CT_DelaunayVertex *v1t, CT_DelaunayVertex *v2t, CT_DelaunayVertex *v3t)
{
    _v1 = v1t;
    _v2 = v2t;
    _v3 = v3t;

    _n12 = NULL;
    _n23 = NULL;
    _n31 = NULL;

    _toRemove = false;

    calculateCircle();
}

CT_DelaunayTriangle::~CT_DelaunayTriangle()
{
    _v1 = NULL;
    _v2 = NULL;
    _v3 = NULL;
    _n12 = NULL;
    _n23 = NULL;
    _n31 = NULL;
}

void CT_DelaunayTriangle::init(CT_DelaunayVertex *v1t, CT_DelaunayVertex *v2t, CT_DelaunayVertex *v3t)
{
    _v1 = v1t;
    _v2 = v2t;
    _v3 = v3t;

    _n12 = NULL;
    _n23 = NULL;
    _n31 = NULL;

    _toRemove = false;

    calculateCircle();
}

void CT_DelaunayTriangle::setNeighbor(CT_DelaunayVertex *vt1, CT_DelaunayVertex *vt2, CT_DelaunayTriangle *ngb)
{

    if        (vt1 == _v1) {

        if      (vt2 == _v2) {_n12 = ngb;}
        else if (vt2 == _v3) {_n31 = ngb;}

    } else if (vt1 == _v2) {

        if      (vt2 == _v1) {_n12 = ngb;}
        else if (vt2 == _v3) {_n23 = ngb;}

    } else if (vt1 == _v3) {

        if      (vt2 == _v1) {_n31 = ngb;}
        else if (vt2 == _v2) {_n23 = ngb;}

    }
}

void CT_DelaunayTriangle::calculateCircle()
{

    double dx12 = _v2->x() - _v1->x();
    double dy12 = _v2->y() - _v1->y();
    double dx23 = _v3->x() - _v2->x();
    double dy23 = _v3->y() - _v2->y();
    double x12  = (_v1->x() + _v2->x())/2.0;
    double y12  = (_v1->y() + _v2->y())/2.0;
    double x23  = (_v2->x() + _v3->x())/2.0;
    double y23  = (_v2->y() + _v3->y())/2.0;

    // vertices are aligned :
    // should not happen during adding delaunay triangulation method
    // (if two point never have same position: verified by DelaunayTriangulation.doInsertion ())
    if ((dx12*dy23 - dy12*dx23) == 0)
    {
        qDebug() << "Vertices " << _v1 << ", " << _v2 << ", " << _v3 << " aligned";
        _ccX = NAN;
        _ccY = NAN;
        _r = NAN;
        return;
    }

    // (v1v2) is horizontal
    if (dy12 == 0)
    {
        _ccX = x12;
        _ccY = y23 - (_ccX - x23)*dx23/dy23;

        //  (v1v2) is oblique
    } else {
        if (dy23 == 0)
        {
            _ccX = x23;
            _ccY = y12 - (_ccX - x12)*dx12/dy12;
        } else
        {
            _ccX = (y23 - y12 + dx23*x23/dy23 - dx12*x12/dy12) / (dx23/dy23 - dx12/dy12);
            _ccY = y23 - (_ccX - x23)*dx23/dy23;
        }
    }

    // circle radius calculation
    _r = sqrt (pow(_ccX - _v1->x(),2) + pow(_ccY - _v1->y(),2));
}

bool CT_DelaunayTriangle::circleContains(double xt, double yt)
{
    return ( _r >= (sqrt (pow(_ccX - xt,2) + pow(_ccY - yt,2))) );
}

bool CT_DelaunayTriangle::contains(double xt, double yt) const
{
    CT_DelaunayVertex* max;
    CT_DelaunayVertex* med;
    CT_DelaunayVertex* min;
    double x1, x2;

    // 1) sorting vertices by y coordinate

    if (_v1->y() > _v2->y())
    {
        max = _v1;
        med = _v2;
    } else
    {
        max = _v2;
        med = _v1;
    }

    if (_v3->y() > max->y())
    {
        min = med;
        med = max;
        max = _v3;
    } else if (_v3->y() > med->y())
    {
        min = med;
        med = _v3;
    } else
    {
        min = _v3;
    }

    // test if point is in the good y coordinate range
    if ((yt > max->y()) || (yt < min->y()))
    {
        return false;
    } else
    {

        // 2) caclulation of the point on max-min which split the triangle : (xs, ys)
        double ys = med->y();
        double xs = (ys - max->y()) * (max->x() - min->x()) / (max->y() - min->y()) + max->x();


        // 3) top subTriangle analysis
        if (yt > ys)
        {
            // x-bounds calculations for xt to be in subTriangle, with y=yt=constant
            x1 = (yt - max->y()) * (max->x() - med->x()) / (max->y() - med->y()) + max->x();
            x2 = (yt - max->y()) * (max->x() - xs) / (max->y() - ys) + max->x();

            // 4) bottom subTriangle analysis
        } else if (yt < ys)
        {
            // x-bounds calculations for xt to be in subTriangle, with y=yt=constant
            x1 = (yt - med->y()) * (med->x() - min->x()) / (med->y() - min->y()) + med->x();
            x2 = (yt - ys) * (xs - min->x()) / (ys - min->y()) + xs;

            // 5) yt==ys
        } else
        {
            // x-bounds calculations for xt to be in subTriangle, with y=yt=constant
            x1 = xs;
            x2 = med->x();
        }

        // test if xt is included between x1 and x2 bounds
        if (x1 < x2)
        {
            return ( (xt >= x1) && (xt <= x2) );
        } else
        {
            return ( (xt>=x2) && (xt<=x1) );
        }

    }
}

double CT_DelaunayTriangle::getBaryX()
{
    return (_v1->x() + _v2->x() + _v3->x()) / 3.0;
}

double CT_DelaunayTriangle::getBaryY()
{
    return (_v1->y() + _v2->y() + _v3->y()) / 3.0;
}

CT_DelaunayVertex *CT_DelaunayTriangle::getThirdVertex(CT_DelaunayVertex *vt1, CT_DelaunayVertex *vt2)
{

    if      ((vt1 == _v1) && (vt2 == _v2)) {return _v3;}
    else if ((vt1 == _v2) && (vt2 == _v3)) {return _v1;}
    else if ((vt1 == _v1) && (vt2 == _v3)) {return _v2;}
    else {return NULL;}
}

CT_DelaunayVertex *CT_DelaunayTriangle::getThirdVertex(const QVector<CT_DelaunayVertex *> &vt)
{
    if      ((vt.at(0) == _v1) && (vt.at(1) == _v2)) {return _v3;}
    else if ((vt.at(0) == _v2) && (vt.at(1) == _v3)) {return _v1;}
    else if ((vt.at(0) == _v1) && (vt.at(1) == _v3)) {return _v2;}
    else                                           {return NULL;}
}

CT_DelaunayTriangle* CT_DelaunayTriangle::getNextTriangleTo(double xt, double yt)
{

    double dx, dy;
    double a, b;
    bool n12ok = false;
    bool n23ok = false;
    bool n31ok = false;

    if (_n12 != NULL)
    {
        dx = _v2->x() - _v1->x();
        dy = _v2->y() - _v1->y();

        if (dx == 0) {
            if (((_v3->x() < _v1->x()) && (xt > _v1->x())) || ((_v3->x() > _v1->x()) && (xt < _v1->x())))
            {
                n12ok = true;
            }

        } else  {
            a = (_v1->y() + (dy/dx)*(_v3->x() - _v1->x()));
            b = (_v1->y() + (dy/dx)*(xt - _v1->x()));

            if (((_v3->y() < a) && (yt > b)) || ((_v3->y() > a) && (yt < b)))
            {
                n12ok = true;
            }
        }
    }

    if (_n23 != NULL)
    {
        dx = _v3->x()-_v2->x();
        dy = _v3->y()-_v2->y();

        if (dx == 0)
        {
            if (((_v1->x() < _v2->x()) && (xt > _v2->x())) || ((_v1->x() > _v2->x()) && (xt < _v2->x())))
            {
                n23ok = true;
            }

        } else
        {
            a = (_v2->y() + (dy/dx)*(_v1->x() - _v2->x()));
            b = (_v2->y() + (dy/dx)*(xt - _v2->x()));

            if (((_v1->y() < a) && (yt > b)) || ((_v1->y() > a) && (yt < b)))
            {
                n23ok = true;
            }
        }
    }

    if (_n31 != NULL)
    {
        dx = _v1->x()-_v3->x();
        dy = _v1->y()-_v3->y();

        if (dx == 0) {
            if (((_v2->x() < _v3->x()) && (xt > _v3->x())) || ((_v2->x() > _v3->x()) && (xt < _v3->x())))
            {
                n31ok = true;
            }

        } else
        {
            a = (_v3->y() + (dy/dx)*(_v2->x() - _v3->x()));
            b = (_v3->y() + (dy/dx)*(xt - _v3->x()));

            if (((_v2->y() < a) && (yt > b)) || ((_v2->y() > a) && (yt < b)))
            {
                n31ok = true;
            }
        }
    }


    if ((n12ok)  && (!n23ok) && (!n31ok)) {return _n12;}
    if ((!n12ok) && (n23ok)  && (!n31ok)) {return _n23;}
    if ((!n12ok) && (!n23ok) && (n31ok))  {return _n31;}


    // stochastic process
    // (to avoid cycles, which normally sould not happen in delaunay triangulation)

    if ((n12ok) && (n23ok)) {if ((std::rand() % 2) == 1) {return _n12;} else {return _n23;}}
    if ((n23ok) && (n31ok)) {if ((std::rand() % 2) == 1) {return _n23;} else {return _n31;}}
    if ((n12ok) && (n31ok)) {if ((std::rand() % 2) == 1) {return _n12;} else {return _n31;}}

    // if (n12ok==false) && (n23ok==false) && (n31ok==false) :
    // then point (xt,yt) is in the this triangle
    return this;

}

CT_Polygon2DData *CT_DelaunayTriangle::getShape()
{
    QVector<Eigen::Vector2d *> vertices;

    vertices.append(new Eigen::Vector2d(_v1->x(), _v1->y()));
    vertices.append(new Eigen::Vector2d(_v2->x(), _v2->y()));
    vertices.append(new Eigen::Vector2d(_v3->x(), _v3->y()));

    return new CT_Polygon2DData(vertices, false);
}