Graph2DTool.cxx

Go to the documentation of this file.

#include "NAGASH/Graph2DTool.h"
#include "Eigen/Dense"
 
using namespace NAGASH;
using namespace std;
using namespace Eigen;
 
Circle Triangle::CircumscribedCircle()
{
    Matrix3d Bottom;
    Matrix3d XTop;
    Matrix3d YTop;
 
    for (int i = 0; i < 3; i++)
    {
        Bottom(i, 0) = P[i].X();
        Bottom(i, 1) = P[i].Y();
        Bottom(i, 2) = 1;
 
        XTop(i, 0) = pow(P[i].X(), 2) + pow(P[i].Y(), 2);
        XTop(i, 1) = P[i].Y();
        XTop(i, 2) = 1;
 
        YTop(i, 0) = P[i].X();
        YTop(i, 1) = pow(P[i].X(), 2) + pow(P[i].Y(), 2);
        YTop(i, 2) = 1;
    }
 
    Point cent(XTop.determinant() / Bottom.determinant() / 2,
               YTop.determinant() / Bottom.determinant() / 2);
 
    double r = cent.Distance(P[0]);
 
    return Circle(cent, r);
}
 
bool Point::IsInside(Circle c)
{
    return (c.Radius() > Distance(c.Center()));
}
 
bool Point::IsInside(Triangle t)
{
    Point PA(t.Vertex(0).X() - x, t.Vertex(0).Y() - y);
    Point PB(t.Vertex(1).X() - x, t.Vertex(1).Y() - y);
    Point PC(t.Vertex(2).X() - x, t.Vertex(2).Y() - y);
 
    double t1 = PA.X() * PB.Y() - PA.Y() * PB.X();
    double t2 = PB.X() * PC.Y() - PB.Y() * PC.X();
    double t3 = PC.X() * PA.Y() - PC.Y() * PA.X();
 
    return ((t1 * t2 >= 0) && (t1 * t3 >= 0));
}
 
Graph2DTool::Graph2DTool(std::shared_ptr<MSGTool> msg, double prec) : Tool(msg)
{
    PRECISION = prec;
}
 
std::vector<Triangle> Graph2DTool::DelaunayTriangulation(std::vector<Point> p_vec)
{
    vector<Triangle> tri_vec;
    vector<Triangle> temp_tri_vec;
 
    if (p_vec.size() < 3)
    {
        MSGUser()->MSG_ERROR("Need at least 3 points for triangulation!");
        return tri_vec;
    }
 
    vector<int> index_vec;
 
    for (int i = 0; i < p_vec.size(); i++)
        index_vec.push_back(i);
 
    for (int i = 0; i < index_vec.size() - 1; i++)
    {
        for (int j = i + 1; j < index_vec.size(); j++)
        {
            if (p_vec[index_vec[i]].X() > p_vec[index_vec[j]].X())
            {
                int temp = index_vec[i];
                index_vec[i] = index_vec[j];
                index_vec[j] = temp;
            }
        }
    }
 
    double XMax = p_vec[0].X();
    double XMin = p_vec[0].X();
    double YMax = p_vec[0].Y();
    double YMin = p_vec[0].Y();
 
    for (int i = 1; i < p_vec.size(); i++)
    {
        if (p_vec[i].X() > XMax)
            XMax = p_vec[i].X();
        if (p_vec[i].X() < XMin)
            XMin = p_vec[i].X();
        if (p_vec[i].Y() > YMax)
            YMax = p_vec[i].Y();
        if (p_vec[i].Y() < YMin)
            YMin = p_vec[i].Y();
    }
 
    double DX = XMax - XMin;
    double DY = YMax - YMin;
 
    Triangle super(Point((XMax + XMin) / 2, YMax + 2 * DY),
                   Point((XMax + XMin) / 2 + 2.5 * DX, YMin - 2 * DY),
                   Point((XMax + XMin) / 2 - 2.5 * DX, YMin - 2 * DY));
 
    // Check
    for (int i = 0; i < p_vec.size(); i++)
    {
        if (!p_vec[i].IsInside(super))
            MSGUser()->MSG_ERROR("Point ", i, " is outside the super triangle!");
    }
 
    temp_tri_vec.push_back(super);
 
    for (int i = 0; i < index_vec.size(); i++)
    {
        vector<Edge> edg_vec;
        Point P = p_vec[index_vec[i]];
        for (int j = 0; j < temp_tri_vec.size(); j++)
        {
            Circle outcir = temp_tri_vec[j].CircumscribedCircle();
 
            if (P.X() > (outcir.Center().X() + outcir.Radius()))
            {
                tri_vec.push_back(temp_tri_vec[j]);
                temp_tri_vec.erase(temp_tri_vec.begin() + j);
                j--;
                continue;
            }
            else if (P.IsInside(outcir))
            {
                edg_vec.push_back(Edge(temp_tri_vec[j].Vertex(0), temp_tri_vec[j].Vertex(1)));
                edg_vec.push_back(Edge(temp_tri_vec[j].Vertex(1), temp_tri_vec[j].Vertex(2)));
                edg_vec.push_back(Edge(temp_tri_vec[j].Vertex(0), temp_tri_vec[j].Vertex(2)));
 
                temp_tri_vec.erase(temp_tri_vec.begin() + j);
                j--;
                continue;
            }
            else
            {
                continue;
            }
        }
 
        if (edg_vec.size() > 1)
        {
            for (int j = 0; j < edg_vec.size() - 1; j++)
            {
                for (int k = j + 1; k < edg_vec.size(); k++)
                {
                    if ((edg_vec[j].Start().Distance(edg_vec[k].Start()) < PRECISION) &&
                        (edg_vec[j].End().Distance(edg_vec[k].End()) < PRECISION))
                    {
                        edg_vec.erase(edg_vec.begin() + k);
                        k--;
                        continue;
                    }
                    if ((edg_vec[j].Start().Distance(edg_vec[k].End()) < PRECISION) &&
                        (edg_vec[j].End().Distance(edg_vec[k].Start()) < PRECISION))
                    {
                        edg_vec.erase(edg_vec.begin() + k);
                        k--;
                        continue;
                    }
                }
            }
        }
 
        if (edg_vec.size() > 0)
        {
            for (int j = 0; j < edg_vec.size(); j++)
                temp_tri_vec.push_back(Triangle(edg_vec[j], P));
        }
    }
 
    tri_vec.insert(tri_vec.end(), temp_tri_vec.begin(), temp_tri_vec.end());
 
    for (int i = 0; i < tri_vec.size(); i++)
    {
        bool toremove = false;
        for (int j = 0; j < 3; j++)
        {
            for (int k = 0; k < 3; k++)
            {
                if (tri_vec[i].Vertex(j).Distance(super.Vertex(k)) < PRECISION)
                    toremove = true;
            }
        }
        for (int j = 0; j < p_vec.size(); j++)
        {
            bool isOnEdge = false;
            for (int k = 0; k < 3; k++)
            {
                if (p_vec[j].Distance(tri_vec[i].Vertex(k)) < PRECISION)
                    isOnEdge = true;
            }
            if (isOnEdge)
                continue;
            if (p_vec[j].IsInside(tri_vec[i].CircumscribedCircle()))
                toremove = true;
        }
        if (toremove)
        {
            tri_vec.erase(tri_vec.begin() + i);
            i--;
        }
    }
 
    return tri_vec;
}