LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
voronoi2d::VoronoiDiagram Class Reference

VoronoiDiagram class definiton. More...

#include "Voronoi.h"

Public Types

using PointPair = std::pair< dcel2d::Point, dcel2d::Point >
 
using MinMaxPointPair = std::pair< PointPair, PointPair >
 

Public Member Functions

 VoronoiDiagram (dcel2d::HalfEdgeList &, dcel2d::VertexList &, dcel2d::FaceList &)
 Constructor. More...
 
virtual ~VoronoiDiagram ()
 Destructor. More...
 
void buildVoronoiDiagram (const dcel2d::PointList &)
 Given an input set of 2D points construct a 2D voronoi diagram. More...
 
void buildVoronoiDiagramBoost (const dcel2d::PointList &)
 
const dcel2d::FaceListgetFaceList () const
 Recover the list of faces. More...
 
const dcel2d::VertexListgetVertexList () const
 Recover the list of vertices. More...
 
double getVoronoiDiagramArea () const
 recover the area of the convex hull More...
 
const dcel2d::PointListgetConvexHull () const
 recover the point list representing the convex hull More...
 
PointPair getExtremePoints () const
 Given an input Point, find the nearest edge. More...
 
PointPair findNearestEdge (const dcel2d::Point &, double &) const
 Given an input Point, find the nearest edge. More...
 
double findNearestDistance (const dcel2d::Point &) const
 Given an input point, find the distance to the nearest edge/point. More...
 

Private Types

using EventQueue = std::priority_queue< IEvent *, std::vector< IEvent * >, bool(*)(const IEvent *, const IEvent *)>
 
using BoostEdgeToEdgeMap = std::map< const boost::polygon::voronoi_edge< double > *, dcel2d::HalfEdge * >
 Translate boost to dcel. More...
 
using BoostVertexToVertexMap = std::map< const boost::polygon::voronoi_vertex< double > *, dcel2d::Vertex * >
 
using BoostCellToFaceMap = std::map< const boost::polygon::voronoi_cell< double > *, dcel2d::Face * >
 

Private Member Functions

void handleSiteEvents (BeachLine &, EventQueue &, IEvent *)
 There are two types of events in the queue, here we handle site events. More...
 
void handleCircleEvents (BeachLine &, EventQueue &, IEvent *)
 There are two types of events in the queue, here we handle circle events. More...
 
void makeLeftCircleEvent (EventQueue &, BSTNode *, double)
 
void makeRightCircleEvent (EventQueue &, BSTNode *, double)
 
IEventmakeCircleEvent (BSTNode *, BSTNode *, BSTNode *, double)
 There are two types of events in the queue, here we handle circle events. More...
 
bool computeCircleCenter (const dcel2d::Coords &, const dcel2d::Coords &, const dcel2d::Coords &, dcel2d::Coords &, double &, double &) const
 
bool computeCircleCenter2 (const dcel2d::Coords &, const dcel2d::Coords &, const dcel2d::Coords &, dcel2d::Coords &, double &, double &) const
 
bool computeCircleCenter3 (const dcel2d::Coords &, const dcel2d::Coords &, const dcel2d::Coords &, dcel2d::Coords &, double &, double &) const
 
void getConvexHull (const BSTNode *)
 this function recovers the convex hull More...
 
void terminateInfiniteEdges (BeachLine &, double)
 this aims to process remaining elements in the beachline after the event queue has been depleted More...
 
bool isInsideConvexHull (const dcel2d::Vertex &) const
 Is a vertex inside the convex hull - meant to be a fast check. More...
 
bool isOutsideConvexHull (const dcel2d::Vertex &, dcel2d::PointList::const_iterator, dcel2d::Coords &, double &) const
 is vertex outside the convex hull and if so return some useful information More...
 
void findBoundingBox (const dcel2d::VertexList &)
 Find the min/max values in x-y to use as a bounding box. More...
 
void boostTranslation (const dcel2d::PointList &, const boost::polygon::voronoi_edge< double > *, const boost::polygon::voronoi_edge< double > *, BoostEdgeToEdgeMap &, BoostVertexToVertexMap &, BoostCellToFaceMap &)
 
void mergeDegenerateVertices ()
 merge degenerate vertices (found by zero length edges) More...
 
double ComputeFaceArea ()
 Compute the area of the faces. More...
 
double crossProduct (const dcel2d::Point &p0, const dcel2d::Point &p1, const dcel2d::Point &p2) const
 Gets the cross product of line from p0 to p1 and p0 to p2. More...
 
double Area () const
 Computes the area of the enclosed convext hull. More...
 
bool isLeft (const dcel2d::Point &p0, const dcel2d::Point &p1, const dcel2d::Point &pCheck) const
 Determines whether a point is to the left, right or on line specifiec by p0 and p1. More...
 

Private Attributes

dcel2d::HalfEdgeListfHalfEdgeList
 
dcel2d::VertexListfVertexList
 
dcel2d::FaceListfFaceList
 
dcel2d::PointList fPointList
 
SiteEventList fSiteEventList
 
CircleEventList fCircleEventList
 
BSTNodeList fCircleNodeList
 
dcel2d::PointList fConvexHullList
 
dcel2d::Coords fConvexHullCenter
 
double fXMin
 
double fXMax
 
double fYMin
 
double fYMax
 
int fNumBadCircles
 
double fVoronoiDiagramArea
 
EventUtilities fUtilities
 

Detailed Description

VoronoiDiagram class definiton.

Definition at line 33 of file Voronoi.h.

Member Typedef Documentation

using voronoi2d::VoronoiDiagram::BoostCellToFaceMap = std::map<const boost::polygon::voronoi_cell<double>*, dcel2d::Face*>
private

Definition at line 153 of file Voronoi.h.

using voronoi2d::VoronoiDiagram::BoostEdgeToEdgeMap = std::map<const boost::polygon::voronoi_edge<double>*, dcel2d::HalfEdge*>
private

Translate boost to dcel.

Definition at line 151 of file Voronoi.h.

using voronoi2d::VoronoiDiagram::BoostVertexToVertexMap = std::map<const boost::polygon::voronoi_vertex<double>*, dcel2d::Vertex*>
private

Definition at line 152 of file Voronoi.h.

using voronoi2d::VoronoiDiagram::EventQueue = std::priority_queue<IEvent*, std::vector<IEvent*>, bool(*)(const IEvent*,const IEvent*)>
private

Definition at line 97 of file Voronoi.h.

Definition at line 37 of file Voronoi.h.

Definition at line 36 of file Voronoi.h.

Constructor & Destructor Documentation

voronoi2d::VoronoiDiagram::VoronoiDiagram ( dcel2d::HalfEdgeList halfEdgeList,
dcel2d::VertexList vertexList,
dcel2d::FaceList faceList 
)

Constructor.

Parameters
pset

Definition at line 53 of file Voronoi.cxx.

References Area(), fCircleEventList, fCircleNodeList, fConvexHullList, fFaceList, fHalfEdgeList, fPointList, fSiteEventList, fVertexList, and fVoronoiDiagramArea.

53  :
54  fHalfEdgeList(halfEdgeList),
55  fVertexList(vertexList),
56  fFaceList(faceList),
57  fXMin(0.),
58  fXMax(0.),
59  fYMin(0.),
60  fYMax(0.),
62 {
63  fHalfEdgeList.clear();
64  fVertexList.clear();
65  fFaceList.clear();
66  fPointList.clear();
67  fSiteEventList.clear();
68  fCircleEventList.clear();
69  fCircleNodeList.clear();
70  fConvexHullList.clear();
71 
72  // And the area
74 }
dcel2d::HalfEdgeList & fHalfEdgeList
Definition: Voronoi.h:187
double Area() const
Computes the area of the enclosed convext hull.
Definition: Voronoi.cxx:106
dcel2d::FaceList & fFaceList
Definition: Voronoi.h:189
BSTNodeList fCircleNodeList
Definition: Voronoi.h:194
CircleEventList fCircleEventList
Definition: Voronoi.h:193
SiteEventList fSiteEventList
Definition: Voronoi.h:192
dcel2d::VertexList & fVertexList
Definition: Voronoi.h:188
dcel2d::PointList fPointList
Definition: Voronoi.h:191
dcel2d::PointList fConvexHullList
Definition: Voronoi.h:196
voronoi2d::VoronoiDiagram::~VoronoiDiagram ( )
virtual

Destructor.

Definition at line 78 of file Voronoi.cxx.

79 {
80 }

Member Function Documentation

double voronoi2d::VoronoiDiagram::Area ( ) const
private

Computes the area of the enclosed convext hull.

Definition at line 106 of file Voronoi.cxx.

References crossProduct(), and fPointList.

Referenced by VoronoiDiagram().

107 {
108  double area(0.);
109 
110  // Compute the area by taking advantage of
111  // 1) the ability to decompose a convex hull into triangles,
112  // 2) the ability to use the cross product to calculate the area
113  // So, the technique is to pick a point (for technical reasons we use 0,0)
114  // and then sum the signed area of triangles formed from this point to two adjecent
115  // vertices on the convex hull.
116  dcel2d::Point center(0.,0.,NULL);
117  dcel2d::Point lastPoint = fPointList.front();
118 
119  for(const auto& point : fPointList)
120  {
121  if (point != lastPoint) area += 0.5 * crossProduct(center,lastPoint,point);
122 
123  lastPoint = point;
124  }
125 
126  return area;
127 }
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:34
dcel2d::PointList fPointList
Definition: Voronoi.h:191
double crossProduct(const dcel2d::Point &p0, const dcel2d::Point &p1, const dcel2d::Point &p2) const
Gets the cross product of line from p0 to p1 and p0 to p2.
Definition: Voronoi.cxx:93
void voronoi2d::VoronoiDiagram::boostTranslation ( const dcel2d::PointList pointList,
const boost::polygon::voronoi_edge< double > *  edge,
const boost::polygon::voronoi_edge< double > *  twin,
BoostEdgeToEdgeMap boostEdgeToEdgeMap,
BoostVertexToVertexMap boostVertexToVertexMap,
BoostCellToFaceMap boostCellToFaceMap 
)
private

Definition at line 323 of file Voronoi.cxx.

References fFaceList, fHalfEdgeList, fVertexList, dcel2d::HalfEdge::setFace(), dcel2d::HalfEdge::setLastHalfEdge(), dcel2d::HalfEdge::setNextHalfEdge(), dcel2d::HalfEdge::setTargetVertex(), and dcel2d::HalfEdge::setTwinHalfEdge().

Referenced by buildVoronoiDiagramBoost().

329 {
330  dcel2d::HalfEdge* halfEdge = NULL;
331  dcel2d::HalfEdge* twinEdge = NULL;
332 
333  if (boostEdgeToEdgeMap.find(edge) != boostEdgeToEdgeMap.end()) halfEdge = boostEdgeToEdgeMap.at(edge);
334  else
335  {
336  fHalfEdgeList.emplace_back();
337 
338  halfEdge = &fHalfEdgeList.back();
339 
340  boostEdgeToEdgeMap[edge] = halfEdge;
341  }
342 
343  if (boostEdgeToEdgeMap.find(twin) != boostEdgeToEdgeMap.end()) twinEdge = boostEdgeToEdgeMap.at(twin);
344  else
345  {
346  fHalfEdgeList.emplace_back();
347 
348  twinEdge = &fHalfEdgeList.back();
349 
350  boostEdgeToEdgeMap[twin] = twinEdge;
351  }
352 
353  // Do the primary half edge first
354  const boost::polygon::voronoi_vertex<double>* boostVertex = edge->vertex1();
355  dcel2d::Vertex* vertex = NULL;
356 
357  // note we can have a null vertex (infinite edge)
358  if (boostVertex)
359  {
360  if (boostVertexToVertexMap.find(boostVertex) == boostVertexToVertexMap.end())
361  {
362  dcel2d::Coords coords(boostVertex->y(),boostVertex->x(),0.);
363 
364  fVertexList.emplace_back(coords, halfEdge);
365 
366  vertex = &fVertexList.back();
367 
368  boostVertexToVertexMap[boostVertex] = vertex;
369  }
370  else vertex = boostVertexToVertexMap.at(boostVertex);
371  }
372 
373  const boost::polygon::voronoi_cell<double>* boostCell = edge->cell();
374  dcel2d::Face* face = NULL;
375 
376  if (boostCellToFaceMap.find(boostCell) == boostCellToFaceMap.end())
377  {
378  dcel2d::PointList::const_iterator pointItr = pointList.begin();
379  int pointIdx = boostCell->source_index();
380 
381  std::advance(pointItr, pointIdx);
382 
383  const dcel2d::Point& point = *pointItr;
384  dcel2d::Coords coords(std::get<0>(point),std::get<1>(point),0.);
385 
386  fFaceList.emplace_back(halfEdge,coords,std::get<2>(point));
387 
388  face = &fFaceList.back();
389 
390  boostCellToFaceMap[boostCell] = face;
391  }
392 
393  halfEdge->setTargetVertex(vertex);
394  halfEdge->setFace(face);
395  halfEdge->setTwinHalfEdge(twinEdge);
396 
397  // For the prev/next half edges we can have two cases, so check:
398  if (boostEdgeToEdgeMap.find(edge->next()) != boostEdgeToEdgeMap.end())
399  {
400  dcel2d::HalfEdge* nextEdge = boostEdgeToEdgeMap.at(edge->next());
401 
402  halfEdge->setNextHalfEdge(nextEdge);
403  nextEdge->setLastHalfEdge(halfEdge);
404  }
405 
406  if (boostEdgeToEdgeMap.find(edge->prev()) != boostEdgeToEdgeMap.end())
407  {
408  dcel2d::HalfEdge* lastEdge = boostEdgeToEdgeMap.at(edge->prev());
409 
410  halfEdge->setLastHalfEdge(lastEdge);
411  lastEdge->setNextHalfEdge(halfEdge);
412  }
413 
414  return;
415 }
dcel2d::HalfEdgeList & fHalfEdgeList
Definition: Voronoi.h:187
dcel2d::FaceList & fFaceList
Definition: Voronoi.h:189
void setLastHalfEdge(HalfEdge *last)
Definition: DCEL.h:167
void setTwinHalfEdge(HalfEdge *twin)
Definition: DCEL.h:165
intermediate_table::const_iterator const_iterator
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:34
dcel2d::VertexList & fVertexList
Definition: Voronoi.h:188
void setFace(Face *face)
Definition: DCEL.h:164
void setNextHalfEdge(HalfEdge *next)
Definition: DCEL.h:166
Eigen::Vector3f Coords
Definition: DCEL.h:36
void setTargetVertex(Vertex *vertex)
Definition: DCEL.h:163
vertex reconstruction
void voronoi2d::VoronoiDiagram::buildVoronoiDiagram ( const dcel2d::PointList pointList)

Given an input set of 2D points construct a 2D voronoi diagram.

Parameters
PointListThe input list of 2D points

Definition at line 135 of file Voronoi.cxx.

References voronoi2d::compareSiteEventPtrs(), voronoi2d::BeachLine::countLeaves(), fCircleEventList, fCircleNodeList, fFaceList, fHalfEdgeList, findBoundingBox(), fNumBadCircles, fPointList, fSiteEventList, fVertexList, fXMax, fXMin, fYMax, fYMin, dcel2d::HalfEdge::getFace(), dcel2d::HalfEdge::getLastHalfEdge(), dcel2d::HalfEdge::getNextHalfEdge(), handleCircleEvents(), handleSiteEvents(), terminateInfiniteEdges(), and voronoi2d::BeachLine::traverseBeach().

Referenced by lar_cluster3d::ClusterPathFinder::buildVoronoiDiagram(), and lar_cluster3d::VoronoiPathFinder::buildVoronoiDiagram().

136 {
137  // Insure all the local data structures have been cleared
138  fHalfEdgeList.clear();
139  fVertexList.clear();
140  fFaceList.clear();
141  fPointList.clear();
142  fSiteEventList.clear();
143  fCircleEventList.clear();
144  fCircleNodeList.clear();
145  fNumBadCircles = 0;
146 
147  std::cout << "******************************************************************************************************************" << std::endl;
148  std::cout << "******************************************************************************************************************" << std::endl;
149  std::cout << "******************************************************************************************************************" << std::endl;
150  std::cout << "==> # input points: " << pointList.size() << std::endl;
151 
152  // Define the priority queue to contain our events
153  EventQueue eventQueue(compareSiteEventPtrs);
154 
155  // Now populate the event queue with site events
156  for(const auto& point : pointList)
157  {
158  fSiteEventList.emplace_back(point);
159  IEvent* iEvent = &fSiteEventList.back();
160  eventQueue.push(iEvent);
161  }
162 
163  // Declare the beachline which will contain the BSTNode objects for site events
164  BeachLine beachLine;
165 
166  // Finally, we need a container for our circle event BSTNodes
167  BSTNodeList circleNodeList;
168 
169  // Now process the queue
170  while(!eventQueue.empty())
171  {
172  IEvent* event = eventQueue.top();
173 
174  eventQueue.pop();
175 
176  // If a site or circle event then handle appropriately
177  if (event->isSite()) handleSiteEvents(beachLine, eventQueue, event);
178  else if (event->isValid()) handleCircleEvents(beachLine, eventQueue, event);
179  }
180 
181  std::cout << "*******> # input points: " << pointList.size() << ", remaining leaves: " << beachLine.countLeaves() << ", " << beachLine.traverseBeach() << ", # bad circles: " << fNumBadCircles << std::endl;
182  std::cout << " Faces: " << fFaceList.size() << ", Vertices: " << fVertexList.size() << ", # half edges: " << fHalfEdgeList.size() << std::endl;
183 
184  // Get the bounding box
186 
187  std::cout << " Range min/maxes, x: " << fXMin << ", " << fXMax << ", y: " << fYMin << ", " << fYMax << std::endl;
188 
189  // Terminate the infinite edges
190  terminateInfiniteEdges(beachLine, std::get<0>(pointList.front()));
191 
192  // Look for open faces
193  int nOpenFaces(0);
194 
195  std::map<int,int> edgeCountMap;
196  std::map<int,int> openCountMap;
197 
198  for (const auto& face : fFaceList)
199  {
200  int nEdges(1);
201  bool closed(false);
202  const dcel2d::HalfEdge* startEdge = face.getHalfEdge();
203 
204  // Count forwards
205  for(const dcel2d::HalfEdge* halfEdge = startEdge->getNextHalfEdge(); halfEdge && !closed;)
206  {
207  if (halfEdge->getFace() != &face)
208  {
209  std::cout << " ===> halfEdge does not match face: " << halfEdge << ", face: " << halfEdge->getFace() << ", base: " << &face << std::endl;
210  }
211 
212  if (halfEdge == startEdge)
213  {
214  closed = true;
215  break;
216  }
217  nEdges++;
218  halfEdge = halfEdge->getNextHalfEdge();
219  }
220 
221  // Count backwards. Note if face is closed then nothing to do here
222  if (!closed)
223  {
224  for(const dcel2d::HalfEdge* halfEdge = startEdge->getLastHalfEdge(); halfEdge && !closed;)
225  {
226  if (halfEdge->getFace() != &face)
227  {
228  std::cout << " ===> halfEdge does not match face: " << halfEdge << ", face: " << halfEdge->getFace() << ", base: " << &face << std::endl;
229  }
230 
231  if (halfEdge == startEdge)
232  {
233  closed = true;
234  break;
235  }
236  nEdges++;
237  halfEdge = halfEdge->getLastHalfEdge();
238  }
239  }
240 
241  if (!closed)
242  {
243  nOpenFaces++;
244  openCountMap[nEdges]++;
245  }
246 
247  edgeCountMap[nEdges]++;
248  }
249 
250  std::cout << "==> Found " << nOpenFaces << " open faces from total of " << fFaceList.size() << std::endl;
251  for(const auto& edgeCount : edgeCountMap) std::cout << " -> all edges, # edges: " << edgeCount.first << ", count: " << edgeCount.second << std::endl;
252  for(const auto& edgeCount : openCountMap) std::cout << " -> open edges, # edges: " << edgeCount.first << ", count: " << edgeCount.second << std::endl;
253  std::cout << "******************************************************************************************************************" << std::endl;
254  std::cout << "******************************************************************************************************************" << std::endl;
255  std::cout << "******************************************************************************************************************" << std::endl;
256 
257  // Clear internal containers that are no longer useful
258  fSiteEventList.clear();
259  fCircleEventList.clear();
260  fCircleNodeList.clear();
261 
262  return;
263 }
HalfEdge * getLastHalfEdge() const
Definition: DCEL.h:161
HalfEdge * getNextHalfEdge() const
Definition: DCEL.h:160
dcel2d::HalfEdgeList & fHalfEdgeList
Definition: Voronoi.h:187
Face * getFace() const
Definition: DCEL.h:158
dcel2d::FaceList & fFaceList
Definition: Voronoi.h:189
std::list< BSTNode > BSTNodeList
Definition: BeachLine.h:122
BSTNodeList fCircleNodeList
Definition: Voronoi.h:194
void handleCircleEvents(BeachLine &, EventQueue &, IEvent *)
There are two types of events in the queue, here we handle circle events.
Definition: Voronoi.cxx:474
CircleEventList fCircleEventList
Definition: Voronoi.h:193
SiteEventList fSiteEventList
Definition: Voronoi.h:192
std::priority_queue< IEvent *, std::vector< IEvent * >, bool(*)(const IEvent *, const IEvent *)> EventQueue
Definition: Voronoi.h:97
bool compareSiteEventPtrs(const IEvent *left, const IEvent *right)
Definition: Voronoi.cxx:131
void findBoundingBox(const dcel2d::VertexList &)
Find the min/max values in x-y to use as a bounding box.
Definition: Voronoi.cxx:1348
dcel2d::VertexList & fVertexList
Definition: Voronoi.h:188
dcel2d::PointList fPointList
Definition: Voronoi.h:191
void terminateInfiniteEdges(BeachLine &, double)
this aims to process remaining elements in the beachline after the event queue has been depleted ...
Definition: Voronoi.cxx:831
Event finding and building.
void handleSiteEvents(BeachLine &, EventQueue &, IEvent *)
There are two types of events in the queue, here we handle site events.
Definition: Voronoi.cxx:417
void voronoi2d::VoronoiDiagram::buildVoronoiDiagramBoost ( const dcel2d::PointList pointList)

Definition at line 267 of file Voronoi.cxx.

References boostTranslation(), fCircleEventList, fCircleNodeList, fFaceList, fHalfEdgeList, fNumBadCircles, fPointList, fSiteEventList, and fVertexList.

268 {
269  // Insure all the local data structures have been cleared
270  fHalfEdgeList.clear();
271  fVertexList.clear();
272  fFaceList.clear();
273  fPointList.clear();
274  fSiteEventList.clear();
275  fCircleEventList.clear();
276  fCircleNodeList.clear();
277  fNumBadCircles = 0;
278 
279  std::cout << "******************************************************************************************************************" << std::endl;
280  std::cout << "******************************************************************************************************************" << std::endl;
281  std::cout << "******************************************************************************************************************" << std::endl;
282  std::cout << "==> # input points: " << pointList.size() << std::endl;
283 
284  // Construct out voronoi diagram
285  boost::polygon::voronoi_diagram<double> vd;
286  boost::polygon::construct_voronoi(pointList.begin(),pointList.end(),&vd);
287 
288  // Make maps for translating from boost to me (or we can rewrite our code in boost... for now maps)
289  using BoostEdgeToEdgeMap = std::map<const boost::polygon::voronoi_edge<double>*, dcel2d::HalfEdge*>;
290  using BoostVertexToVertexMap = std::map<const boost::polygon::voronoi_vertex<double>*, dcel2d::Vertex*>;
291  using BoostCellToFaceMap = std::map<const boost::polygon::voronoi_cell<double>*, dcel2d::Face*>;
292 
293  BoostEdgeToEdgeMap boostEdgeToEdgeMap;
294  BoostVertexToVertexMap boostVertexToVertexMap;
295  BoostCellToFaceMap boostCellToFaceMap;
296 
297  // Loop over the edges
298  for(const auto& edge : vd.edges())
299  {
300  const boost::polygon::voronoi_edge<double>* twin = edge.twin();
301 
302  boostTranslation(pointList, &edge, twin, boostEdgeToEdgeMap, boostVertexToVertexMap, boostCellToFaceMap);
303  boostTranslation(pointList, twin, &edge, boostEdgeToEdgeMap, boostVertexToVertexMap, boostCellToFaceMap);
304  }
305 
306  //std::cout << "==> Found " << nOpenFaces << " open faces from total of " << fFaceList.size() << std::endl;
307  //for(const auto& edgeCount : edgeCountMap) std::cout << " -> all edges, # edges: " << edgeCount.first << ", count: " << edgeCount.second << std::endl;
308  //for(const auto& edgeCount : openCountMap) std::cout << " -> open edges, # edges: " << edgeCount.first << ", count: " << edgeCount.second << std::endl;
309  std::cout << "==> Boost returns " << fHalfEdgeList.size() << " edges, " << fFaceList.size() << " faces, " << fVertexList.size() << " vertices " << std::endl;
310  std::cout << " " << vd.edges().size() << " edges, " << vd.cells().size() << " faces, " << vd.vertices().size() << " vertices" << std::endl;
311  std::cout << "******************************************************************************************************************" << std::endl;
312  std::cout << "******************************************************************************************************************" << std::endl;
313  std::cout << "******************************************************************************************************************" << std::endl;
314 
315  // Clear internal containers that are no longer useful
316  fSiteEventList.clear();
317  fCircleEventList.clear();
318  fCircleNodeList.clear();
319 
320  return;
321 }
dcel2d::HalfEdgeList & fHalfEdgeList
Definition: Voronoi.h:187
dcel2d::FaceList & fFaceList
Definition: Voronoi.h:189
BSTNodeList fCircleNodeList
Definition: Voronoi.h:194
CircleEventList fCircleEventList
Definition: Voronoi.h:193
SiteEventList fSiteEventList
Definition: Voronoi.h:192
dcel2d::VertexList & fVertexList
Definition: Voronoi.h:188
dcel2d::PointList fPointList
Definition: Voronoi.h:191
void boostTranslation(const dcel2d::PointList &, const boost::polygon::voronoi_edge< double > *, const boost::polygon::voronoi_edge< double > *, BoostEdgeToEdgeMap &, BoostVertexToVertexMap &, BoostCellToFaceMap &)
Definition: Voronoi.cxx:323
std::map< const boost::polygon::voronoi_edge< double > *, dcel2d::HalfEdge * > BoostEdgeToEdgeMap
Translate boost to dcel.
Definition: Voronoi.h:151
std::map< const boost::polygon::voronoi_vertex< double > *, dcel2d::Vertex * > BoostVertexToVertexMap
Definition: Voronoi.h:152
std::map< const boost::polygon::voronoi_cell< double > *, dcel2d::Face * > BoostCellToFaceMap
Definition: Voronoi.h:153
bool voronoi2d::VoronoiDiagram::computeCircleCenter ( const dcel2d::Coords p1,
const dcel2d::Coords p2,
const dcel2d::Coords p3,
dcel2d::Coords center,
double &  radius,
double &  delta 
) const
private

Definition at line 692 of file Voronoi.cxx.

References e, fNumBadCircles, max, x2, and y2.

Referenced by makeCircleEvent().

698 {
699  // The method is to translate the three points to a system where the first point is at the origin. Then we
700  // are looking for a circle that passes through the origin and the two remaining (translated) points. In
701  // this mode the circle radius is the distance from the origin to the circle center which simplifies the
702  // calculation.
703  double xCoord = p1[0];
704  double yCoord = p1[1];
705 
706  double x2 = p2[0] - xCoord;
707  double y2 = p2[1] - yCoord;
708  double x3 = p3[0] - xCoord;
709  double y3 = p3[1] - yCoord;
710 
711  double det = x2 * y3 - x3 * y2;
712 
713  // Points are colinear so cannot make a circle
714  // Or, if negative then the midpoint is "left" of line between first and third points meaning the
715  // circle curvature is the "wrong" way for making a circle event
716  if (det <= 0.) // std::numeric_limits<double>::epsilon())
717  //if (!(std::abs(det) > 0.)) // std::numeric_limits<double>::epsilon())
718  {
719  if (det > -std::numeric_limits<double>::epsilon())
720  std::cout << " --->Circle failure, det: " << det << ", mid x: " << p2[0] << ", y: " << p2[1]
721  << ", l x: " << p1[0] << ", y: " << p1[1]
722  << ", r x: " << p3[0] << ", y: " << p3[1] << std::endl;
723 
724  return false;
725  }
726 
727  double p2sqr = x2 * x2 + y2 * y2;
728  double p3sqr = x3 * x3 + y3 * y3;
729 
730  double cxpr = 0.5 * (y3 * p2sqr - y2 * p3sqr) / det;
731  double cypr = 0.5 * (x2 * p3sqr - x3 * p2sqr) / det;
732 
733  radius = std::sqrt(cxpr * cxpr + cypr * cypr);
734  center[0] = cxpr + xCoord;
735  center[1] = cypr + yCoord;
736  center[2] = 0.;
737 
738  // So... roundoff error can cause degenerate circles to get missed
739  // We need to attack that problem by calculating the radius all possible
740  // ways and then take the largest...
741  dcel2d::Coords p1Rad = p1 - center;
742  dcel2d::Coords p2Rad = p2 - center;
743  dcel2d::Coords p3Rad = p3 - center;
744 
745  std::vector<float> radSqrVec;
746 
747  radSqrVec.push_back(p1Rad.norm());
748  radSqrVec.push_back(p2Rad.norm());
749  radSqrVec.push_back(p3Rad.norm());
750 
751  double maxRadius = *std::max_element(radSqrVec.begin(),radSqrVec.end());
752 
753  delta = std::max(5.e-7, maxRadius - radius);
754 
755  if (radius > 1000.)
756  {
757 // std::cout << " ***> Radius = " << radius << ", circ x,y: " << center[0] << ", " << center[1] << ", p1: " << xCoord << "," << yCoord << ", p2: " << p2[0] << "," << p2[1] << ", p3: " << p3[0] << "," << p3[1] << ", det: " << det << std::endl;
758  fNumBadCircles++;
759  }
760 
761  return true;
762 }
Float_t y2[n_points_geant4]
Definition: compare.C:26
Int_t max
Definition: plot.C:27
Double_t radius
Eigen::Vector3f Coords
Definition: DCEL.h:36
Float_t x2[n_points_geant4]
Definition: compare.C:26
Float_t e
Definition: plot.C:34
bool voronoi2d::VoronoiDiagram::computeCircleCenter2 ( const dcel2d::Coords p1,
const dcel2d::Coords p2,
const dcel2d::Coords p3,
dcel2d::Coords center,
double &  radius,
double &  delta 
) const
private

Definition at line 764 of file Voronoi.cxx.

770 {
771  // Compute the circle center as the intersection of the two perpendicular bisectors of rays between the points
772  double slope12 = (p2[1] - p1[1]) / (p2[0] - p1[0]);
773  double slope32 = (p3[1] - p2[1]) / (p3[0] - p2[0]);
774 
775  if (std::abs(slope12 - slope32) <= std::numeric_limits<double>::epsilon())
776  {
777  std::cout << " >>>> Matching slopes! points: (" << p1[0] << "," << p1[1] << "), ("<< p2[0] << "," << p2[1] << "), ("<< p3[0] << "," << p3[1] << ")" << std::endl;
778 
779  return false;
780  }
781 
782  center[0] = (slope12 * slope32 * (p3[1] - p1[1]) + slope12 * (p2[0] + p3[0]) - slope32 * (p1[0] + p2[0])) / (2. * (slope12 - slope32));
783  center[1] = 0.5 * (p1[1] + p2[1]) - (center[0] - 0.5 * (p1[0] + p2[0])) / slope12;
784  center[2] = 0.;
785  radius = std::sqrt(std::pow((p2[0] - center[0]),2) + std::pow((p2[1] - center[1]),2));
786 
787  if (radius > 100.)
788  {
789  std::cout << " ***> Rad2 = " << radius << ", circ x,y: " << center[0] << "," << center[1] << ", p1: " << p1[0] << "," << p1[1] << ", p2: " << p2[0] << "," << p2[1] << ", p3: " << p3[0] << ", " << p3[1] << std::endl;
790  }
791 
792  return true;
793 }
Double_t radius
bool voronoi2d::VoronoiDiagram::computeCircleCenter3 ( const dcel2d::Coords p1,
const dcel2d::Coords p2,
const dcel2d::Coords p3,
dcel2d::Coords center,
double &  radius,
double &  delta 
) const
private

Definition at line 795 of file Voronoi.cxx.

801 {
802  // Yet another bisector method to calculate the circle center...
803  double temp = p2[0] * p2[0] + p2[1] * p2[1];
804  double p1p2 = (p1[0] * p1[0] + p1[1] * p1[1] - temp) / 2.0;
805  double p2p3 = (temp - p3[0] * p3[0] - p3[1] * p3[1]) / 2.0;
806  double det = (p1[0] - p2[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p2[1]);
807 
808  if (std::abs(det) < 10. * std::numeric_limits<double>::epsilon())
809  {
810  std::cout << " >>>> Determinant zero! points: (" << p1[0] << "," << p1[1] << "), ("<< p2[0] << "," << p2[1] << "), ("<< p3[0] << "," << p3[1] << ")" << std::endl;
811 
812  return false;
813  }
814 
815  det = 1. / det;
816 
817  center[0] = (p1p2 * (p2[1] - p3[1]) - p2p3 * (p1[1] - p2[1])) * det;
818  center[1] = (p2p3 * (p1[0] - p2[0]) - p1p2 * (p2[0] - p3[0])) * det;
819  center[2] = 0.;
820 
821  radius = std::sqrt(std::pow((p1[0] - center[0]),2) + std::pow((p1[1] - center[1]),2));
822 
823  if (radius > 100.)
824  {
825  std::cout << " ***> Rad3 = " << radius << ", circ x,y: " << center[0] << "," << center[1] << ", p1: " << p1[0] << "," << p1[1] << ", p2: " << p2[0] << "," << p2[1] << ", p3: " << p3[0] << ", " << p3[1] << std::endl;
826  }
827 
828  return true;
829 }
Double_t radius
double voronoi2d::VoronoiDiagram::ComputeFaceArea ( )
private

Compute the area of the faces.

Definition at line 1217 of file Voronoi.cxx.

References crossProduct(), e, fFaceList, dcel2d::Vertex::getCoords(), dcel2d::HalfEdge::getNextHalfEdge(), dcel2d::HalfEdge::getTargetVertex(), dcel2d::HalfEdge::getTwinHalfEdge(), art::left(), max, min, and art::right().

Referenced by terminateInfiniteEdges().

1218 {
1219  // Compute the area by taking advantage of
1220  // 1) the ability to decompose a convex hull into triangles,
1221  // 2) the ability to use the cross product to calculate the area
1222  // So the idea is to loop through all the faces and then follow the edges
1223  // around the face to compute the area of the face.
1224  // Note that a special case are the "infinite faces" which lie on the
1225  // convex hull of the event. Skip these for now...
1226 
1227  double totalArea(0.);
1228  int nNonInfiniteFaces(0);
1229  double smallestArea(std::numeric_limits<double>::max());
1230  double largestArea(0.);
1231 
1232  std::vector<std::pair<double,const dcel2d::Face*>> areaFaceVec;
1233 
1234  areaFaceVec.reserve(fFaceList.size());
1235 
1236  for(auto& face : fFaceList)
1237  {
1238 // const dcel2d::Coords& faceCoords = face.getCoords();
1239  const dcel2d::HalfEdge* halfEdge = face.getHalfEdge();
1240  double faceArea(0.);
1241  int numEdges(0);
1242  bool doNext = true;
1243 
1244  dcel2d::Coords faceCenter(0.,0.,0.);
1245 
1246  while(doNext)
1247  {
1248  if (halfEdge->getTargetVertex())
1249  faceCenter += halfEdge->getTargetVertex()->getCoords();
1250 
1251  numEdges++;
1252 
1253  halfEdge = halfEdge->getNextHalfEdge();
1254 
1255  if (!halfEdge)
1256  {
1257  faceArea = std::numeric_limits<double>::max();
1258  doNext = false;
1259  }
1260 
1261  if (halfEdge == face.getHalfEdge()) doNext = false;
1262  }
1263 
1264  faceCenter /= numEdges;
1265 
1266  halfEdge = face.getHalfEdge();
1267  doNext = true;
1268 
1269  while(doNext)
1270  {
1271  const dcel2d::HalfEdge* twinEdge = halfEdge->getTwinHalfEdge();
1272 
1273  if (!halfEdge->getTargetVertex() || !twinEdge->getTargetVertex())
1274  {
1275  faceArea = std::numeric_limits<double>::max();
1276  break;
1277  }
1278 
1279  // Recover the two vertex points
1280  const dcel2d::Coords& p1 = halfEdge->getTargetVertex()->getCoords();
1281  const dcel2d::Coords& p2 = twinEdge->getTargetVertex()->getCoords();
1282 
1283  // Define a quick 2D cross product here since it will used quite a bit!
1284  double dp1p0X = p1[0] - faceCenter[0];
1285  double dp1p0Y = p1[1] - faceCenter[1];
1286  double dp2p0X = p2[0] - faceCenter[0];
1287  double dp2p0Y = p2[1] - faceCenter[1];
1288 
1289  //faceArea += dp1p0X * dp2p0Y - dp1p0Y * dp2p0X;
1290  double crossProduct = dp1p0X * dp2p0Y - dp1p0Y * dp2p0X;
1291 
1292  faceArea += crossProduct;
1293 
1294  if (crossProduct < 0.)
1295  {
1296  dcel2d::Coords edgeVec = p1 - p2;
1297  std::cout << "--- negative cross: " << crossProduct << ", edgeLen: " << edgeVec.norm() << ", x/y: " << edgeVec[0] << "/" << edgeVec[1] << std::endl;
1298  }
1299 
1300 // numEdges++;
1301 
1302  halfEdge = halfEdge->getNextHalfEdge();
1303 
1304  if (!halfEdge)
1305  {
1306  faceArea = std::numeric_limits<double>::max();
1307  break;
1308  }
1309 
1310  if (halfEdge == face.getHalfEdge()) doNext = false;
1311  }
1312 
1313  areaFaceVec.emplace_back(faceArea,&face);
1314 
1315  if (faceArea < std::numeric_limits<double>::max() && faceArea > 0.)
1316  {
1317  nNonInfiniteFaces++;
1318  totalArea += faceArea;
1319  smallestArea = std::min(faceArea,smallestArea);
1320  largestArea = std::max(faceArea,largestArea);
1321  }
1322 
1323  if (faceArea < 1.e-4) std::cout << "---> face area <= 0: " << faceArea << ", with " << numEdges << " edges" << std::endl;
1324 
1325  face.setFaceArea(faceArea);
1326  }
1327 
1328  // Calculate a truncated mean...
1329  std::sort(areaFaceVec.begin(),areaFaceVec.end(),[](const auto& left, const auto& right){return left.first < right.first;});
1330 
1331  std::vector<std::pair<double,const dcel2d::Face*>>::iterator firstItr = std::find_if(areaFaceVec.begin(),areaFaceVec.end(),[](const auto& val){return val.first > 0.;});
1332  std::vector<std::pair<double,const dcel2d::Face*>>::iterator lastItr = std::find_if(areaFaceVec.begin(),areaFaceVec.end(),[](const auto& val){return !(val.first < std::numeric_limits<double>::max());});
1333 
1334  size_t nToKeep = 0.8 * std::distance(firstItr,lastItr);
1335 
1336  std::cout << ">>>>> nToKeep: " << nToKeep << ", last non infinite entry: " << std::distance(areaFaceVec.begin(),lastItr) << std::endl;
1337 
1338  double totalTruncArea = std::accumulate(firstItr,firstItr+nToKeep,0.,[](auto sum, const auto& val){return sum+val.first;});
1339  double aveTruncArea = totalTruncArea / double(nToKeep);
1340 
1341  if (nNonInfiniteFaces > 0) std::cout << ">>>> Face area for " << nNonInfiniteFaces << ", ave area: " << totalArea / nNonInfiniteFaces << ", ave trunc area: " << aveTruncArea << ", ratio: " << totalTruncArea / totalArea << ", smallest: " << smallestArea << ", largest: " << largestArea << std::endl;
1342  else std::cout << ">>>>> there are no non infinite faces" << std::endl;
1343 
1344  return totalArea;
1345 }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
HalfEdge * getNextHalfEdge() const
Definition: DCEL.h:160
HalfEdge * getTwinHalfEdge() const
Definition: DCEL.h:159
dcel2d::FaceList & fFaceList
Definition: Voronoi.h:189
Int_t max
Definition: plot.C:27
double crossProduct(const dcel2d::Point &p0, const dcel2d::Point &p1, const dcel2d::Point &p2) const
Gets the cross product of line from p0 to p1 and p0 to p2.
Definition: Voronoi.cxx:93
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
Int_t min
Definition: plot.C:26
Vertex * getTargetVertex() const
Definition: DCEL.h:157
Eigen::Vector3f Coords
Definition: DCEL.h:36
Float_t e
Definition: plot.C:34
const Coords & getCoords() const
Definition: DCEL.h:61
double voronoi2d::VoronoiDiagram::crossProduct ( const dcel2d::Point p0,
const dcel2d::Point p1,
const dcel2d::Point p2 
) const
private

Gets the cross product of line from p0 to p1 and p0 to p2.

Definition at line 93 of file Voronoi.cxx.

Referenced by Area(), ComputeFaceArea(), and isLeft().

94 {
95  // Define a quick 2D cross product here since it will used quite a bit!
96  double deltaX = std::get<0>(p1) - std::get<0>(p0);
97  double deltaY = std::get<1>(p1) - std::get<1>(p0);
98  double dCheckX = std::get<0>(p2) - std::get<0>(p0);
99  double dCheckY = std::get<1>(p2) - std::get<1>(p0);
100 
101  return ((deltaX * dCheckY) - (deltaY * dCheckX));
102 }
void voronoi2d::VoronoiDiagram::findBoundingBox ( const dcel2d::VertexList vertexList)
private

Find the min/max values in x-y to use as a bounding box.

Definition at line 1348 of file Voronoi.cxx.

References fXMax, fXMin, fYMax, fYMin, art::left(), and art::right().

Referenced by buildVoronoiDiagram().

1349 {
1350  // Find extremes in x to start
1351  std::pair<dcel2d::VertexList::const_iterator,dcel2d::VertexList::const_iterator> minMaxItrX = std::minmax_element(vertexList.begin(),vertexList.end(),[](const auto& left, const auto& right){return left.getCoords()[0] < right.getCoords()[0];});
1352 
1353  fXMin = minMaxItrX.first->getCoords()[0];
1354  fXMax = minMaxItrX.second->getCoords()[0];
1355 
1356  // To get the extremes in y we need to make a pass through the list
1357  std::pair<dcel2d::VertexList::const_iterator,dcel2d::VertexList::const_iterator> minMaxItrY = std::minmax_element(vertexList.begin(),vertexList.end(),[](const auto& left, const auto& right){return left.getCoords()[1] < right.getCoords()[1];});
1358 
1359  fYMin = minMaxItrY.first->getCoords()[1];
1360  fYMax = minMaxItrY.second->getCoords()[1];
1361 
1362  return;
1363 }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
double voronoi2d::VoronoiDiagram::findNearestDistance ( const dcel2d::Point point) const

Given an input point, find the distance to the nearest edge/point.

Definition at line 1426 of file Voronoi.cxx.

References findNearestEdge().

Referenced by getConvexHull().

1427 {
1428  double closestDistance;
1429 
1430  findNearestEdge(point,closestDistance);
1431 
1432  return closestDistance;
1433 }
PointPair findNearestEdge(const dcel2d::Point &, double &) const
Given an input Point, find the nearest edge.
Definition: Voronoi.cxx:1365
VoronoiDiagram::PointPair voronoi2d::VoronoiDiagram::findNearestEdge ( const dcel2d::Point point,
double &  closestDistance 
) const

Given an input Point, find the nearest edge.

Definition at line 1365 of file Voronoi.cxx.

References fPointList, isLeft(), and max.

Referenced by findNearestDistance(), and getConvexHull().

1366 {
1367  // The idea is to find the nearest edge of the convex hull, defined by
1368  // two adjacent vertices of the hull, to the input point.
1369  // As near as I can tell, the best way to do this is to apply brute force...
1370  // Idea will be to iterate over pairs of points
1371  dcel2d::PointList::const_iterator curPointItr = fPointList.begin();
1372  dcel2d::Point prevPoint = *curPointItr++;
1373  dcel2d::Point curPoint = *curPointItr;
1374 
1375  // Set up the winner
1376  PointPair closestEdge(prevPoint,curPoint);
1377 
1378  closestDistance = std::numeric_limits<double>::max();
1379 
1380  // curPointItr is meant to point to the second point
1381  while(curPointItr != fPointList.end())
1382  {
1383  if (curPoint != prevPoint)
1384  {
1385  // Dereference some stuff
1386  double xPrevToPoint = (std::get<0>(point) - std::get<0>(prevPoint));
1387  double yPrevToPoint = (std::get<1>(point) - std::get<1>(prevPoint));
1388  double xPrevToCur = (std::get<0>(curPoint) - std::get<0>(prevPoint));
1389  double yPrevToCur = (std::get<1>(curPoint) - std::get<1>(prevPoint));
1390  double edgeLength = std::sqrt(xPrevToCur * xPrevToCur + yPrevToCur * yPrevToCur);
1391 
1392  // Find projection onto convex hull edge
1393  double projection = ((xPrevToPoint * xPrevToCur) + (yPrevToPoint * yPrevToCur)) / edgeLength;
1394 
1395  // DOCA point
1396  dcel2d::Point docaPoint(std::get<0>(prevPoint) + projection * xPrevToCur / edgeLength,
1397  std::get<1>(prevPoint) + projection * yPrevToCur / edgeLength, 0);
1398 
1399  if (projection > edgeLength) docaPoint = curPoint;
1400  if (projection < 0) docaPoint = prevPoint;
1401 
1402  double xDocaDist = std::get<0>(point) - std::get<0>(docaPoint);
1403  double yDocaDist = std::get<1>(point) - std::get<1>(docaPoint);
1404  double docaDist = xDocaDist * xDocaDist + yDocaDist * yDocaDist;
1405 
1406  if (docaDist < closestDistance)
1407  {
1408  closestEdge = PointPair(prevPoint,curPoint);
1409  closestDistance = docaDist;
1410  }
1411  }
1412 
1413  prevPoint = curPoint;
1414  curPoint = *curPointItr++;
1415  }
1416 
1417  closestDistance = std::sqrt(closestDistance);
1418 
1419  // Convention is convex hull vertices sorted in counter clockwise fashion so if the point
1420  // lays to the left of the nearest edge then it must be an interior point
1421  if (isLeft(closestEdge.first,closestEdge.second,point)) closestDistance = -closestDistance;
1422 
1423  return closestEdge;
1424 }
Int_t max
Definition: plot.C:27
intermediate_table::const_iterator const_iterator
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:34
dcel2d::PointList fPointList
Definition: Voronoi.h:191
std::pair< dcel2d::Point, dcel2d::Point > PointPair
Definition: Voronoi.h:36
bool isLeft(const dcel2d::Point &p0, const dcel2d::Point &p1, const dcel2d::Point &pCheck) const
Determines whether a point is to the left, right or on line specifiec by p0 and p1.
Definition: Voronoi.cxx:84
const dcel2d::PointList& voronoi2d::VoronoiDiagram::getConvexHull ( ) const
inline

recover the point list representing the convex hull

Definition at line 78 of file Voronoi.h.

References fConvexHullList, findNearestDistance(), findNearestEdge(), and getExtremePoints().

Referenced by terminateInfiniteEdges().

78 {return fConvexHullList;}
dcel2d::PointList fConvexHullList
Definition: Voronoi.h:196
void voronoi2d::VoronoiDiagram::getConvexHull ( const BSTNode topNode)
private

this function recovers the convex hull

Definition at line 965 of file Voronoi.cxx.

References fConvexHullCenter, fConvexHullList, lar_cluster3d::ConvexHull::getConvexHull(), voronoi2d::IEvent::getCoords(), voronoi2d::BSTNode::getEvent(), voronoi2d::BSTNode::getFace(), voronoi2d::BSTNode::getLeftChild(), voronoi2d::IEvent::getPoint(), voronoi2d::BSTNode::getPredecessor(), voronoi2d::BSTNode::getRightChild(), art::left(), art::right(), and dcel2d::Face::setOnConvexHull().

966 {
967  // Assume the input node is the top of the binary search tree and represents
968  // the beach line at the end of the sweep algorithm
969  // The convex hull will then be the leaf elements along the beach line.
970  // Note that this algorithm will give you the convex hull in a clockwise manner
971  if (topNode)
972  {
973  const BSTNode* node = topNode;
974 
975  // Initialize the center
976  fConvexHullCenter = dcel2d::Coords(0.,0.,0.);
977 
978  // Run down the left branches to find the right most leaf
979  // We do it this way to make sure we get a CCW convex hull
980  while(node->getRightChild()) node = node->getRightChild();
981 
982  // Include a check on convexity...
983  Eigen::Vector3f prevVec(0.,0.,0.);
984  dcel2d::Coords lastPoint(0.,0.,0.);
985 
986  // Now we are going to traverse across the beach line and process the remaining leaves
987  // This will identify the convex hull
988  while(node)
989  {
990  if (!node->getLeftChild() && !node->getRightChild())
991  {
992  // Add point to the convex hull list
993  fConvexHullList.emplace_back(node->getEvent()->getPoint());
994 
995  // Mark the face as on the convex hull
996  node->getFace()->setOnConvexHull();
997 
998  dcel2d::Coords nextPoint = node->getEvent()->getCoords();
999  Eigen::Vector3f curVec = nextPoint - lastPoint;
1000 
1001  curVec.normalize();
1002 
1003  double dotProd = prevVec.dot(curVec);
1004 
1005  std::cout << "--> lastPoint: " << lastPoint[0] << "/" << lastPoint[1] << ", tan: " << std::atan2(lastPoint[1],lastPoint[0]) << ", curPoint: " << nextPoint[0] << "/" << nextPoint[1] << ", tan: " << std::atan2(nextPoint[1],nextPoint[0]) << ", dot: " << dotProd << std::endl;
1006 
1007  prevVec = curVec;
1008  lastPoint = nextPoint;
1009  }
1010 
1011  node = node->getPredecessor();
1012  }
1013 
1014  // Annoyingly, our algorithm does not contain only the convex hull points and so we need to skim out the renegades...
1016 
1017  for(const auto& edgePoint : fConvexHullList) localList.emplace_back(std::get<0>(edgePoint),std::get<1>(edgePoint),std::get<2>(edgePoint));
1018 
1019  // Sort the point vec by increasing x, then increase y
1020  localList.sort([](const auto& left, const auto& right){return (std::abs(std::get<0>(left) - std::get<0>(right)) > std::numeric_limits<float>::epsilon()) ? std::get<0>(left) < std::get<0>(right) : std::get<1>(left) < std::get<1>(right);});
1021 
1022  // Why do I need to do this?
1023  lar_cluster3d::ConvexHull convexHull(localList);
1024 
1025  // Clear the convex hull list...
1026 // fConvexHullList.clear();
1027 
1028  std::cout << "~~~>> there are " << convexHull.getConvexHull().size() << " convex hull points and " << fConvexHullList.size() << " infinite cells" << std::endl;
1029 
1030  // Now rebuild it
1031  for(const auto& hullPoint : convexHull.getConvexHull())
1032  {
1033 // fConvexHullList.emplace_back(std::get<0>(hullPoint),std::get<1>(hullPoint),std::get<2>(hullPoint));
1034 
1035  std::cout << "~~~ Convex hull Point: " << std::get<0>(hullPoint) << ", " << std::get<1>(hullPoint) << std::endl;
1036  }
1037  }
1038 
1039  return;
1040 }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
std::list< Point > PointList
The list of the projected points.
Definition: ConvexHull.h:32
dcel2d::Coords fConvexHullCenter
Definition: Voronoi.h:197
ConvexHull class definiton.
Definition: ConvexHull.h:25
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
dcel2d::PointList fConvexHullList
Definition: Voronoi.h:196
Eigen::Vector3f Coords
Definition: DCEL.h:36
VoronoiDiagram::PointPair voronoi2d::VoronoiDiagram::getExtremePoints ( ) const

Given an input Point, find the nearest edge.

Definition at line 1042 of file Voronoi.cxx.

References fConvexHullList.

Referenced by getConvexHull().

1043 {
1044  dcel2d::PointList::const_iterator nextPointItr = fConvexHullList.begin();
1045  dcel2d::PointList::const_iterator firstPointItr = nextPointItr++;
1046 
1047  float maxSeparation(0.);
1048 
1049  PointPair extremePoints(dcel2d::Point(0,0,NULL),dcel2d::Point(0,0,NULL));
1050 
1051  while(nextPointItr != fConvexHullList.end())
1052  {
1053  // Get a vector representing the edge from the first to the current next point
1054  dcel2d::Point firstPoint = *firstPointItr++;
1055  dcel2d::Point nextPoint = *nextPointItr;
1056  Eigen::Vector2f firstEdge(std::get<0>(*firstPointItr) - std::get<0>(firstPoint), std::get<1>(*firstPointItr) - std::get<1>(firstPoint));
1057 
1058  // normalize it
1059  firstEdge.normalize();
1060 
1061  dcel2d::PointList::const_iterator endPointItr = nextPointItr;
1062 
1063  while(++endPointItr != fConvexHullList.end())
1064  {
1065  dcel2d::Point endPoint = *endPointItr;
1066  Eigen::Vector2f nextEdge(std::get<0>(endPoint) - std::get<0>(nextPoint), std::get<1>(endPoint) - std::get<1>(nextPoint));
1067 
1068  // normalize it
1069  nextEdge.normalize();
1070 
1071  // Have we found the turnaround point?
1072  if (firstEdge.dot(nextEdge) < 0.)
1073  {
1074  Eigen::Vector2f separation(std::get<0>(nextPoint) - std::get<0>(firstPoint), std::get<1>(nextPoint) - std::get<1>(firstPoint));
1075  float separationDistance = separation.norm();
1076 
1077  if (separationDistance > maxSeparation)
1078  {
1079  extremePoints.first = firstPoint;
1080  extremePoints.second = nextPoint;
1081  maxSeparation = separationDistance;
1082  }
1083 
1084  // Regardless of thise being the maximum distance we have hit a turnaround point so
1085  // we need to break out of this loop
1086  break;
1087  }
1088 
1089  nextPointItr = endPointItr;
1090  nextPoint = endPoint;
1091  }
1092 
1093  // If we have hit the end of the convex hull without finding a turnaround point then we are not
1094  // going to find one so break out of the main loop
1095  if (endPointItr == fConvexHullList.end()) break;
1096 
1097  // Need to make sure we don't overrun the next point
1098  if (firstPointItr == nextPointItr) nextPointItr++;
1099  }
1100 
1101  return extremePoints;
1102 }
intermediate_table::const_iterator const_iterator
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:34
dcel2d::PointList fConvexHullList
Definition: Voronoi.h:196
std::pair< dcel2d::Point, dcel2d::Point > PointPair
Definition: Voronoi.h:36
const dcel2d::FaceList& voronoi2d::VoronoiDiagram::getFaceList ( ) const
inline

Recover the list of faces.

Definition at line 63 of file Voronoi.h.

References fFaceList.

63 {return fFaceList;}
dcel2d::FaceList & fFaceList
Definition: Voronoi.h:189
const dcel2d::VertexList& voronoi2d::VoronoiDiagram::getVertexList ( ) const
inline

Recover the list of vertices.

Definition at line 68 of file Voronoi.h.

References fVertexList.

68 {return fVertexList;}
dcel2d::VertexList & fVertexList
Definition: Voronoi.h:188
double voronoi2d::VoronoiDiagram::getVoronoiDiagramArea ( ) const
inline

recover the area of the convex hull

Definition at line 73 of file Voronoi.h.

References fVoronoiDiagramArea.

73 {return fVoronoiDiagramArea;}
void voronoi2d::VoronoiDiagram::handleCircleEvents ( BeachLine beachLine,
EventQueue eventQueue,
IEvent circleEvent 
)
private

There are two types of events in the queue, here we handle circle events.

Definition at line 474 of file Voronoi.cxx.

References voronoi2d::IEvent::circleCenter(), fHalfEdgeList, fVertexList, voronoi2d::BSTNode::getAssociated(), voronoi2d::IEvent::getBSTNode(), voronoi2d::BSTNode::getFace(), dcel2d::HalfEdge::getFace(), voronoi2d::BSTNode::getHalfEdge(), voronoi2d::BSTNode::getPredecessor(), voronoi2d::BSTNode::getSuccessor(), dcel2d::HalfEdge::getTwinHalfEdge(), makeLeftCircleEvent(), makeRightCircleEvent(), voronoi2d::BeachLine::removeLeaf(), dcel2d::HalfEdge::setFace(), dcel2d::Vertex::setHalfEdge(), voronoi2d::BSTNode::setHalfEdge(), dcel2d::HalfEdge::setLastHalfEdge(), dcel2d::HalfEdge::setNextHalfEdge(), dcel2d::HalfEdge::setTargetVertex(), dcel2d::HalfEdge::setTwinHalfEdge(), and voronoi2d::IEvent::xPos().

Referenced by buildVoronoiDiagram().

477 {
478  BSTNode* circleNode = circleEvent->getBSTNode();
479  BSTNode* arcNode = circleNode->getAssociated();
480 
481  // Recover the half edges for the breakpoints either side of the arc we're removing
482  dcel2d::HalfEdge* leftHalfEdge = arcNode->getPredecessor()->getHalfEdge();
483  dcel2d::HalfEdge* rightHalfEdge = arcNode->getSuccessor()->getHalfEdge();
484 
485  if (leftHalfEdge->getTwinHalfEdge()->getFace() != rightHalfEdge->getFace())
486  {
487  std::cout << ">>>> Face mismatch in circle handling, left face: " << leftHalfEdge->getFace() << ", left twin: " << leftHalfEdge->getTwinHalfEdge()->getFace() << ", right: " << rightHalfEdge->getFace() << ", right twin: " << rightHalfEdge->getTwinHalfEdge()->getFace() << ", arc face: " << arcNode->getFace() << std::endl;
488  }
489 
490  // Create the new vertex point
491  // Don't forget we need to swap coordinates when returning to the outside world
492  dcel2d::Coords vertexPos(circleEvent->circleCenter());
493 
494  fVertexList.push_back(dcel2d::Vertex(vertexPos,leftHalfEdge->getTwinHalfEdge()));
495 
496  dcel2d::Vertex* vertex = &fVertexList.back();
497 
498  // The edges we obtained "emanate" from the new vertex point,
499  // their twins will point to it
500  leftHalfEdge->getTwinHalfEdge()->setTargetVertex(vertex);
501  rightHalfEdge->getTwinHalfEdge()->setTargetVertex(vertex);
502 
503  // Remove the arc which will return the new breakpoint
504  // NOTE: invalidation of any possible circle events occurs in the call to this routine
505  BSTNode* newBreakPoint = beachLine.removeLeaf(arcNode);
506 
507  // Now create edges associated with the new break point
508  fHalfEdgeList.push_back(dcel2d::HalfEdge());
509 
510  dcel2d::HalfEdge* halfEdgeOne = &fHalfEdgeList.back();
511 
512  fHalfEdgeList.push_back(dcel2d::HalfEdge());
513 
514  dcel2d::HalfEdge* halfEdgeTwo = &fHalfEdgeList.back();
515 
516  // Associate the first edge with the new break point
517  newBreakPoint->setHalfEdge(halfEdgeTwo);
518 
519  // These are twins
520  halfEdgeOne->setTwinHalfEdge(halfEdgeTwo);
521  halfEdgeTwo->setTwinHalfEdge(halfEdgeOne);
522 
523  // The second points to the vertex we just made
524  halfEdgeTwo->setTargetVertex(vertex);
525  vertex->setHalfEdge(halfEdgeOne);
526 
527  // halfEdgeTwo points to the vertex, face to left, so pairs with
528  // the leftHalfEdge (leaving the vertex, face to left)
529  halfEdgeTwo->setNextHalfEdge(leftHalfEdge);
530  leftHalfEdge->setLastHalfEdge(halfEdgeTwo);
531  halfEdgeTwo->setFace(leftHalfEdge->getFace());
532 
533  // halfEdgeOne emanates from the vertex, face to left, so pairs with
534  // the twin of the rightHalfEdge (pointing to vertex, face to left)
535  halfEdgeOne->setLastHalfEdge(rightHalfEdge->getTwinHalfEdge());
536  rightHalfEdge->getTwinHalfEdge()->setNextHalfEdge(halfEdgeOne);
537  halfEdgeOne->setFace(rightHalfEdge->getTwinHalfEdge()->getFace());
538 
539  // Finally, the twin of the leftHalfEdge points to the vertex (face to left)
540  // and will pair with the rightHalfEdge emanating from the vertex
541  // In this case they should already be sharing the same face
542  rightHalfEdge->setLastHalfEdge(leftHalfEdge->getTwinHalfEdge());
543  leftHalfEdge->getTwinHalfEdge()->setNextHalfEdge(rightHalfEdge);
544 
545  // We'll try to make circle events with the remnant arcs so get the pointers now
546  // Note that we want the former left and right arcs to be the middle arcs in this case
547  BSTNode* leftArc = newBreakPoint->getPredecessor();
548  BSTNode* rightArc = newBreakPoint->getSuccessor();
549 
550  // Look for a new circle candidates
551  // In this case, we'd like the right arc to be the middle of the right circle,
552  // the left arc to be the middle of the left circle, hence the logic below
553  makeRightCircleEvent(eventQueue, leftArc, circleEvent->xPos());
554  makeLeftCircleEvent(eventQueue, rightArc, circleEvent->xPos());
555 
556  return;
557 }
dcel2d::HalfEdgeList & fHalfEdgeList
Definition: Voronoi.h:187
void setHalfEdge(HalfEdge *half)
Definition: DCEL.h:73
HalfEdge * getTwinHalfEdge() const
Definition: DCEL.h:159
Face * getFace() const
Definition: DCEL.h:158
void makeLeftCircleEvent(EventQueue &, BSTNode *, double)
Definition: Voronoi.cxx:559
void setLastHalfEdge(HalfEdge *last)
Definition: DCEL.h:167
void setTwinHalfEdge(HalfEdge *twin)
Definition: DCEL.h:165
dcel2d::VertexList & fVertexList
Definition: Voronoi.h:188
void makeRightCircleEvent(EventQueue &, BSTNode *, double)
Definition: Voronoi.cxx:608
void setFace(Face *face)
Definition: DCEL.h:164
void setNextHalfEdge(HalfEdge *next)
Definition: DCEL.h:166
Eigen::Vector3f Coords
Definition: DCEL.h:36
void setTargetVertex(Vertex *vertex)
Definition: DCEL.h:163
vertex reconstruction
void voronoi2d::VoronoiDiagram::handleSiteEvents ( BeachLine beachLine,
EventQueue eventQueue,
IEvent siteEvent 
)
private

There are two types of events in the queue, here we handle site events.

Definition at line 417 of file Voronoi.cxx.

References fFaceList, fHalfEdgeList, voronoi2d::IEvent::getCoords(), voronoi2d::BSTNode::getFace(), dcel2d::Face::getHalfEdge(), voronoi2d::IEvent::getPoint(), voronoi2d::BSTNode::getPredecessor(), voronoi2d::BSTNode::getSuccessor(), voronoi2d::BeachLine::insertNewLeaf(), makeLeftCircleEvent(), makeRightCircleEvent(), voronoi2d::BSTNode::setFace(), dcel2d::HalfEdge::setFace(), voronoi2d::BSTNode::setHalfEdge(), dcel2d::Face::setHalfEdge(), dcel2d::HalfEdge::setTwinHalfEdge(), and voronoi2d::IEvent::xPos().

Referenced by buildVoronoiDiagram().

420 {
421  // Insert the new site event into the beach line and recover the leaf for the
422  // new arc in the beach line
423  // NOTE: invalidation of any possible circle events occurs in the call to this routine
424  BSTNode* newLeaf = beachLine.insertNewLeaf(siteEvent);
425 
426  // Create a new vertex for each site event
427  fFaceList.emplace_back(dcel2d::Face(NULL,siteEvent->getCoords(),std::get<2>(siteEvent->getPoint())));
428 
429  dcel2d::Face* face = &fFaceList.back();
430 
431  newLeaf->setFace(face);
432 
433  // If we are the first site added to the BST then we are done
434  if (!(newLeaf->getPredecessor() || newLeaf->getSuccessor())) return;
435 
436  // So, now we deal with creating the edges that will be mapped out by the two breakpoints as they
437  // move with the beachline. Note that this is a single edge pair (edge and its twin) between the
438  // two breakpoints that have been created.
439  fHalfEdgeList.push_back(dcel2d::HalfEdge());
440 
441  dcel2d::HalfEdge* halfEdge = &fHalfEdgeList.back();
442 
443  fHalfEdgeList.push_back(dcel2d::HalfEdge());
444 
445  dcel2d::HalfEdge* halfTwin = &fHalfEdgeList.back();
446 
447  // Each half edge is the twin of the other
448  halfEdge->setTwinHalfEdge(halfTwin);
449  halfTwin->setTwinHalfEdge(halfEdge);
450 
451  // Point the breakpoint to the new edge
452  newLeaf->getPredecessor()->setHalfEdge(halfEdge);
453  newLeaf->getSuccessor()->setHalfEdge(halfTwin);
454 
455  // The second half edge corresponds to our face for the left breakpoint...
456  face->setHalfEdge(halfTwin);
457  halfTwin->setFace(face);
458 
459  // Update for the left leaf
460  BSTNode* leftLeaf = newLeaf->getPredecessor()->getPredecessor();
461 
462  // This needs to be checked since the first face will not have an edge set to it
463  if (!leftLeaf->getFace()->getHalfEdge()) leftLeaf->getFace()->setHalfEdge(halfEdge);
464 
465  halfEdge->setFace(leftLeaf->getFace());
466 
467  // Try to make circle events either side of this leaf
468  makeLeftCircleEvent( eventQueue, newLeaf, siteEvent->xPos());
469  makeRightCircleEvent(eventQueue, newLeaf, siteEvent->xPos());
470 
471  return;
472 }
dcel2d::HalfEdgeList & fHalfEdgeList
Definition: Voronoi.h:187
dcel2d::FaceList & fFaceList
Definition: Voronoi.h:189
void makeLeftCircleEvent(EventQueue &, BSTNode *, double)
Definition: Voronoi.cxx:559
void setTwinHalfEdge(HalfEdge *twin)
Definition: DCEL.h:165
void makeRightCircleEvent(EventQueue &, BSTNode *, double)
Definition: Voronoi.cxx:608
void setFace(Face *face)
Definition: DCEL.h:164
void setHalfEdge(HalfEdge *half)
Definition: DCEL.h:110
bool voronoi2d::VoronoiDiagram::isInsideConvexHull ( const dcel2d::Vertex vertex) const
private

Is a vertex inside the convex hull - meant to be a fast check.

Definition at line 1104 of file Voronoi.cxx.

References fConvexHullList, dcel2d::Vertex::getCoords(), and isLeft().

Referenced by terminateInfiniteEdges().

1105 {
1106  bool insideHull(true);
1107  dcel2d::Point vertexPos(vertex.getCoords()[0],vertex.getCoords()[1],NULL);
1108 
1109  // Now check to see if the vertex is inside the convex hull
1111  dcel2d::Point firstPoint = *hullItr++;
1112 
1113  // We assume here the convex hull is stored in a CCW order
1114  while(hullItr != fConvexHullList.end())
1115  {
1116  dcel2d::Point secondPoint = *hullItr++;
1117 
1118  // CCW order means we check to see if this point lies to left of line from first to second point
1119  // in order for the point to be "inside" the convex hull
1120  // Note that we are looking to reject points that are outside...
1121  if (!isLeft(firstPoint,secondPoint,vertexPos))
1122  {
1123  insideHull = false;
1124  break;
1125  }
1126 
1127  firstPoint = secondPoint;
1128  }
1129 
1130  return insideHull;
1131 }
intermediate_table::const_iterator const_iterator
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:34
dcel2d::PointList fConvexHullList
Definition: Voronoi.h:196
const Coords & getCoords() const
Definition: DCEL.h:61
bool isLeft(const dcel2d::Point &p0, const dcel2d::Point &p1, const dcel2d::Point &pCheck) const
Determines whether a point is to the left, right or on line specifiec by p0 and p1.
Definition: Voronoi.cxx:84
bool voronoi2d::VoronoiDiagram::isLeft ( const dcel2d::Point p0,
const dcel2d::Point p1,
const dcel2d::Point pCheck 
) const
private

Determines whether a point is to the left, right or on line specifiec by p0 and p1.

Definition at line 84 of file Voronoi.cxx.

References crossProduct().

Referenced by findNearestEdge(), isInsideConvexHull(), and isOutsideConvexHull().

85 {
86  // Use the cross product to determine if the check point lies to the left, on or right
87  // of the line defined by points p0 and p1
88  return crossProduct(p0, p1, pCheck) > 0;
89 }
double crossProduct(const dcel2d::Point &p0, const dcel2d::Point &p1, const dcel2d::Point &p2) const
Gets the cross product of line from p0 to p1 and p0 to p2.
Definition: Voronoi.cxx:93
bool voronoi2d::VoronoiDiagram::isOutsideConvexHull ( const dcel2d::Vertex vertex,
dcel2d::PointList::const_iterator  firstHullPointItr,
dcel2d::Coords intersection,
double &  distToConvexHull 
) const
private

is vertex outside the convex hull and if so return some useful information

Definition at line 1133 of file Voronoi.cxx.

References fConvexHullList, dcel2d::Vertex::getCoords(), isLeft(), and max.

1137 {
1138  bool outsideHull(false);
1139  dcel2d::Point vertexPos(vertex.getCoords()[0],vertex.getCoords()[1],NULL);
1140 
1141  distToConvexHull = std::numeric_limits<double>::max();
1142  firstHullPointItr = fConvexHullList.begin();
1143 
1144  // Now check to see if the vertex is inside the convex hull
1145  dcel2d::PointList::const_iterator hullItr = firstHullPointItr;
1146  dcel2d::PointList::const_iterator firstPointItr = hullItr++;
1147 
1148  while(hullItr != fConvexHullList.end())
1149  {
1150  dcel2d::PointList::const_iterator secondPointItr = hullItr++;
1151 
1152  // Dereference some stuff
1153  double xPrevToPoint = (std::get<0>(vertexPos) - std::get<0>(*firstPointItr));
1154  double yPrevToPoint = (std::get<1>(vertexPos) - std::get<1>(*firstPointItr));
1155  double xPrevToCur = (std::get<0>(*secondPointItr) - std::get<0>(*firstPointItr));
1156  double yPrevToCur = (std::get<1>(*secondPointItr) - std::get<1>(*firstPointItr));
1157  double edgeLength = std::sqrt(xPrevToCur * xPrevToCur + yPrevToCur * yPrevToCur);
1158 
1159  // Find projection onto convex hull edge
1160  double projection = ((xPrevToPoint * xPrevToCur) + (yPrevToPoint * yPrevToCur)) / edgeLength;
1161 
1162  // DOCA point
1163  dcel2d::Point docaPoint(std::get<0>(*firstPointItr) + projection * xPrevToCur / edgeLength,
1164  std::get<1>(*firstPointItr) + projection * yPrevToCur / edgeLength, 0);
1165 
1166  if (projection > edgeLength) docaPoint = *secondPointItr;
1167  if (projection < 0) docaPoint = *firstPointItr;
1168 
1169  double xDocaDist = std::get<0>(vertexPos) - std::get<0>(docaPoint);
1170  double yDocaDist = std::get<1>(vertexPos) - std::get<1>(docaPoint);
1171  double docaDist = xDocaDist * xDocaDist + yDocaDist * yDocaDist;
1172 
1173  if (docaDist < distToConvexHull)
1174  {
1175  firstHullPointItr = firstPointItr;
1176  intersection = dcel2d::Coords(std::get<0>(docaPoint),std::get<1>(docaPoint),0.);
1177  distToConvexHull = docaDist;
1178  }
1179 
1180  // Check to see if this point is outside the convex hull
1181  if (isLeft(*firstPointItr,*secondPointItr,vertexPos)) outsideHull = true;
1182 
1183  firstPointItr = secondPointItr;
1184  }
1185 
1186  return outsideHull;
1187 }
Int_t max
Definition: plot.C:27
intermediate_table::const_iterator const_iterator
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:34
dcel2d::PointList fConvexHullList
Definition: Voronoi.h:196
Eigen::Vector3f Coords
Definition: DCEL.h:36
const Coords & getCoords() const
Definition: DCEL.h:61
bool isLeft(const dcel2d::Point &p0, const dcel2d::Point &p1, const dcel2d::Point &pCheck) const
Determines whether a point is to the left, right or on line specifiec by p0 and p1.
Definition: Voronoi.cxx:84
IEvent * voronoi2d::VoronoiDiagram::makeCircleEvent ( BSTNode arc1,
BSTNode arc2,
BSTNode arc3,
double  beachLinePos 
)
private

There are two types of events in the queue, here we handle circle events.

Definition at line 655 of file Voronoi.cxx.

References computeCircleCenter(), e, fCircleEventList, voronoi2d::IEvent::getCoords(), voronoi2d::BSTNode::getEvent(), and radius.

Referenced by makeLeftCircleEvent(), and makeRightCircleEvent().

656 {
657  // It might be that we don't create a new circle
658  IEvent* circle = 0;
659 
660  // First step is to calculate the center and radius of the circle determined by the three input site events
661  dcel2d::Coords center;
662  double radius;
663  double deltaR;
664 
665  IEvent* p1 = arc1->getEvent();
666  IEvent* p2 = arc2->getEvent();
667  IEvent* p3 = arc3->getEvent();
668 
669  // Compute the center of the circle. Note that this will also automagically check that breakpoints are
670  // converging so that this is a circle we want
671  if (computeCircleCenter( p1->getCoords(), p2->getCoords(), p3->getCoords(), center, radius, deltaR))
672  {
673  double circleBottomX = center[0] - radius;
674 
675  // Now check if the bottom of this circle lies below the beach line
676  if (beachLinePos >= circleBottomX - 10. * deltaR)
677  {
678  // Making a circle event!
679  dcel2d::Point circleBottom(circleBottomX, center[1], NULL);
680 
681  fCircleEventList.emplace_back(circleBottom, center);
682 
683  circle = &fCircleEventList.back();
684  }
685  else if (circleBottomX - beachLinePos < 1.e-4)
686  std::cout << "==> Circle close, beachLine: " << beachLinePos << ", circleBottomX: " << circleBottomX << ", deltaR: " << deltaR << ", d: " << circleBottomX - beachLinePos << std::endl;
687  }
688 
689  return circle;
690 }
bool computeCircleCenter(const dcel2d::Coords &, const dcel2d::Coords &, const dcel2d::Coords &, dcel2d::Coords &, double &, double &) const
Definition: Voronoi.cxx:692
CircleEventList fCircleEventList
Definition: Voronoi.h:193
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:34
Double_t radius
Eigen::Vector3f Coords
Definition: DCEL.h:36
Float_t e
Definition: plot.C:34
void voronoi2d::VoronoiDiagram::makeLeftCircleEvent ( EventQueue eventQueue,
BSTNode leaf,
double  beachLine 
)
private

Definition at line 559 of file Voronoi.cxx.

References fCircleNodeList, voronoi2d::BSTNode::getAssociated(), voronoi2d::BSTNode::getEvent(), voronoi2d::IEvent::getPoint(), voronoi2d::BSTNode::getPredecessor(), makeCircleEvent(), voronoi2d::BSTNode::setAssociated(), and voronoi2d::IEvent::setInvalid().

Referenced by handleCircleEvents(), and handleSiteEvents().

562 {
563 
564  // Check status of triplet of site events to the left of this new leaf
565  if (leaf->getPredecessor())
566  {
567  BSTNode* midLeaf = leaf->getPredecessor()->getPredecessor();
568 
569  if (midLeaf && midLeaf->getPredecessor())
570  {
571  BSTNode* edgeLeaf = midLeaf->getPredecessor()->getPredecessor();
572 
573  // edge leaves must be different
574  if (leaf->getEvent()->getPoint() == edgeLeaf->getEvent()->getPoint()) return;
575 
576  if (edgeLeaf)
577  {
578  IEvent* circleEvent = makeCircleEvent(edgeLeaf, midLeaf, leaf, beachLine);
579 
580  // Did we succeed in making a circle event?
581  if (circleEvent)
582  {
583  // Add to the circle node list
584  fCircleNodeList.emplace_back(circleEvent);
585 
586  BSTNode* circleNode = &fCircleNodeList.back();
587 
588  // If there was an associated circle event to this node, invalidate it
589  if (midLeaf->getAssociated())
590  {
591  midLeaf->getAssociated()->setAssociated(NULL);
592  midLeaf->getAssociated()->getEvent()->setInvalid();
593  }
594 
595  // Now reset to point at the new one
596  midLeaf->setAssociated(circleNode);
597  circleNode->setAssociated(midLeaf);
598 
599  eventQueue.push(circleEvent);
600  }
601  }
602  }
603  }
604 
605  return;
606 }
BSTNodeList fCircleNodeList
Definition: Voronoi.h:194
IEvent * makeCircleEvent(BSTNode *, BSTNode *, BSTNode *, double)
There are two types of events in the queue, here we handle circle events.
Definition: Voronoi.cxx:655
void voronoi2d::VoronoiDiagram::makeRightCircleEvent ( EventQueue eventQueue,
BSTNode leaf,
double  beachLine 
)
private

Definition at line 608 of file Voronoi.cxx.

References fCircleNodeList, voronoi2d::BSTNode::getAssociated(), voronoi2d::BSTNode::getEvent(), voronoi2d::IEvent::getPoint(), voronoi2d::BSTNode::getSuccessor(), makeCircleEvent(), voronoi2d::BSTNode::setAssociated(), and voronoi2d::IEvent::setInvalid().

Referenced by handleCircleEvents(), and handleSiteEvents().

611 {
612  // Check status of triplet of site events to the left of this new leaf
613  if (leaf->getSuccessor())
614  {
615  BSTNode* midLeaf = leaf->getSuccessor()->getSuccessor();
616 
617  if (midLeaf && midLeaf->getSuccessor())
618  {
619  BSTNode* edgeLeaf = midLeaf->getSuccessor()->getSuccessor();
620 
621  if (leaf->getEvent()->getPoint() == edgeLeaf->getEvent()->getPoint()) return;
622 
623  if (edgeLeaf)
624  {
625  IEvent* circleEvent = makeCircleEvent(leaf, midLeaf, edgeLeaf, beachLine);
626 
627  // Did we succeed in making a circle event?
628  if (circleEvent)
629  {
630  // Add to the circle node list
631  fCircleNodeList.emplace_back(circleEvent);
632 
633  BSTNode* circleNode = &fCircleNodeList.back();
634 
635  // If there was an associated circle event to this node, invalidate it
636  if (midLeaf->getAssociated())
637  {
638  midLeaf->getAssociated()->setAssociated(NULL);
639  midLeaf->getAssociated()->getEvent()->setInvalid();
640  }
641 
642  // Now reset to point at the new one
643  midLeaf->setAssociated(circleNode);
644  circleNode->setAssociated(midLeaf);
645 
646  eventQueue.push(circleEvent);
647  }
648  }
649  }
650  }
651 
652  return;
653 }
BSTNodeList fCircleNodeList
Definition: Voronoi.h:194
IEvent * makeCircleEvent(BSTNode *, BSTNode *, BSTNode *, double)
There are two types of events in the queue, here we handle circle events.
Definition: Voronoi.cxx:655
void voronoi2d::VoronoiDiagram::mergeDegenerateVertices ( )
private

merge degenerate vertices (found by zero length edges)

Definition at line 1189 of file Voronoi.cxx.

References fHalfEdgeList, dcel2d::Vertex::getCoords(), dcel2d::HalfEdge::getTargetVertex(), and dcel2d::HalfEdge::getTwinHalfEdge().

Referenced by terminateInfiniteEdges().

1190 {
1192 
1193  while(edgeItr != fHalfEdgeList.end())
1194  {
1195  dcel2d::HalfEdge* halfEdge = &(*edgeItr);
1196  dcel2d::HalfEdge* twinEdge = halfEdge->getTwinHalfEdge();
1197 
1198  // Make sure we are not looking at an infinite edge
1199  if (halfEdge->getTargetVertex() && twinEdge->getTargetVertex())
1200  {
1201  dcel2d::Coords vtxPosDiff = halfEdge->getTargetVertex()->getCoords() - twinEdge->getTargetVertex()->getCoords();
1202 
1203  if (vtxPosDiff.norm() < 1.e-3)
1204  {
1205  std::cout << "***>> found a degenerate vertex! " << halfEdge->getTargetVertex()->getCoords()[0] << "/" << halfEdge->getTargetVertex()->getCoords()[1] << ", d: " << vtxPosDiff.norm() << std::endl;
1206  edgeItr++;
1207  }
1208  else edgeItr++;
1209  }
1210  else edgeItr++;
1211  }
1212 
1213  return;
1214 }
dcel2d::HalfEdgeList & fHalfEdgeList
Definition: Voronoi.h:187
intermediate_table::iterator iterator
HalfEdge * getTwinHalfEdge() const
Definition: DCEL.h:159
Vertex * getTargetVertex() const
Definition: DCEL.h:157
Eigen::Vector3f Coords
Definition: DCEL.h:36
const Coords & getCoords() const
Definition: DCEL.h:61
void voronoi2d::VoronoiDiagram::terminateInfiniteEdges ( BeachLine beachLine,
double  beachLinePos 
)
private

this aims to process remaining elements in the beachline after the event queue has been depleted

Definition at line 831 of file Voronoi.cxx.

References voronoi2d::EventUtilities::computeArcVal(), voronoi2d::EventUtilities::computeBreak(), ComputeFaceArea(), fConvexHullList, fUtilities, fVertexList, voronoi2d::BSTNode::getAssociated(), getConvexHull(), voronoi2d::IEvent::getCoords(), dcel2d::Vertex::getCoords(), voronoi2d::BSTNode::getEvent(), voronoi2d::BSTNode::getFace(), dcel2d::HalfEdge::getFace(), voronoi2d::BSTNode::getHalfEdge(), dcel2d::HalfEdge::getLastHalfEdge(), voronoi2d::BSTNode::getLeftChild(), dcel2d::HalfEdge::getNextHalfEdge(), voronoi2d::BSTNode::getPredecessor(), voronoi2d::BSTNode::getRightChild(), voronoi2d::BSTNode::getSuccessor(), dcel2d::HalfEdge::getTargetVertex(), voronoi2d::BeachLine::getTopNode(), dcel2d::HalfEdge::getTwinHalfEdge(), voronoi2d::IEvent::isCircle(), isInsideConvexHull(), voronoi2d::IEvent::isValid(), mergeDegenerateVertices(), dcel2d::HalfEdge::setTargetVertex(), voronoi2d::IEvent::xPos(), and voronoi2d::IEvent::yPos().

Referenced by buildVoronoiDiagram().

832 {
833  // Need to complete processing of the beachline, the remaning leaves represent the site points with "infinite"
834  // edges which we need to terminate at our bounding box.
835  // For now let's just do step one which is to find the convex hull
836  const BSTNode* node = beachLine.getTopNode();
837 
838  // Do the convex hull independently but before terminating edges
839  getConvexHull(node);
840 
841  if (node)
842  {
843  // Run down the left branches to find the left most leaf
844  while(node->getLeftChild()) node = node->getLeftChild();
845 
846  // Things to keep track of...
847  int nodeCount(0);
848  int nBadBreaks(0);
849  double lastBreakPoint = std::numeric_limits<double>::lowest();
850 
851  // Now we are going to traverse across the beach line and process the remaining leaves
852  // This will identify the infinite edges and find the convex hull
853  while(node)
854  {
855  std::cout << "--------------------------------------------------------------------------------------------" << std::endl;
856 
857  // Do we have a leaf?
858  if (!node->getLeftChild() && !node->getRightChild())
859  {
860  std::cout << "node " << nodeCount << " has leaf: " << node << ", face: " << node->getFace() << ", pos: " << node->getEvent()->getCoords()[0] << "/" << node->getEvent()->getCoords()[1] << std::endl;
861  }
862  // Otherwise it is a break point
863  else
864  {
865  RootsPair roots;
866  double breakPoint = fUtilities.computeBreak(beachLinePos, node->getPredecessor()->getEvent(), node->getSuccessor()->getEvent(), roots);
867  double leftArcVal = fUtilities.computeArcVal(beachLinePos, breakPoint, node->getPredecessor()->getEvent());
868  double rightArcVal = fUtilities.computeArcVal(beachLinePos, breakPoint, node->getSuccessor()->getEvent());
869 
870  dcel2d::HalfEdge* halfEdge = node->getHalfEdge();
871  dcel2d::HalfEdge* halfTwin = halfEdge->getTwinHalfEdge();
872  dcel2d::Vertex* vertex = halfEdge->getTargetVertex();
873 
874  std::cout << "node " << nodeCount << " has break: " << breakPoint << ", beachPos: " << beachLinePos << " (last: " << lastBreakPoint << "), leftArcVal: " << leftArcVal << ", rightArcVal: " << rightArcVal << std::endl;
875  if (lastBreakPoint > breakPoint) std::cout << " ***** lastBreakPoint larger than current break point *****" << std::endl;
876  std::cout << " left arc x/y: " << node->getPredecessor()->getEvent()->xPos() << "/" << node->getPredecessor()->getEvent()->yPos() << ", right arc x/y: " << node->getSuccessor()->getEvent()->xPos() << "/" << node->getSuccessor()->getEvent()->yPos() << std::endl;
877  std::cout << " halfEdge: " << halfEdge << ", target vtx: " << vertex << ", face: " << halfEdge->getFace() << ", twin: " << halfTwin << ", twin tgt: " << halfTwin->getTargetVertex() << ", twin face: " << halfTwin->getFace() << ", next: " << halfEdge->getNextHalfEdge() << ", last: " << halfEdge->getLastHalfEdge() << std::endl;
878 
879  const dcel2d::Coords& leftLeafPos = node->getPredecessor()->getEvent()->getCoords();
880  const dcel2d::Coords& rightLeafPos = node->getSuccessor()->getEvent()->getCoords();
881  Eigen::Vector3f lrPosDiff = rightLeafPos - leftLeafPos;
882  Eigen::Vector3f lrPosSum = rightLeafPos + leftLeafPos;
883 
884  lrPosDiff.normalize();
885  lrPosSum *= 0.5;
886 
887  dcel2d::Coords leafPos = node->getSuccessor()->getEvent()->getCoords();
888  dcel2d::Coords leafDir = lrPosDiff;
889 
890  if (vertex)
891  {
892  dcel2d::Coords breakDir(leafDir[1],-leafDir[0],0.);
893  dcel2d::Coords breakPos = lrPosSum;
894  dcel2d::Coords vertexPos = vertex->getCoords();
895 
896  // Get the arclength from the break position to the intersection of the two lines
897  dcel2d::Coords breakToLeafPos = leafPos - vertexPos;
898  double arcLenToLine = breakToLeafPos.dot(breakDir);
899 
900  // Now get position
901  dcel2d::Coords breakVertexPos = vertexPos + arcLenToLine * breakDir;
902 
903  std::cout << " halfEdge position: " << breakPos[0] << "/" << breakPos[1] << ", vertex: " << vertexPos[0] << "/" << vertexPos[1] << ", end: " << breakVertexPos[0] << "/" << breakVertexPos[1] << ", dir: " << breakDir[0] << "/" << breakDir[1] << ", arclen: " << arcLenToLine << std::endl;
904 
905  // At this point we need to create a vertex and then terminate the half edges here...
906  fVertexList.push_back(dcel2d::Vertex(breakVertexPos,halfEdge));
907 
908  dcel2d::Vertex* breakVertex = &fVertexList.back();
909 
910  halfTwin->setTargetVertex(breakVertex);
911  }
912  else
913  {
914  std::cout << "****** null vertex!!! Skipping to next node... *********" << std::endl;
915  }
916 
917  if (lastBreakPoint > breakPoint) nBadBreaks++;
918 
919  lastBreakPoint = breakPoint;
920  }
921 
922  if (node->getAssociated())
923  {
924  BSTNode* associated = node->getAssociated();
925  IEvent* iEvent = associated->getEvent();
926 
927  std::cout << " -> associated circle: " << iEvent->isCircle() << ", is valid: " << iEvent->isValid() << std::endl;
928  }
929 
930  node = node->getSuccessor();
931  nodeCount++;
932  }
933 
934  // Now that we have the convex hull, loop over vertices to see if they are inside or outside of the convex hull
935  dcel2d::VertexList::iterator curVertexItr = fVertexList.begin();
936 
937  size_t nVerticesInitial = fVertexList.size();
938 
939  // Loop over all vertices to begin with
940  while(curVertexItr != fVertexList.end())
941  {
942  // Dereference vertex
943  dcel2d::Vertex& vertex = *curVertexItr;
944 
945  bool outsideHull = !isInsideConvexHull(vertex);
946 
947  // Do we need to drop this vertex?
948  if (outsideHull)
949  {
950  curVertexItr = fVertexList.erase(curVertexItr);
951  }
952  else curVertexItr++;
953  }
954 
956  ComputeFaceArea();
957 
958  std::cout << "Loop over beachline done, saved " << fConvexHullList.size() << " arcs, encountered " << nBadBreaks << " bad break points" << std::endl;
959  std::cout << "-- started with " << nVerticesInitial << " vertices, found " << fVertexList.size() << " inside convex hull" << std::endl;
960  }
961 
962  return;
963 }
HalfEdge * getLastHalfEdge() const
Definition: DCEL.h:161
HalfEdge * getNextHalfEdge() const
Definition: DCEL.h:160
intermediate_table::iterator iterator
HalfEdge * getTwinHalfEdge() const
Definition: DCEL.h:159
Face * getFace() const
Definition: DCEL.h:158
bool isInsideConvexHull(const dcel2d::Vertex &) const
Is a vertex inside the convex hull - meant to be a fast check.
Definition: Voronoi.cxx:1104
double computeArcVal(const double, const double, const IEvent *) const
double ComputeFaceArea()
Compute the area of the faces.
Definition: Voronoi.cxx:1217
EventUtilities fUtilities
Definition: Voronoi.h:206
dcel2d::VertexList & fVertexList
Definition: Voronoi.h:188
void mergeDegenerateVertices()
merge degenerate vertices (found by zero length edges)
Definition: Voronoi.cxx:1189
dcel2d::PointList fConvexHullList
Definition: Voronoi.h:196
double computeBreak(const double, const IEvent *, const IEvent *, RootsPair &) const
const dcel2d::PointList & getConvexHull() const
recover the point list representing the convex hull
Definition: Voronoi.h:78
std::pair< double, double > RootsPair
Vertex * getTargetVertex() const
Definition: DCEL.h:157
Eigen::Vector3f Coords
Definition: DCEL.h:36
const Coords & getCoords() const
Definition: DCEL.h:61
void setTargetVertex(Vertex *vertex)
Definition: DCEL.h:163
vertex reconstruction

Member Data Documentation

CircleEventList voronoi2d::VoronoiDiagram::fCircleEventList
private
BSTNodeList voronoi2d::VoronoiDiagram::fCircleNodeList
private
dcel2d::Coords voronoi2d::VoronoiDiagram::fConvexHullCenter
private

Definition at line 197 of file Voronoi.h.

Referenced by getConvexHull().

dcel2d::PointList voronoi2d::VoronoiDiagram::fConvexHullList
private
dcel2d::FaceList& voronoi2d::VoronoiDiagram::fFaceList
private
int voronoi2d::VoronoiDiagram::fNumBadCircles
mutableprivate

Definition at line 203 of file Voronoi.h.

Referenced by buildVoronoiDiagram(), buildVoronoiDiagramBoost(), and computeCircleCenter().

dcel2d::PointList voronoi2d::VoronoiDiagram::fPointList
private
SiteEventList voronoi2d::VoronoiDiagram::fSiteEventList
private

Definition at line 192 of file Voronoi.h.

Referenced by buildVoronoiDiagram(), buildVoronoiDiagramBoost(), and VoronoiDiagram().

EventUtilities voronoi2d::VoronoiDiagram::fUtilities
private

Definition at line 206 of file Voronoi.h.

Referenced by terminateInfiniteEdges().

double voronoi2d::VoronoiDiagram::fVoronoiDiagramArea
private

Definition at line 204 of file Voronoi.h.

Referenced by getVoronoiDiagramArea(), and VoronoiDiagram().

double voronoi2d::VoronoiDiagram::fXMax
private

Definition at line 200 of file Voronoi.h.

Referenced by buildVoronoiDiagram(), and findBoundingBox().

double voronoi2d::VoronoiDiagram::fXMin
private

Definition at line 199 of file Voronoi.h.

Referenced by buildVoronoiDiagram(), and findBoundingBox().

double voronoi2d::VoronoiDiagram::fYMax
private

Definition at line 202 of file Voronoi.h.

Referenced by buildVoronoiDiagram(), and findBoundingBox().

double voronoi2d::VoronoiDiagram::fYMin
private

Definition at line 201 of file Voronoi.h.

Referenced by buildVoronoiDiagram(), and findBoundingBox().


The documentation for this class was generated from the following files: