PolygonUtils.cpp

Go to the documentation of this file.

/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
#include "PolygonUtils.h"

#include "util/container/circular_iterator.h"

#include <QPainterPath>
#include <QVector>
#include <cassert>
#include <cmath>
#include <list>
#include <map>

using namespace SimPT_Sim::Container;

namespace SimPT_Editor {

bool PolygonUtils::IsClockwise(const QPolygonF& polygon)
{
        //The algorithm used to determine this has been adapted from:
        //      [http://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order]

        assert(!polygon.isClosed() && "The polygon isn't open.");
        assert(polygon.size() >= 3 && "The polygon isn't valid.");
        assert(IsSimplePolygon(polygon) && "The polygon isn't simple.");

        double double_area = 0;
        auto polygon_std = polygon.toStdVector();
        auto cit = make_const_circular(polygon_std);
        do {
                double_area += (next(cit)->x() - cit->x()) * (next(cit)->y() + cit->y());
        } while (++cit != polygon_std.begin());

        return (double_area > 0);
}

bool PolygonUtils::IsSimplePolygon(const QPolygonF& polygon)
{
        assert(!polygon.isClosed() && "The polygon isn't open.");
        assert(polygon.size() >= 3 && "The polygon isn't valid.");

        auto polygon_std = polygon.toStdVector();
        auto cit = make_const_circular(polygon_std);

        std::list<QLineF> edges;
        do {
                edges.push_back(QLineF((*cit), *next(cit)));
        } while (++cit != polygon_std.begin());

        for (auto e1 = edges.begin(); e1 != edges.end(); e1++) {
                for (auto e2 = std::next(e1); e2 != edges.end(); e2++) {
                        QPointF intersection;
                        if (e1->intersect(*e2, &intersection) == QLineF::BoundedIntersection
                                && intersection != e1->p1()
                                && intersection != e1->p2()
                                && intersection != e2->p1()
                                && intersection != e2->p2()) {
                                return false;
                        }
                }
        }

        return true;
}

QPolygonF PolygonUtils::OpenPolygon(const QPolygonF &polygon)
{
        assert(polygon.size() >= 3 && "The polygon isn't valid.");

        if (polygon.isClosed()) {
                QPolygonF opened = polygon;
                opened.pop_back();
                return opened;
        } else {
                return polygon;
        }
}

QPolygonF PolygonUtils::Counterclockwise(const QPolygonF& polygon)
{
        assert(!polygon.isClosed() && "The polygon isn't open.");
        assert(polygon.size() >= 3 && "The polygon isn't valid.");
        assert(IsSimplePolygon(polygon) && "The polygon isn't simple.");

        if (PolygonUtils::IsClockwise(polygon)) {
                QPolygonF reversed;
                std::copy(polygon.begin(), polygon.end(), std::front_inserter(reversed));
                return reversed;
        } else {
                return polygon;
        }
}

double PolygonUtils::CalculateArea(const QPolygonF& polygon)
{
        assert(!polygon.isClosed() && "The polygon isn't open.");
        assert(polygon.size() >= 3 && "The polygon isn't valid.");
        assert(IsSimplePolygon(polygon) && "The polygon isn't simple.");

        double area = 0;

        auto polygonStd = polygon.toStdVector();
        auto cit = make_const_circular(polygonStd);
        do {
                area += (*cit).x() * (*next(cit)).y() - (*next(cit)).x() * (*cit).y();
        } while (++cit != polygonStd.begin());

        area /= 2;

        return fabs(area);
}

std::list<QPolygonF> PolygonUtils::ClipPolygon(const QPolygonF& polygon1, const QPolygonF& polygon2)
{
        assert(!polygon1.isClosed() && !polygon2.isClosed() && "The polygon aren't open.");
        assert(polygon1.size() >= 3 && polygon2.size() >= 3 && "The polygon isn't valid.");
        assert(IsSimplePolygon(polygon1) && IsSimplePolygon(polygon2) && "The polygon isn't simple.");



        QPainterPath path1;
        path1.addPolygon(polygon1);
        QPainterPath path2;
        path2.addPolygon(polygon2);

        std::list<QPolygonF> subpathPolygons = path1.intersected(path2).toSubpathPolygons().toStdList();
        std::transform(subpathPolygons.begin(), subpathPolygons.end(), subpathPolygons.begin(), &OpenPolygon);



        // Lambda to check whether to points are equal with a small perturbation.
        auto is_perturbed = [](const QPointF& p1, const QPointF& p2) {
                return fabs(p1.x() - p2.x()) < g_accuracy && fabs(p1.y() - p2.y()) < g_accuracy;
        };

        // Lambda to return a perturbed line.
        auto perturb_line = [](QLineF line) {
                line.setLength(line.length() + g_accuracy);
                QLineF reversed(line.p2(), line.p1());
                reversed.setLength(reversed.length() + g_accuracy);
                return QLineF(reversed.p2(), reversed.p1());
        };

        // Lambda to check whether a given point is located on the given line.
        auto point_on_line = [perturb_line](const QPointF& p, const QLineF& l) {
                QLineF normal = l.normalVector();
                QLineF perpendicular(p, normal.p2() - normal.p1() + p);

                QPointF intersection;
                if (perturb_line(perpendicular).intersect(perturb_line(l), &intersection) == QLineF::BoundedIntersection) {
                        return QLineF(intersection, p).length() < g_accuracy;
                }
                else {
                        return false;
                }
        };



        //Add omitted (during intersection) points to the subpathPolygons.
        for (auto polygon_it = subpathPolygons.begin(); polygon_it != subpathPolygons.end(); polygon_it++) {
                std::list<QPointF> polygonStd = polygon_it->toList().toStdList();

                auto c_prev = [&polygonStd](const decltype(polygonStd.begin())& it) {
                        return (it == polygonStd.begin() ? std::prev(polygonStd.end()) : std::prev(it));
                };

                //Can't use circular iterator for insertion at the start (probably due to the inheritance of prev).
                for (auto cit = polygonStd.begin(); cit != polygonStd.end(); cit++) {
                        QLineF line(*c_prev(cit), *cit);
                        for (const QPointF& p : polygon1) {
                                if (point_on_line(p, line) && !is_perturbed(p, *c_prev(cit)) && !is_perturbed(p, *cit)) {
                                        polygonStd.insert(cit, p);
                                        cit--;
                                        line = QLineF(*c_prev(cit), *cit);
                                }
                        }
                        for (const QPointF& p : polygon2) {
                                if (point_on_line(p, line) && !is_perturbed(p, *c_prev(cit)) && !is_perturbed(p, *cit)) {
                                        polygonStd.insert(cit, p);
                                        cit--;
                                        line = QLineF(*c_prev(cit), *cit);
                                }
                        }
                }
                *polygon_it = QPolygonF::fromList(QList<QPointF>::fromStdList(polygonStd));
        }



        //Check for floating-point errors.
        for (QPolygonF& polygon : subpathPolygons) {
                for (int i = 0; i < polygon.size(); i++) {
                        for (const QPointF& p2 : polygon1) {
                                if (is_perturbed(polygon[i], p2)) {
                                        polygon.replace(i, p2);
                                }
                        }
                        for (const QPointF& p2 : polygon2) {
                                if (is_perturbed(polygon[i], p2)) {
                                        polygon.replace(i, p2);
                                }
                        }
                }
        }

        return subpathPolygons;
}

std::list<QPolygonF> PolygonUtils::SlicePolygon(const QPolygonF& polygon, QLineF cut)
{
        assert(!polygon.isClosed() && "The polygon isn't open.");
        assert(polygon.size() >= 3 && "The polygon isn't valid.");
        assert(IsSimplePolygon(polygon) && "The polygon isn't simple.");
        assert((!polygon.containsPoint(cut.p1(), Qt::OddEvenFill) || polygon.contains(cut.p1()))
                && (!polygon.containsPoint(cut.p2(), Qt::OddEvenFill) || polygon.contains(cut.p2()))
                && "One of the points of the cut is located inside the polygon.");
        assert(cut.p1() != cut.p2() && "The cut isn't a line.");


        QPolygonF initialPolygon;
        std::list<QPointF> intersections;
        auto edges = polygon.toStdVector();


        // Lambda to check whether to points are equal with a small perturbation.
        auto is_perturbed = [](const QPointF& p1, const QPointF& p2) {
                return fabs(p1.x() - p2.x()) < g_accuracy && fabs(p1.y() - p2.y()) < g_accuracy;
        };

        // Lambda to return a perturbed line.
        auto perturb_line = [](QLineF line) {
                line.setLength(line.length() + g_accuracy);
                QLineF reversed(line.p2(), line.p1());
                reversed.setLength(reversed.length() + g_accuracy);
                return QLineF(reversed.p2(), reversed.p1());
        };

        // Lambda to check whether a given point is located on the given line.
        auto point_on_line = [perturb_line](const QPointF& p, const QLineF& l) {
                QLineF normal = l.normalVector();
                QLineF perpendicular(p, normal.p2() - normal.p1() + p);

                QPointF intersection;
                if (perturb_line(perpendicular).intersect(perturb_line(l), &intersection) == QLineF::BoundedIntersection) {
                        return QLineF(intersection, p).length() < g_accuracy;
                }
                else {
                        return false;
                }
        };

        // Lambda to check whether two given lines fall on top of each other (with one of the lines completely overlapping the other).
        auto line_on_line = [point_on_line](const QLineF l1, const QLineF l2) {
                return ((point_on_line(l1.p1(), l2) && point_on_line(l1.p2(), l2)) || (point_on_line(l2.p1(), l1) && point_on_line(l2.p2(), l1)));
        };

        // Find the intersections with the cut.
        auto cit = make_const_circular(edges);
        do {
                QPointF intersection;
                QLineF edge(*cit, *next(cit));
                bool onLine = line_on_line(cut, edge);

                initialPolygon.append(*cit);

                if (perturb_line(cut).intersect(perturb_line(edge), &intersection) == QLineF::BoundedIntersection
                        && !onLine) {
                        //The last check is necessary, because the intersection points could be *prev(cit) and *next(cit), so *cit itself won't be.


                        if (is_perturbed(*cit, intersection)) {
                                intersection = *cit;
                        }
                        else if (is_perturbed(*next(cit), intersection)) {
                                intersection = *next(cit);
                        }
                        else {
                                initialPolygon.append(intersection);
                        }

                        if (intersection != *next(cit)) {
                                // Only consider *cit as an intersection, as *next(cit) will be found in the next iteration.
                                intersections.push_back(intersection);
                        }
                }
                else if (onLine) {
                        //Necessary for polygon partition check in the iteration over the intersections.
                        intersections.push_back(*cit);
                }
        } while (++cit != edges.begin());

        // If there are no intersections (or one), return the original polygon.
        if (intersections.size() <= 1) {
                return std::list<QPolygonF>({polygon});
        }

        // Sort the intersection points by their distance to the first point in the cut.
        auto cmp = [cut](const QPointF& p1, const QPointF& p2) {
                return pow(cut.p1().x() - p1.x(), 2) + pow(cut.p1().y() - p1.y(), 2)
                        < pow(cut.p1().x() - p2.x(), 2) + pow(cut.p1().y() - p2.y(), 2);
        };
        intersections.sort(cmp);


        // Create the subpolygons.
        std::list<QPolygonF> subpolygons({initialPolygon});
        for (auto it = intersections.begin(); it != std::prev(intersections.end()); it++) {
                QLineF line(*it, *std::next(it));

                //Check whether this part of the cut is part of the polygon.
                //In this case, do nothing with it, as it's not really separating two polygons.
                auto firstPoint = std::find(edges.begin(), edges.end(), line.p1());
                if (firstPoint != edges.end()) {
                        std::rotate(edges.begin(), firstPoint, edges.end());
                        if (*std::next(edges.begin()) == line.p2() || *std::prev(edges.end()) == line.p2()) {
                                continue;
                        }
                }

                //If the center of this part of the cut falls outside of the polygon, do nothing,
                //as the cut goes outside of the polygon.
                QPointF center = (line.p1() + line.p2()) / 2;
                if (!polygon.containsPoint(center, Qt::OddEvenFill) && !polygon.contains(center)) {
                        continue;
                }


                auto currentPolygon = std::find_if(subpolygons.begin(), subpolygons.end(), [line](const QPolygonF& p){
                                auto first = std::find(p.begin(), p.end(), line.p1());
                                auto second = std::find(p.begin(), p.end(), line.p2());
                                return first != p.end() && second != p.end();
                        });


                QPolygonF polygon1;
                QPolygonF polygon2;
                QPolygonF* currentAppending = &polygon1;
                for (const QPointF& p : *currentPolygon) {
                        currentAppending->append(p);

                        if (p == line.p1() || p == line.p2()) {
                                currentAppending = (currentAppending == &polygon1 ? &polygon2 : &polygon1);
                                currentAppending->append(p);
                        }
                }
                subpolygons.erase(currentPolygon);
                subpolygons.push_back(polygon1);
                subpolygons.push_back(polygon2);
        }


        // If there are no (useful) intersections, return the original polygon.
        if (subpolygons.size() == 0) {
                return std::list<QPolygonF>({polygon});
        }
        else {
                return subpolygons;
        }
}

int PolygonUtils::Turn(const QLineF& line1, const QLineF& line2)
{
        QPointF u = line1.p2() - line1.p1();
        QPointF v = line2.p2() - line2.p1();
        double z = u.x() * v.y() - u.y() * v.x();
        return (z > 0) - (z < 0); // Sign of z
}

} // namespace