VoronoiTesselation.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 "VoronoiTesselation.h"

#include "editor/PolygonUtils.h"
#include "util/container/circular_iterator.h"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <functional>
#include <QGraphicsLineItem>
#include <QGraphicsScene>
#include <QGraphicsSceneMouseEvent>
#include <unordered_set>
#include <utility>
#include <vector>

using namespace SimPT_Sim::Container;
using boost::polygon::point_data;
using boost::polygon::voronoi_diagram;

namespace SimPT_Editor {

const QRectF VoronoiTesselation::g_point_rect(-0.5, -0.5, 1.0, 1.0);
const QBrush VoronoiTesselation::g_point_brush(Qt::black);
const double VoronoiTesselation::g_points_precision = 1.0e3;

VoronoiTesselation::VoronoiTesselation(const QPolygonF &boundaryPolygon, QGraphicsItem *parent)
        : QGraphicsPolygonItem(PolygonUtils::Counterclockwise(PolygonUtils::OpenPolygon(boundaryPolygon)), parent)
{
        setFlag(QGraphicsItem::ItemClipsChildrenToShape, true);
}

VoronoiTesselation::~VoronoiTesselation()
{
}

std::list<QPolygonF> VoronoiTesselation::GetCellPolygons()
{
        std::list<QPolygonF> clippedPolygons;
        for (const auto &cellPolygon : m_cell_polygons) {
                clippedPolygons.splice(clippedPolygons.end(), PolygonUtils::ClipPolygon(cellPolygon, polygon()));
        }

        return clippedPolygons;
}

void VoronoiTesselation::mousePressEvent(QGraphicsSceneMouseEvent *event)
{
        QGraphicsPolygonItem::mousePressEvent(event);

        if (event->button() == Qt::RightButton) {
                QGraphicsEllipseItem *point = dynamic_cast<QGraphicsEllipseItem*>(scene()->itemAt(event->scenePos(), QTransform()));
                if (point) {
                        auto it = std::find(m_points.begin(), m_points.end(), point->pos());
                        assert(it != m_points.end() && "m_points should contain the point that has to be removed");
                        m_points.erase(it);

                        delete point;

                        UpdateTesselation();
                }
        }
}

void VoronoiTesselation::mouseDoubleClickEvent(QGraphicsSceneMouseEvent *event)
{
        QGraphicsPolygonItem::mouseDoubleClickEvent(event);

        if (event->button() == Qt::LeftButton) {
                m_points.push_back(event->pos());

                QGraphicsEllipseItem *point = new QGraphicsEllipseItem(g_point_rect, this);
                point->setPos(event->pos());
                point->setBrush(g_point_brush);
                point->setPen(Qt::NoPen);

                UpdateTesselation();
        }
}

void VoronoiTesselation::UpdateTesselation()
{
        std::vector<point_data<long>> points;
        std::transform(m_points.begin(), m_points.end(), std::back_inserter(points), &ToVoronoiPoint);

        voronoi_diagram<double> voronoi;
        construct_voronoi(points.begin(), points.end(), &voronoi);

        m_cell_polygons.clear();

        std::map<const voronoi_diagram<double>::edge_type*, ClippedEdge> clippingCache;

        for (auto &cell : voronoi.cells()) {
                auto edge = cell.incident_edge();
                if (!edge) // No edges
                        continue;

                QPolygonF polygon;
                do {
                        AddEdge(edge, polygon, clippingCache);
                        edge = edge->next();
                } while (edge != cell.incident_edge());
                m_cell_polygons.push_back(polygon);
        }

        RedrawLines();
}

void VoronoiTesselation::RedrawLines()
{
        std::for_each(m_lines.begin(), m_lines.end(), [] (QGraphicsLineItem *item) { delete item; });
        m_lines.clear();

         // Set of line already drawn, so we wont draw two line at the same place (which is ugly, when antialiassed)
        std::hash<double> hash_qreal;
        auto hash_qlinef = [hash_qreal] (const QLineF &line) { return 29 * hash_qreal(line.p1().x()) + 53 * hash_qreal(line.p1().y()) + 29 * hash_qreal(line.p2().x()) + 53 * hash_qreal(line.p2().y()); };
        auto equal_qlinef = [] (const QLineF &line1, const QLineF &line2) { return (line1.p1() == line2.p1() && line1.p2() == line2.p2()) || (line1.p1() == line2.p2() && line1.p2() == line2.p1()); };
        std::unordered_set<QLineF, decltype(hash_qlinef), decltype(equal_qlinef)> drawnLines(10, hash_qlinef, equal_qlinef);

        for (const auto &polygon : m_cell_polygons) {
                std::vector<QPointF> polygonVector = polygon.toStdVector();

                auto cit = make_const_circular(polygonVector.begin(), polygonVector.end(), polygonVector.begin());
                do {
                        QLineF line(*cit, *next(cit));
                        if (drawnLines.count(line) == 0) {
                                drawnLines.insert(line);
                                QGraphicsLineItem *lineItem = new QGraphicsLineItem(line, this);
                                lineItem->setZValue(-1);
                                m_lines.push_back(lineItem);
                        }
                } while (++cit != polygonVector.begin());
        }
}

void VoronoiTesselation::AddEdge(const voronoi_diagram<double>::edge_type *edge, QPolygonF &cellPolygon, std::map<const voronoi_diagram<double>::edge_type*, ClippedEdge> &clippingCache)
{
        ClippedEdge clippedEdge = GetClippedEdge(edge, clippingCache);

        // Check if edge not completely outside of the bounding rect, else this should be handled by a previous or future call to AddEdge
        if (clippedEdge.clippedP1 || clippedEdge.clippedP2 || boundingRect().contains(clippedEdge.edge.p1())) {
                cellPolygon.append(clippedEdge.edge.p1());

                if (clippedEdge.clippedP2) { // If the endpoint was clipped, the voronoi edge leaves the bounding rect, so we find the edge that reenters the bounding rect and append all points of the boundaryRect polygon in between
                        cellPolygon.append(clippedEdge.edge.p2());

                        auto reenterEdge = edge->next();
                        ClippedEdge reenterClipped;
                        do {
                                reenterClipped = GetClippedEdge(reenterEdge, clippingCache);
                                assert((reenterClipped.clippedP1 || reenterEdge != edge) && "Can't find edge reentering the bounding rect");
                                reenterEdge = reenterEdge->next();
                        } while (!reenterClipped.clippedP1);

                        // Add points in between of exiting and reentering of edge
                        std::vector<QPointF> boundingPolygon = PolygonUtils::OpenPolygon(QPolygonF(boundingRect())).toStdVector(); // Qt Documentation says: clockwise order, but this is as drawn on the screen (negative y-axis = up, so lefthanded), so already counterclockwise in a standard right-handed coordinate system (which is good, that's how we need it)
                        assert(!PolygonUtils::IsClockwise(QPolygonF::fromStdVector(boundingPolygon)) && "Qt, get your documentation right!");

                        auto nextBoundaryPoint = std::find(boundingPolygon.begin(), boundingPolygon.end(), clippedEdge.nextBoundaryP2);
                        auto reenterBoundaryPoint = std::find(boundingPolygon.begin(), boundingPolygon.end(), reenterClipped.nextBoundaryP1);
                        assert(nextBoundaryPoint != boundingPolygon.end() && "nextBoundaryPoint not found in boundingRect polygon");
                        assert(reenterBoundaryPoint != boundingPolygon.end() && "reenterBoundaryPoint not found in boundingRect polygon");

                        auto cit = make_const_circular(boundingPolygon.begin(), boundingPolygon.end(), nextBoundaryPoint);
                        while (cit != reenterBoundaryPoint) {
                                cellPolygon.append(*cit);
                                ++cit;
                        }
                }
        }
}

VoronoiTesselation::ClippedEdge VoronoiTesselation::GetClippedEdge(const boost::polygon::voronoi_diagram<double>::edge_type *edge, std::map<const boost::polygon::voronoi_diagram<double>::edge_type*, ClippedEdge> &clippingCache)
{
        if (clippingCache.find(edge) != clippingCache.end()) {
                return clippingCache[edge];
        } else if (clippingCache.find(edge->twin()) != clippingCache.end()) {
                auto reverse = clippingCache[edge->twin()];
                return ClippedEdge{QLineF(reverse.edge.p2(), reverse.edge.p1()), reverse.clippedP2, reverse.nextBoundaryP2, reverse.clippedP1, reverse.nextBoundaryP1};
        } else {
                auto clipped = ClipEdge(edge);
                clippingCache[edge] = clipped;
                return clipped;
        }
}

VoronoiTesselation::ClippedEdge VoronoiTesselation::ClipEdge(const boost::polygon::voronoi_diagram<double>::edge_type *edge)
{
        QPointF c1 = m_points[edge->cell()->source_index()];
        QPointF c2 = m_points[edge->twin()->cell()->source_index()];
        QPointF direction = QLineF(QPointF(), QPointF(c1.y() - c2.y(), c2.x() - c1.x())).unitVector().p2(); // Vector perpendicular to the direction of the line between the centers of the two adjacent cells, in the direction from vertex1 to vertex2
        direction *= 1 + 1e-12; // Wachting out for fp errors (if the line happens to start in a corner and goes in the direction of the other

        double boundingDiagonal = QLineF(boundingRect().topLeft(), boundingRect().bottomRight()).length();

        ClippedEdge clipped{};
        if (edge->vertex0() && edge->vertex1()) { // both vertex0 and vertex
                clipped.edge.setP1(FromVoronoiPoint(*edge->vertex0()));
                clipped.edge.setP2(FromVoronoiPoint(*edge->vertex1()));
        } else if (edge->vertex0()) { // vertex1 is located at infinity
                clipped.edge.setP1(FromVoronoiPoint(*edge->vertex0()));
                clipped.edge.setP2(clipped.edge.p1() + (boundingDiagonal + QLineF(clipped.edge.p1(), boundingRect().topLeft()).length()) * direction);
        } else if (edge->vertex1()) { // vertex0 is located at infinity
                clipped.edge.setP2(FromVoronoiPoint(*edge->vertex1()));
                clipped.edge.setP1(clipped.edge.p2() - (boundingDiagonal + QLineF(clipped.edge.p2(), boundingRect().topLeft()).length()) * direction);
        } else { // vertex0 and vertex1 are both located at infinity
                QPointF center = (c1 + c2) / 2;
                clipped.edge.setP1(center - boundingDiagonal * direction);
                clipped.edge.setP2(center + boundingDiagonal * direction);
        }

        std::vector<std::pair<QPointF, QPointF>> intersections;
        std::vector<QPointF> boundingPolygon = PolygonUtils::OpenPolygon(QPolygonF(boundingRect())).toStdVector(); // Qt Documentation says: clockwise order, but this is as drawn on the screen (negative y-axis = up, so lefthanded), so already counterclockwise in a standard right-handed coordinate system (which is good, that's how we need it)
        assert(!PolygonUtils::IsClockwise(QPolygonF::fromStdVector(boundingPolygon)) && "Qt, get your documentation right!");

        auto cit = make_const_circular(boundingPolygon.begin(), boundingPolygon.end(), boundingPolygon.begin());
        do {
                QPointF p;
                if (QLineF(*cit, *next(cit)).intersect(clipped.edge, &p) == QLineF::BoundedIntersection) {
                        intersections.push_back(std::make_pair(p, *next(cit))); // *cit as cit is the next point in counterclockwise order, as the polygon has its points in clockwise order
                }
        } while (++cit != boundingPolygon.begin());

        if (intersections.size() > 0) {
                std::sort(intersections.begin(), intersections.end(), [clipped] (const std::pair<QPointF, QPointF> &p, const std::pair<QPointF, QPointF> &q) { return QLineF(clipped.edge.p1(), p.first).length() < QLineF(clipped.edge.p1(), q.first).length(); });
                auto itUnique = std::unique(intersections.begin(), intersections.end());
                size_t uniqueIntersections = std::distance(intersections.begin(), itUnique);

                assert(uniqueIntersections <= 2 && "More than two intersections of a line with a rectangle found");
                assert(uniqueIntersections >= static_cast<size_t>(!edge->vertex0() + !edge->vertex1()) && "An edge should at least be clipped at every point at infinity");

                if (uniqueIntersections == 1) {
                        if (!boundingRect().contains(clipped.edge.p1())) { // p2 completely inside of the boudningRect (not at the edges, where it could be intersected with)
                                clipped.clippedP1 = true;
                                clipped.edge.setP1(intersections[0].first);
                                clipped.nextBoundaryP1 = intersections[0].second;
                        } else if (!boundingRect().contains(clipped.edge.p2())) {
                                clipped.clippedP2 = true;
                                clipped.edge.setP2(intersections[0].first);
                                clipped.nextBoundaryP2 = intersections[0].second;
                        }
                } else if (uniqueIntersections == 2) {
                        clipped.clippedP1 = clipped.edge.p1() != intersections[0].first;
                        clipped.clippedP2 = clipped.edge.p2() != intersections[1].first;
                        clipped.edge.setP1(intersections[0].first);
                        clipped.edge.setP2(intersections[1].first);
                        clipped.nextBoundaryP1 = intersections[0].second;
                        clipped.nextBoundaryP2 = intersections[1].second;

                }
        }

        return clipped;
}

point_data<long> VoronoiTesselation::ToVoronoiPoint(const QPointF &point)
{
        return point_data<long>(std::roundl(point.x() * g_points_precision), std::roundl(point.y() * g_points_precision));
}

QPointF VoronoiTesselation::FromVoronoiPoint(const voronoi_diagram<double>::vertex_type &point)
{
        return QPointF(point.x() / g_points_precision, point.y() / g_points_precision);
}

} // namespace