poly_utils.cpp

Go to the documentation of this file.

#include <poly_utils.h>
#include <frechet/poly/poly_utils.h>


using namespace frechet;
using namespace poly;

PolygonUtilities::PolygonUtilities()
 : cgp()
{ }

PolygonUtilities::PolygonUtilities(const Curve& p)
 : cgp()
{
    Q_ASSERT(p.isClosed());

    // we index vertices in ascending order; so we can detecte diagonals later
    auto v = p.cbegin();
    //  with QPolygonF, the first vertex is duplicated as the last vertex
    //  with CGAL::Polygon_2 is IS NOT !
    auto vend = p.cend()-1;

    for(int i=0; v!=vend; ++v,++i)
        cgp.push_back(CgalPoint(*v,i));

    Q_ASSERT(assertStrictlySorted(cgp));
    Q_ASSERT(cgp.size()+1 == p.size());
}

bool PolygonUtilities::is_simple() const
{    
    return cgp.is_simple();
}

bool PolygonUtilities::is_convex() const
{
    return cgp.is_convex();
}

bool PolygonUtilities::is_counter_clockwise()
{
    return cgp.is_counterclockwise_oriented();
}

int PolygonUtilities::decompose(DecompositionAlgorithm algo, bool verify)
{
    convex.clear();
    Q_ASSERT(assertStrictlySorted(cgp));

    try {
        PartitionTraits traits;
        switch(algo)
        {
        case APPROX_GREENE:
            //  Greene approximation, O(n log n) (approx. ratio 4, but rarely)
            CGAL::greene_approx_convex_partition_2(
               cgp.vertices_begin(), cgp.vertices_end(),
               std::back_inserter(convex), traits);
            break;
        case APPROX_HERTEL:
            //  4-approximation of Hertel and Mehlhorn
            //  time O(n)
            CGAL::approx_convex_partition_2(
               cgp.vertices_begin(), cgp.vertices_end(),
               std::back_inserter(convex),
               traits);
            break;
        case OPTIMAL_GREENE:
            //  optimal, Greene, O(n^4)
            CGAL::optimal_convex_partition_2(
               cgp.vertices_begin(), cgp.vertices_end(),
               std::back_inserter(convex),
               traits);
            break;
        }

        if (verify) {
            bool ok = CGAL::convex_partition_is_valid_2(
                    cgp.vertices_begin(), cgp.vertices_end(),
                    convex.begin(), convex.end(),
                    traits);
            //Q_ASSERT(ok);
            if (!ok) return -1;
        }
    } catch (std::exception& error) {
        std::cerr << error.what() << std::endl;
        return -1;
    }

    Q_ASSERT(assertStrictlySorted(cgp));
    return convex.size();
}

int PolygonUtilities::partition(Partition& partition) const
{    
    //  assumes decompose has already been called
    Q_ASSERT(! convex.empty());
    partition.clear();

    int n = cgp.size();
    //  copy to partition
    for(const auto& part : convex)
    {
        partition.push_back(Polygon());

        Q_ASSERT(! part.is_empty());
        Q_ASSERT(part.size() >= 3);

        auto v1 = part.vertices_begin();
        auto v2 = v1+1;
        auto vend = part.vertices_end();

        for( ; v2!=vend; ++v1,++v2)
            create1Diagonal(v1->i,v2->i, n, partition);

        v2 = part.vertices_begin();
        create1Diagonal(v1->i,v2->i, n, partition);
    }

    return (int) partition.size();   //  ?
}

bool PolygonUtilities::create1Diagonal(int u, int v, int n, Partition& partition)
{
    Segment seg(u,v);
    if (! seg.is_diagonal(n)) return false;

    partition.back().push_back(u);
    partition.back().push_back(v);
    return true;
}

bool PolygonUtilities::assertStrictlySorted(const CgalPolygon& cgp)
{
    auto p = cgp.vertices_begin();
    auto q = p+1;
    auto end = cgp.vertices_end();

    int j=0;
    for( ; q!=end; ++p,++q,++j) {
        Q_ASSERT(p->i >= 0);
        Q_ASSERT(q->i >= 0);
        Q_ASSERT(p->i==j);
        Q_ASSERT(p->i+1 == q->i);
    }
    return true;
}


int PolygonUtilities::simplify()
{
    Q_ASSERT(!convex.empty());
    CgalPolygonList::iterator i = convex.begin();
    CgalPolygonList::iterator iend = convex.end();

//    std::cout << convex << std::endl;

    while(i != iend) {
        CgalPolygon& poly = *i;
        if (! simplify(poly))
            i = convex.erase(i);
        else
            ++i;
    }
    return convex.size();
//    std::cout << convex << std::endl;
}

bool PolygonUtilities::simplify(CgalPolygon& poly)
{
    if (poly.size() <= 2)
        return false;

    //  contract all "real" edges
    int n = cgp.size();
    CgalPolygon::Vertex_circulator pbegin = poly.vertices_circulator();
    CgalPolygon::Vertex_circulator p = pbegin;

    //  mark for erase
    do {
        if (Segment((p-1)->i, p->i).is_edge(n)
            && Segment(p->i,(p+1)->i).is_edge(n))
            p->erase=true;
    } while(++p != pbegin);

    //  do erase
    auto q = poly.vertices_begin();
    while(q != poly.vertices_end())
        if (q->erase)
            q = poly.erase(q);
        else
            ++q;

    return poly.size() > 2;
}

static int distance(int a, int b, int n)
{
    if (b < a) b += n;
    return b-a;
}

static bool is_edge(int a, int b, int n)
{
    return distance(a,b,n) == 1;
}

static bool is_diagonal(int a, int b, int n)
{
    return distance(a,b,n) > 1;
}

void PolygonUtilities::fanTriangulate(const CgalPolygon& poly, Triangulation& tri) const
{
    auto p0 = poly.vertices_begin();
    auto p1 = p0+1;

    if (poly.size() < 3) {
        //  degenerated poly.
        //  Diagonals are actually outer edges and have no neighbor.
        Q_ASSERT(false);
//        if (poly.size()==2)
//            tri.outerEdge(Segment(p0->i,p1->i));
        return;
    }

    //  Faces erstellen
    auto p2 = p0+2;
    auto pend = poly.vertices_end();

    tri.createVertex((Point)*p0, p0->i);
    tri.createVertex((Point)*p1, p1->i);

    for( ; p2!=pend; p1=p2++) {
        tri.createVertex((Point)*p2, p2->i);
        tri.createFace(p0->i,p1->i,p2->i);
    }
}

//  erstellt ein Fan-Triangulierung mehrerer (konvexer!) Polygone
void PolygonUtilities::fanTriangulate(Triangulation& tri) const
{
    Q_ASSERT(!convex.empty());
    for(const CgalPolygon& poly : convex)
        fanTriangulate(poly,tri);

    tri.complement();
    //  all vertices must be defined
    Q_ASSERT(tri.isComplete());    
    //Q_ASSERT(tri.size()==cgp.size());
    //Q_ASSERT(tri.countFaces() > convex.size());
    //Q_ASSERT(tri.countVertexes()>=0 && tri.countVertexes()<=cgp.size());
#ifdef QT_DEBUG
    for (int i=0; i < cgp.size(); ++i)
    {
        Triangulation::Vertex_handle v = tri[i];
        if (v==Triangulation::Vertex_handle()) continue;
        Q_ASSERT(v->pindex == i);
        Q_ASSERT(v->pindex == cgp[i].i);
        Q_ASSERT(cgp[i].erase==false);
        Q_ASSERT(((Point)*v)==((Point)cgp[i]));
    }
#endif
}

void PolygonUtilities::findReflexVertices(frechet::poly::Polygon& reflex, double tolerance)
{
    CgalPolygon::Vertex_circulator pbegin = cgp.vertices_circulator();
    CgalPolygon::Vertex_circulator p = pbegin;

    do {
        if (is_reflex_vertex(p,tolerance))
            reflex.push_back(p->i);
    } while(++p != pbegin);
}

bool PolygonUtilities::is_reflex_vertex(const CgalPolygon::Vertex_circulator& v, double tolerance)
{
    const CgalPoint& a = *(v-1);
    const CgalPoint& b = *v;
    const CgalPoint& p = *(v+1);

    // TODO code is duplicated from PolygonShortestPaths::halfplaneTest; not necessary
    double A = (p.y()-a.y()) * (b.x()-a.x());
    double B = (p.x()-a.x()) * (b.y()-a.y());
    return A < B+tolerance;
}


CgalPoint::CgalPoint() : super(),i(-1),erase(false)
{ }


CgalPoint::CgalPoint(const QPointF& p, int ai)
: CgalPoint(p.x(),-p.y(),ai)
{ }

CgalPoint::CgalPoint(double x, double y, int ai)
: super(x,y), i(ai), erase(false)
{ }

CgalPoint& CgalPoint::operator= (const QPointF& p)
{
    super::operator=(super(p.x(),-p.y()));
    i = -1;
    return *this;
}

CgalPoint::operator QPointF() const
{
    return QPointF(x(),-y());
}


std::ostream& frechet::poly::operator<< (std::ostream& out, const frechet::poly::CgalPolygon& poly)
{
    out << "[";
    auto pbegin = poly.vertices_begin();
    auto pend = poly.vertices_end();
    for(auto p=pbegin; p != pend; ++p)
    {
        if (p != pbegin) out <<",";
        out << p->i;
    }
    return out << "]";
}

std::ostream& frechet::poly::operator<< (std::ostream& out, const frechet::poly::CgalPolygonList& polys)
{
    out << "{";
    auto pbegin = polys.begin();
    auto pend = polys.end();
    for(auto p=pbegin; p != pend; ++p)
    {
        if (p != pbegin) out << ", ";
        out << *p;
    }
    return out << "}";
}