VPTissue Reference Manual
VoronoiTesselation.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2011-2016 Universiteit Antwerpen
3  *
4  * Licensed under the EUPL, Version 1.1 or as soon they will be approved by
5  * the European Commission - subsequent versions of the EUPL (the "Licence");
6  * You may not use this work except in compliance with the Licence.
7  * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the Licence is distributed on an "AS IS" basis,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the Licence for the specific language governing
13  * permissions and limitations under the Licence.
14  */
20 #include "VoronoiTesselation.h"
21 
22 #include "editor/PolygonUtils.h"
24 
25 #include <algorithm>
26 #include <cassert>
27 #include <cmath>
28 #include <functional>
29 #include <QGraphicsLineItem>
30 #include <QGraphicsScene>
31 #include <QGraphicsSceneMouseEvent>
32 #include <unordered_set>
33 #include <utility>
34 #include <vector>
35 
36 using namespace SimPT_Sim::Container;
37 using boost::polygon::point_data;
38 using boost::polygon::voronoi_diagram;
39 
40 namespace SimPT_Editor {
41 
42 const QRectF VoronoiTesselation::g_point_rect(-0.5, -0.5, 1.0, 1.0);
43 const QBrush VoronoiTesselation::g_point_brush(Qt::black);
44 const double VoronoiTesselation::g_points_precision = 1.0e3;
45 
46 VoronoiTesselation::VoronoiTesselation(const QPolygonF &boundaryPolygon, QGraphicsItem *parent)
47  : QGraphicsPolygonItem(PolygonUtils::Counterclockwise(PolygonUtils::OpenPolygon(boundaryPolygon)), parent)
48 {
49  setFlag(QGraphicsItem::ItemClipsChildrenToShape, true);
50 }
51 
53 {
54 }
55 
57 {
58  std::list<QPolygonF> clippedPolygons;
59  for (const auto &cellPolygon : m_cell_polygons) {
60  clippedPolygons.splice(clippedPolygons.end(), PolygonUtils::ClipPolygon(cellPolygon, polygon()));
61  }
62 
63  return clippedPolygons;
64 }
65 
66 void VoronoiTesselation::mousePressEvent(QGraphicsSceneMouseEvent *event)
67 {
68  QGraphicsPolygonItem::mousePressEvent(event);
69 
70  if (event->button() == Qt::RightButton) {
71  QGraphicsEllipseItem *point = dynamic_cast<QGraphicsEllipseItem*>(scene()->itemAt(event->scenePos(), QTransform()));
72  if (point) {
73  auto it = std::find(m_points.begin(), m_points.end(), point->pos());
74  assert(it != m_points.end() && "m_points should contain the point that has to be removed");
75  m_points.erase(it);
76 
77  delete point;
78 
79  UpdateTesselation();
80  }
81  }
82 }
83 
84 void VoronoiTesselation::mouseDoubleClickEvent(QGraphicsSceneMouseEvent *event)
85 {
86  QGraphicsPolygonItem::mouseDoubleClickEvent(event);
87 
88  if (event->button() == Qt::LeftButton) {
89  m_points.push_back(event->pos());
90 
91  QGraphicsEllipseItem *point = new QGraphicsEllipseItem(g_point_rect, this);
92  point->setPos(event->pos());
93  point->setBrush(g_point_brush);
94  point->setPen(Qt::NoPen);
95 
96  UpdateTesselation();
97  }
98 }
99 
100 void VoronoiTesselation::UpdateTesselation()
101 {
102  std::vector<point_data<long>> points;
103  std::transform(m_points.begin(), m_points.end(), std::back_inserter(points), &ToVoronoiPoint);
104 
105  voronoi_diagram<double> voronoi;
106  construct_voronoi(points.begin(), points.end(), &voronoi);
107 
108  m_cell_polygons.clear();
109 
110  std::map<const voronoi_diagram<double>::edge_type*, ClippedEdge> clippingCache;
111 
112  for (auto &cell : voronoi.cells()) {
113  auto edge = cell.incident_edge();
114  if (!edge) // No edges
115  continue;
116 
117  QPolygonF polygon;
118  do {
119  AddEdge(edge, polygon, clippingCache);
120  edge = edge->next();
121  } while (edge != cell.incident_edge());
122  m_cell_polygons.push_back(polygon);
123  }
124 
125  RedrawLines();
126 }
127 
128 void VoronoiTesselation::RedrawLines()
129 {
130  std::for_each(m_lines.begin(), m_lines.end(), [] (QGraphicsLineItem *item) { delete item; });
131  m_lines.clear();
132 
133  // Set of line already drawn, so we wont draw two line at the same place (which is ugly, when antialiassed)
134  std::hash<double> hash_qreal;
135  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()); };
136  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()); };
137  std::unordered_set<QLineF, decltype(hash_qlinef), decltype(equal_qlinef)> drawnLines(10, hash_qlinef, equal_qlinef);
138 
139  for (const auto &polygon : m_cell_polygons) {
140  std::vector<QPointF> polygonVector = polygon.toStdVector();
141 
142  auto cit = make_const_circular(polygonVector.begin(), polygonVector.end(), polygonVector.begin());
143  do {
144  QLineF line(*cit, *next(cit));
145  if (drawnLines.count(line) == 0) {
146  drawnLines.insert(line);
147  QGraphicsLineItem *lineItem = new QGraphicsLineItem(line, this);
148  lineItem->setZValue(-1);
149  m_lines.push_back(lineItem);
150  }
151  } while (++cit != polygonVector.begin());
152  }
153 }
154 
155 void VoronoiTesselation::AddEdge(const voronoi_diagram<double>::edge_type *edge, QPolygonF &cellPolygon, std::map<const voronoi_diagram<double>::edge_type*, ClippedEdge> &clippingCache)
156 {
157  ClippedEdge clippedEdge = GetClippedEdge(edge, clippingCache);
158 
159  // Check if edge not completely outside of the bounding rect, else this should be handled by a previous or future call to AddEdge
160  if (clippedEdge.clippedP1 || clippedEdge.clippedP2 || boundingRect().contains(clippedEdge.edge.p1())) {
161  cellPolygon.append(clippedEdge.edge.p1());
162 
163  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
164  cellPolygon.append(clippedEdge.edge.p2());
165 
166  auto reenterEdge = edge->next();
167  ClippedEdge reenterClipped;
168  do {
169  reenterClipped = GetClippedEdge(reenterEdge, clippingCache);
170  assert((reenterClipped.clippedP1 || reenterEdge != edge) && "Can't find edge reentering the bounding rect");
171  reenterEdge = reenterEdge->next();
172  } while (!reenterClipped.clippedP1);
173 
174  // Add points in between of exiting and reentering of edge
175  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)
176  assert(!PolygonUtils::IsClockwise(QPolygonF::fromStdVector(boundingPolygon)) && "Qt, get your documentation right!");
177 
178  auto nextBoundaryPoint = std::find(boundingPolygon.begin(), boundingPolygon.end(), clippedEdge.nextBoundaryP2);
179  auto reenterBoundaryPoint = std::find(boundingPolygon.begin(), boundingPolygon.end(), reenterClipped.nextBoundaryP1);
180  assert(nextBoundaryPoint != boundingPolygon.end() && "nextBoundaryPoint not found in boundingRect polygon");
181  assert(reenterBoundaryPoint != boundingPolygon.end() && "reenterBoundaryPoint not found in boundingRect polygon");
182 
183  auto cit = make_const_circular(boundingPolygon.begin(), boundingPolygon.end(), nextBoundaryPoint);
184  while (cit != reenterBoundaryPoint) {
185  cellPolygon.append(*cit);
186  ++cit;
187  }
188  }
189  }
190 }
191 
192 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)
193 {
194  if (clippingCache.find(edge) != clippingCache.end()) {
195  return clippingCache[edge];
196  } else if (clippingCache.find(edge->twin()) != clippingCache.end()) {
197  auto reverse = clippingCache[edge->twin()];
198  return ClippedEdge{QLineF(reverse.edge.p2(), reverse.edge.p1()), reverse.clippedP2, reverse.nextBoundaryP2, reverse.clippedP1, reverse.nextBoundaryP1};
199  } else {
200  auto clipped = ClipEdge(edge);
201  clippingCache[edge] = clipped;
202  return clipped;
203  }
204 }
205 
206 VoronoiTesselation::ClippedEdge VoronoiTesselation::ClipEdge(const boost::polygon::voronoi_diagram<double>::edge_type *edge)
207 {
208  QPointF c1 = m_points[edge->cell()->source_index()];
209  QPointF c2 = m_points[edge->twin()->cell()->source_index()];
210  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
211  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
212 
213  double boundingDiagonal = QLineF(boundingRect().topLeft(), boundingRect().bottomRight()).length();
214 
215  ClippedEdge clipped{};
216  if (edge->vertex0() && edge->vertex1()) { // both vertex0 and vertex
217  clipped.edge.setP1(FromVoronoiPoint(*edge->vertex0()));
218  clipped.edge.setP2(FromVoronoiPoint(*edge->vertex1()));
219  } else if (edge->vertex0()) { // vertex1 is located at infinity
220  clipped.edge.setP1(FromVoronoiPoint(*edge->vertex0()));
221  clipped.edge.setP2(clipped.edge.p1() + (boundingDiagonal + QLineF(clipped.edge.p1(), boundingRect().topLeft()).length()) * direction);
222  } else if (edge->vertex1()) { // vertex0 is located at infinity
223  clipped.edge.setP2(FromVoronoiPoint(*edge->vertex1()));
224  clipped.edge.setP1(clipped.edge.p2() - (boundingDiagonal + QLineF(clipped.edge.p2(), boundingRect().topLeft()).length()) * direction);
225  } else { // vertex0 and vertex1 are both located at infinity
226  QPointF center = (c1 + c2) / 2;
227  clipped.edge.setP1(center - boundingDiagonal * direction);
228  clipped.edge.setP2(center + boundingDiagonal * direction);
229  }
230 
231  std::vector<std::pair<QPointF, QPointF>> intersections;
232  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)
233  assert(!PolygonUtils::IsClockwise(QPolygonF::fromStdVector(boundingPolygon)) && "Qt, get your documentation right!");
234 
235  auto cit = make_const_circular(boundingPolygon.begin(), boundingPolygon.end(), boundingPolygon.begin());
236  do {
237  QPointF p;
238  if (QLineF(*cit, *next(cit)).intersect(clipped.edge, &p) == QLineF::BoundedIntersection) {
239  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
240  }
241  } while (++cit != boundingPolygon.begin());
242 
243  if (intersections.size() > 0) {
244  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(); });
245  auto itUnique = std::unique(intersections.begin(), intersections.end());
246  size_t uniqueIntersections = std::distance(intersections.begin(), itUnique);
247 
248  assert(uniqueIntersections <= 2 && "More than two intersections of a line with a rectangle found");
249  assert(uniqueIntersections >= static_cast<size_t>(!edge->vertex0() + !edge->vertex1()) && "An edge should at least be clipped at every point at infinity");
250 
251  if (uniqueIntersections == 1) {
252  if (!boundingRect().contains(clipped.edge.p1())) { // p2 completely inside of the boudningRect (not at the edges, where it could be intersected with)
253  clipped.clippedP1 = true;
254  clipped.edge.setP1(intersections[0].first);
255  clipped.nextBoundaryP1 = intersections[0].second;
256  } else if (!boundingRect().contains(clipped.edge.p2())) {
257  clipped.clippedP2 = true;
258  clipped.edge.setP2(intersections[0].first);
259  clipped.nextBoundaryP2 = intersections[0].second;
260  }
261  } else if (uniqueIntersections == 2) {
262  clipped.clippedP1 = clipped.edge.p1() != intersections[0].first;
263  clipped.clippedP2 = clipped.edge.p2() != intersections[1].first;
264  clipped.edge.setP1(intersections[0].first);
265  clipped.edge.setP2(intersections[1].first);
266  clipped.nextBoundaryP1 = intersections[0].second;
267  clipped.nextBoundaryP2 = intersections[1].second;
268 
269  }
270  }
271 
272  return clipped;
273 }
274 
275 point_data<long> VoronoiTesselation::ToVoronoiPoint(const QPointF &point)
276 {
277  return point_data<long>(std::roundl(point.x() * g_points_precision), std::roundl(point.y() * g_points_precision));
278 }
279 
280 QPointF VoronoiTesselation::FromVoronoiPoint(const voronoi_diagram<double>::vertex_type &point)
281 {
282  return QPointF(point.x() / g_points_precision, point.y() / g_points_precision);
283 }
284 
285 } // namespace
static QPolygonF OpenPolygon(const QPolygonF &polygon)
Opens the polygon, ie.
Namespace for SimPT tissue editor package.
Definition: Cell.h:32
Combo header for circular iterator.
virtual ~VoronoiTesselation()
Destructor.
static std::list< QPolygonF > ClipPolygon(const QPolygonF &polygon1, const QPolygonF &polygon2)
Clip a polygon with another polygon and return the intersecting subpolygons.
see the online Qt documentation
see the online Qt documentation
static bool IsClockwise(const QPolygonF &polygon)
Checks whether the vertices of a polygon are ordered clockwise or counterclockwise.
virtual void mouseDoubleClickEvent(QGraphicsSceneMouseEvent *event)
Reimplementation of QGraphicsItem::mouseDoubleClickEvent to add a point to the tesselation.
see the online Qt documentation
Namespace for container related classes.
virtual void mousePressEvent(QGraphicsSceneMouseEvent *event)
Reimplementation of QGraphicsItem::mousePressEvent to remove a point from the tesselation.
Interface for VoronoiTesselation.
Functions to manipulate polygons.
Definition: PolygonUtils.h:30
CircularIterator< T > next(CircularIterator< T > i)
Helper yields the position the iterator would have if moved forward (in circular fashion) by 1 positi...
std::list< QPolygonF > GetCellPolygons()
Returns the polygons of cell of the Voronoi tesselation.
ConstCircularIterator< typename T::const_iterator > make_const_circular(const T &c)
Helper produces const circular iterator whose range corresponds to the begin and end iterators of a c...
Interface for PolygonUtils.
see the online Qt documentation