
// ##################################################
// ##    Algoritmus Delaunayovy triangulace        ##
// ##                                              ##
// ##################################################

#pragma once
#ifndef DELAUNAY_TRIANGULATION_H
#define DELAUNAY_TRIANGULATION_H

#include "common.h"

#include <cassert>
#include <cmath>
#include <set>

#include "mathfunc.h"

// ################################

template<typename Site, typename DistanceMeasure>
class delaunay_vertex
{
public:
    typedef vector_2d<DistanceMeasure> point;
private:
    Site* site_;
    point pt_;
public:
    /* konstruktory */
    delaunay_vertex(const delaunay_vertex& v)                         { site_ = v.site_, pt_ = v.pt_;      }
    delaunay_vertex(Site* s)                                          { pt_.x = (double)s->getPosition().x, pt_.y = (double)s->getPosition().y; site_ = s; }
    delaunay_vertex(const point& pt)                                  { pt_ = pt, site_ = NULL;            }
    delaunay_vertex(const DistanceMeasure x, const DistanceMeasure y) { pt_ = point(x,y), site_ = NULL;    }

    /* zjisteni pozice */
    DistanceMeasure getX(void) const     { return pt_.x; };
    DistanceMeasure getY(void) const     { return pt_.y; };
    const point& getPosition(void) const { return pt_;   };
    Site*           getSite(void) const  { return site_; };

    /* porovnavaci operator */
	bool operator<(const delaunay_vertex& v) const
	{
        return (pt_.x == v.pt_.x) ? pt_.y < v.pt_.y : pt_.x < v.pt_.x;
	}

    /* porovnavaci operator */
	bool operator==(const delaunay_vertex& v) const
	{
        return pt_.x == v.pt_.x && pt_.y == v.pt_.y;
	}
	bool operator!=(const delaunay_vertex& v) const
	{
        return !(*this == v);
	}	
};

// ################################

template<typename Site, typename DistanceMeasure>
class delaunay_triangle 
{
public:
    typedef vector_2d<DistanceMeasure> point;
private:
    const delaunay_vertex<Site, DistanceMeasure>* a_;
    const delaunay_vertex<Site, DistanceMeasure>* b_;
    const delaunay_vertex<Site, DistanceMeasure>* c_;
    DistanceMeasure ccRadius_;
    DistanceMeasure ccRadiusSquared_;
    point circumCenter_;
public:
    /* kopirovaci konstruktor */
    delaunay_triangle(const delaunay_triangle& t) 
    {
        a_ = t.a_, b_ = t.b_, c_ = t.c_;
        ccRadius_ = t.ccRadius_;
        ccRadiusSquared_ = t.ccRadiusSquared_;
        circumCenter_ = t.circumCenter_;
    }

    const delaunay_vertex<Site, DistanceMeasure>* getVertex(int i) const
    {
        switch (i)
        {
        case 1: return a_;
        case 2: return b_;
        case 3: return c_;
        default:
            return NULL;
        }
    }

    /* zakladni konstruktor */
    delaunay_triangle(const delaunay_vertex<Site, DistanceMeasure>* a, const delaunay_vertex<Site, DistanceMeasure>* b, 
        const delaunay_vertex<Site, DistanceMeasure>* c)
    {
        a_ = a, b_ = b, c_ = c;
        computeCircumscribedCircle();
    }

    /* porovnavaci operator */
    bool operator<(const delaunay_triangle& t) const
    {
        return circumCenter_.x == t.circumCenter_.x ? circumCenter_.y < t.circumCenter_.y :
            circumCenter_.x < t.circumCenter_.x;
    }

    /* zjisti, zda je delaunay_vertex v kruznici opsane tohoto trojuhelniku */
    bool isInCircumCircle(const delaunay_vertex<Site, DistanceMeasure>& v) const
    {
        return (v.getPosition() - circumCenter_).squaredNorm() <= ccRadiusSquared_;
    }

    /* zjisti, zda je dany vertex napravo od kruznice opsane */
    bool isVertexRightOfCircumCircle(const delaunay_vertex<Site, DistanceMeasure>& v) const
    {
        return (circumCenter_.x + ccRadius_) < v.getX();
    }

    void computeCircumscribedCircle(void)
    {
        const DistanceMeasure ax = a_->getX();
        const DistanceMeasure ay = a_->getY();
        const DistanceMeasure bx = b_->getX();
        const DistanceMeasure by = b_->getY();
        const DistanceMeasure cx = c_->getX();
        const DistanceMeasure cy = c_->getY();

        const DistanceMeasure A = bx - ax;
        const DistanceMeasure B = by - ay;
        const DistanceMeasure C = cx - ax;
        const DistanceMeasure D = cy - ay;
        const DistanceMeasure E = A * (bx + ax) + B * (by + ay);
        const DistanceMeasure F = C * (cx + ax) + D * (cy + ay);
        const DistanceMeasure G = 2 * ( A * (cy - by) - B * (cx - bx) );

        /* pokud G je nulove, jsou to kolinearni body, to by tady ale nastat nemelo */
        assert(G != 0);

        circumCenter_.x = ( D * E - B * F ) / G;
        circumCenter_.y = ( A * F - C * E ) / G;

        /* spocitame vektor polomeru - libovolny, a z nej polomer */
        ccRadiusSquared_         = (a_->getPosition() - circumCenter_).squaredNorm(); 
        ccRadius_                = sqrt(ccRadiusSquared_);
        ccRadiusSquared_        *= 1.000001L; // trochu to zvetsime kvuli zaokrouhlovacim chybam
    }                            
};

// ################################

template<typename Site, typename DistanceMeasure>
class delaunay_edge
{
private:
    const delaunay_vertex<Site, DistanceMeasure>* v1_;
    const delaunay_vertex<Site, DistanceMeasure>* v2_;
public:
    /* konstruktory */
    delaunay_edge(const delaunay_edge& edge)   { v1_ = edge.v1_, v2_ = edge.v2_; }
    delaunay_edge(const delaunay_vertex<Site, DistanceMeasure>* v1, 
        const delaunay_vertex<Site, DistanceMeasure>* v2) : v1_(v1), v2_(v2) {}

    /* porovnavaci operator */
    bool operator<(const delaunay_edge& e) const
    {
        return *v1_ == *e.v1_ ? *v2_ < *e.v2_ : *v1_ < *e.v1_;
    }

    const delaunay_vertex<Site, DistanceMeasure>* firstVertex(void)  { return v1_; }
    const delaunay_vertex<Site, DistanceMeasure>* secondVertex(void) { return v2_; }
};

// ################################

template<typename Site, typename DistanceMeasure>
class delaunay_algorithm
{
public:
    typedef std :: set<delaunay_edge<Site, DistanceMeasure>>          delaunay_edge_set;
    typedef std :: multiset<delaunay_triangle<Site, DistanceMeasure>> delaunay_triangle_set;
    typedef std :: set<delaunay_vertex<Site, DistanceMeasure>>        delaunay_vertex_set;

private:
    void _addEdge(const delaunay_vertex<Site, DistanceMeasure> * v1, 
        const delaunay_vertex<Site, DistanceMeasure> * v2, delaunay_edge_set& es)
    {
        const delaunay_vertex<Site, DistanceMeasure>* pv1 = *v1 < *v2 ? v1 : v2;
        const delaunay_vertex<Site, DistanceMeasure>* pv2 = *v1 < *v2 ? v2 : v1;
        es.insert(delaunay_edge<Site, DistanceMeasure>(pv1, pv2));
    }

    // ################################

    /* funktor, ktery zjisti, zda trojuhelnik obsahuje nektery z bodu a,b,c */
    struct triangleContainsVertex
    {
    private:
        const delaunay_vertex<Site, DistanceMeasure> *a_, *b_, *c_;
    public:
        triangleContainsVertex(const delaunay_vertex<Site, DistanceMeasure>* a,
            const delaunay_vertex<Site, DistanceMeasure>* b, const delaunay_vertex<Site, DistanceMeasure>* c)
        {
            a_ = a, b_ = b, c_ = c;
        }

        bool operator()(const delaunay_triangle<Site, DistanceMeasure>& dt) const
        {
            for (int i = 1; i <= 3; i++)
            {
                const delaunay_vertex<Site, DistanceMeasure>* v = dt.getVertex(i);
                if (a_ == v || b_ == v || c_ == v) return true;
            }
            return false;
        }
    };

    // ################################

    struct triangleFinished
    {
    private:
        delaunay_triangle_set& otp_;
        delaunay_vertex<Site, DistanceMeasure> *a_, *b_, *c_, *currentVertex_;
    public:
        triangleFinished(delaunay_vertex<Site, DistanceMeasure> *a, 
            delaunay_vertex<Site, DistanceMeasure> *b, delaunay_vertex<Site, DistanceMeasure> *c, 
            delaunay_vertex<Site, DistanceMeasure> *currentVertex, delaunay_triangle_set& otp) : otp_(otp)
        {
            a_ = a, b_ = b, c_ = c, currentVertex_ = currentVertex;
        }

	    bool operator()(const delaunay_triangle<Site, DistanceMeasure>& dt) 
        {
            /* je kruh nalevo od aktualniho bodu? */
            bool circleIsLeft = dt.isVertexRightOfCircumCircle(*currentVertex_);

            /* pokud ano, vyrad dany trojuhelnik... a dej ho do zpracovanych, pokud to neni supertrojuhelnik */
            if (circleIsLeft)
            {
                triangleContainsVertex tcv(a_, b_, c_);
                if (!tcv(dt)) otp_.insert(dt);
            }

            return circleIsLeft;
        }
    };

    // ################################
    
    /* vytvari nove trojuhelniky */
    struct vertexIsInCircumCircle
    {
    private:
        delaunay_vertex<Site, DistanceMeasure> * vertex_;
        delaunay_edge_set& edges_;

        void addEdge_(const delaunay_vertex<Site, DistanceMeasure> * v1, 
            const delaunay_vertex<Site, DistanceMeasure> * v2)
        {
            const delaunay_vertex<Site, DistanceMeasure>* pv1 = *v1 < *v2 ? v1 : v2;
            const delaunay_vertex<Site, DistanceMeasure>* pv2 = *v1 < *v2 ? v2 : v1;
            delaunay_edge<Site, DistanceMeasure> e(pv1, pv2);

            delaunay_edge_set :: iterator item = edges_.find(e);
            if (item == edges_.end()) edges_.insert(e); else edges_.erase(item);
        }
    public:
        vertexIsInCircumCircle( delaunay_vertex<Site, DistanceMeasure> * vertex, delaunay_edge_set& es) : edges_(es)
        {
            vertex_ = vertex;
        }

	    bool operator()(const delaunay_triangle<Site, DistanceMeasure>& dt) 
	    {
		    bool vertexIsInCC = dt.isInCircumCircle(*vertex_);

		    if (vertexIsInCC)
		    {
			    addEdge_(dt.getVertex(1), dt.getVertex(2));
			    addEdge_(dt.getVertex(2), dt.getVertex(3));
			    addEdge_(dt.getVertex(3), dt.getVertex(1));
		    }
		    return vertexIsInCC;
	    }
    };

    // ################################

public:
    void triangulation(delaunay_vertex_set& vertices, delaunay_triangle_set& outputTriangles,
        delaunay_edge_set& outputEdges)
    {
        /* pozn. predpokladame, ze vertices.size() > 4 */
        if (vertices.size() < 5) 
            return;

        /* udelame bouding box, ktery obsahuje (0,0) */
        DistanceMeasure xMin = vertices.begin()->getX(), yMin = vertices.begin()->getY();
        DistanceMeasure xMax = xMin, yMax = yMin;
        for (delaunay_vertex_set :: const_iterator iter = vertices.begin(); iter != vertices.end(); iter++)
        {
            xMax = xMax < iter->getX() ? iter->getX() : xMax;
            xMin = xMin > iter->getX() ? iter->getX() : xMin;
            yMax = yMax < iter->getY() ? iter->getY() : yMax;
            yMin = yMin > iter->getY() ? iter->getY() : yMin;
        }

        DistanceMeasure  dx = (xMax - xMin);
        DistanceMeasure  dy = (yMax - yMin);
        DistanceMeasure ddx = dx * 5.0L;
        DistanceMeasure ddy = dy * 5.0L;

        /* Zvetsime bounding box o deset procent */
        xMax += ddx; xMin -= ddx; dx  = xMax - xMin;
        yMax += ddy; yMin -= ddy; dy  = yMax - yMin;

        /* sqrt(3)/3 */
        const DistanceMeasure sqrt3d3 = 0.57735026918962576450914878050196L;
        const DistanceMeasure sqrt3d2 = 0.86602540378443864676372317075294L;

        /* udelame obri trojuhelnik, ktery bude vse obsahovat (ze zacatku) */
        delaunay_vertex<Site, DistanceMeasure> verticesBT_a( xMin - dy * sqrt3d3, yMin);
        delaunay_vertex<Site, DistanceMeasure> verticesBT_b( xMax + dy * sqrt3d3, yMin);
        delaunay_vertex<Site, DistanceMeasure> verticesBT_c((xMin + xMax) / 2.0L, yMax + dx * sqrt3d2);

        delaunay_triangle_set workingSet;
        workingSet.insert(delaunay_triangle<Site, DistanceMeasure>(&verticesBT_a, &verticesBT_b, &verticesBT_c));

        /* prochazej postupne vsechny vrcholy */
        for (delaunay_vertex_set :: iterator iter = vertices.begin(); iter != vertices.end(); iter++)
        {
            delaunay_vertex<Site, DistanceMeasure> * currentV = &(*iter);

            /* smaz trojuhelniky, kterych se tento bod uz netyka */
            delaunay_triangle_set :: iterator last_item = std :: remove_if(workingSet.begin(), workingSet.end(), 
                triangleFinished(&verticesBT_a, &verticesBT_b, &verticesBT_c, currentV, outputTriangles));

            delaunay_edge_set edges;

            /* zjisti, ve kterych kruzich opsanych lezi nas novy vrchol a rozbij dane trojuhelniky,
            tj. smaz je ze seznamu trojuhelniku */
            last_item = remove_if(workingSet.begin(), last_item, vertexIsInCircumCircle(currentV, edges));
            workingSet.erase(last_item, workingSet.end());

            /* ty smazane trojuhelniky nahrad spoustou novych trojuhelniku, se stredem v aktualnim vrcholu */
            for (delaunay_edge_set :: iterator eit = edges.begin(); eit != edges.end(); eit++)
            {
                workingSet.insert(delaunay_triangle<Site, DistanceMeasure>(eit->firstVertex(), 
                    eit->secondVertex(), currentV));
            }	
        }

        /* prekopiruj vysledne trojuhelniky */
        remove_copy_if(workingSet.begin(), workingSet.end(), inserter(outputTriangles, outputTriangles.begin()), triangleContainsVertex(&verticesBT_a, &verticesBT_b, &verticesBT_c));

        for (delaunay_triangle_set :: iterator iter = outputTriangles.begin(); iter != outputTriangles.end(); iter++)
        {
            _addEdge(iter->getVertex(1), iter->getVertex(2), outputEdges);
            _addEdge(iter->getVertex(2), iter->getVertex(3), outputEdges);
            _addEdge(iter->getVertex(3), iter->getVertex(1), outputEdges);
        }
    }
};

// ################################

#endif