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...

	~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::FaceList &	getFaceList () const
	Recover the list of faces. More...

const dcel2d::VertexList &	getVertexList () const
	Recover the list of vertices. More...

double	getVoronoiDiagramArea () const
	recover the area of the convex hull More...

const dcel2d::PointList &	getConvexHull () 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)

IEvent *	makeCircleEvent (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 &) const

bool	computeCircleCenter3 (const dcel2d::Coords &, const dcel2d::Coords &, const dcel2d::Coords &, dcel2d::Coords &, 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::HalfEdgeList &	fHalfEdgeList

dcel2d::VertexList &	fVertexList

dcel2d::FaceList &	fFaceList

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 30 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 167 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 164 of file Voronoi.h.

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

private

Definition at line 166 of file Voronoi.h.

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

private

Definition at line 93 of file Voronoi.h.

using voronoi2d::VoronoiDiagram::MinMaxPointPair = std::pair<PointPair, PointPair>

Definition at line 33 of file Voronoi.h.

using voronoi2d::VoronoiDiagram::PointPair = std::pair<dcel2d::Point, dcel2d::Point>

Definition at line 32 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 47 of file Voronoi.cxx.

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

     : fHalfEdgeList(halfEdgeList)
     , fVertexList(vertexList)
     , fFaceList(faceList)
     , fXMin(0.)
     , fXMax(0.)
     , fYMin(0.)
     , fYMax(0.)
     , fVoronoiDiagramArea(0.)
   {
     fHalfEdgeList.clear();
     fVertexList.clear();
     fFaceList.clear();
     fPointList.clear();
     fSiteEventList.clear();
     fCircleEventList.clear();
     fCircleNodeList.clear();
     fConvexHullList.clear();
 
     // And the area
     fVoronoiDiagramArea = Area();
   }

voronoi2d::VoronoiDiagram::~VoronoiDiagram ( )

Destructor.

Definition at line 74 of file Voronoi.cxx.

74 {}

Member Function Documentation

double voronoi2d::VoronoiDiagram::Area ( ) const

private

Computes the area of the enclosed convext hull.

Definition at line 104 of file Voronoi.cxx.

References crossProduct(), and fPointList.

Referenced by VoronoiDiagram().

   {
     double area(0.);
 
     // Compute the area by taking advantage of
     // 1) the ability to decompose a convex hull into triangles,
     // 2) the ability to use the cross product to calculate the area
     // So, the technique is to pick a point (for technical reasons we use 0,0)
     // and then sum the signed area of triangles formed from this point to two adjecent
     // vertices on the convex hull.
     dcel2d::Point center(0., 0., NULL);
     dcel2d::Point lastPoint = fPointList.front();
 
     for (const auto& point : fPointList) {
       if (point != lastPoint) area += 0.5 * crossProduct(center, lastPoint, point);
 
       lastPoint = point;
     }
 
     return area;
   }

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 355 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().

   {
     dcel2d::HalfEdge* halfEdge = NULL;
     dcel2d::HalfEdge* twinEdge = NULL;
 
     if (boostEdgeToEdgeMap.find(edge) != boostEdgeToEdgeMap.end())
       halfEdge = boostEdgeToEdgeMap.at(edge);
     else {
       fHalfEdgeList.emplace_back();
 
       halfEdge = &fHalfEdgeList.back();
 
       boostEdgeToEdgeMap[edge] = halfEdge;
     }
 
     if (boostEdgeToEdgeMap.find(twin) != boostEdgeToEdgeMap.end())
       twinEdge = boostEdgeToEdgeMap.at(twin);
     else {
       fHalfEdgeList.emplace_back();
 
       twinEdge = &fHalfEdgeList.back();
 
       boostEdgeToEdgeMap[twin] = twinEdge;
     }
 
     // Do the primary half edge first
     const boost::polygon::voronoi_vertex<double>* boostVertex = edge->vertex1();
     dcel2d::Vertex* vertex = NULL;
 
     // note we can have a null vertex (infinite edge)
     if (boostVertex) {
       if (boostVertexToVertexMap.find(boostVertex) == boostVertexToVertexMap.end()) {
         dcel2d::Coords coords(boostVertex->y(), boostVertex->x(), 0.);
 
         fVertexList.emplace_back(coords, halfEdge);
 
         vertex = &fVertexList.back();
 
         boostVertexToVertexMap[boostVertex] = vertex;
       }
       else
         vertex = boostVertexToVertexMap.at(boostVertex);
     }
 
     const boost::polygon::voronoi_cell<double>* boostCell = edge->cell();
     dcel2d::Face* face = NULL;
 
     if (boostCellToFaceMap.find(boostCell) == boostCellToFaceMap.end()) {
       dcel2d::PointList::const_iterator pointItr = pointList.begin();
       int pointIdx = boostCell->source_index();
 
       std::advance(pointItr, pointIdx);
 
       const dcel2d::Point& point = *pointItr;
       dcel2d::Coords coords(std::get<0>(point), std::get<1>(point), 0.);
 
       fFaceList.emplace_back(halfEdge, coords, std::get<2>(point));
 
       face = &fFaceList.back();
 
       boostCellToFaceMap[boostCell] = face;
     }
 
     halfEdge->setTargetVertex(vertex);
     halfEdge->setFace(face);
     halfEdge->setTwinHalfEdge(twinEdge);
 
     // For the prev/next half edges we can have two cases, so check:
     if (boostEdgeToEdgeMap.find(edge->next()) != boostEdgeToEdgeMap.end()) {
       dcel2d::HalfEdge* nextEdge = boostEdgeToEdgeMap.at(edge->next());
 
       halfEdge->setNextHalfEdge(nextEdge);
       nextEdge->setLastHalfEdge(halfEdge);
     }
 
     if (boostEdgeToEdgeMap.find(edge->prev()) != boostEdgeToEdgeMap.end()) {
       dcel2d::HalfEdge* lastEdge = boostEdgeToEdgeMap.at(edge->prev());
 
       halfEdge->setLastHalfEdge(lastEdge);
       lastEdge->setNextHalfEdge(halfEdge);
     }
 
     return;
   }

void voronoi2d::VoronoiDiagram::buildVoronoiDiagram ( const dcel2d::PointList & pointList )

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

Parameters

PointList The 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().

   {
     // Insure all the local data structures have been cleared
     fHalfEdgeList.clear();
     fVertexList.clear();
     fFaceList.clear();
     fPointList.clear();
     fSiteEventList.clear();
     fCircleEventList.clear();
     fCircleNodeList.clear();
     fNumBadCircles = 0;
 
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "==> # input points: " << pointList.size() << std::endl;
 
     // Define the priority queue to contain our events
     EventQueue eventQueue(compareSiteEventPtrs);
 
     // Now populate the event queue with site events
     for (const auto& point : pointList) {
       fSiteEventList.emplace_back(point);
       IEvent* iEvent = &fSiteEventList.back();
       eventQueue.push(iEvent);
     }
 
     // Declare the beachline which will contain the BSTNode objects for site events
     BeachLine beachLine;
 
     // Finally, we need a container for our circle event BSTNodes
     BSTNodeList circleNodeList;
 
     // Now process the queue
     while (!eventQueue.empty()) {
       IEvent* event = eventQueue.top();
 
       eventQueue.pop();
 
       // If a site or circle event then handle appropriately
       if (event->isSite())
         handleSiteEvents(beachLine, eventQueue, event);
       else if (event->isValid())
         handleCircleEvents(beachLine, eventQueue, event);
     }
 
     std::cout << "*******> # input points: " << pointList.size()
               << ", remaining leaves: " << beachLine.countLeaves() << ", "
               << beachLine.traverseBeach() << ", # bad circles: " << fNumBadCircles << std::endl;
     std::cout << "           Faces: " << fFaceList.size() << ", Vertices: " << fVertexList.size()
               << ", # half edges: " << fHalfEdgeList.size() << std::endl;
 
     // Get the bounding box
     findBoundingBox(fVertexList);
 
     std::cout << "           Range min/maxes, x: " << fXMin << ", " << fXMax << ", y: " << fYMin
               << ", " << fYMax << std::endl;
 
     // Terminate the infinite edges
     terminateInfiniteEdges(beachLine, std::get<0>(pointList.front()));
 
     // Look for open faces
     int nOpenFaces(0);
 
     std::map<int, int> edgeCountMap;
     std::map<int, int> openCountMap;
 
     for (const auto& face : fFaceList) {
       int nEdges(1);
       bool closed(false);
       const dcel2d::HalfEdge* startEdge = face.getHalfEdge();
 
       // Count forwards
       for (const dcel2d::HalfEdge* halfEdge = startEdge->getNextHalfEdge(); halfEdge && !closed;) {
         if (halfEdge->getFace() != &face) {
           std::cout << "  ===> halfEdge does not match face: " << halfEdge
                     << ", face: " << halfEdge->getFace() << ", base: " << &face << std::endl;
         }
 
         if (halfEdge == startEdge) {
           closed = true;
           break;
         }
         nEdges++;
         halfEdge = halfEdge->getNextHalfEdge();
       }
 
       // Count backwards. Note if face is closed then nothing to do here
       if (!closed) {
         for (const dcel2d::HalfEdge* halfEdge = startEdge->getLastHalfEdge();
              halfEdge && !closed;) {
           if (halfEdge->getFace() != &face) {
             std::cout << "  ===> halfEdge does not match face: " << halfEdge
                       << ", face: " << halfEdge->getFace() << ", base: " << &face << std::endl;
           }
 
           if (halfEdge == startEdge) {
             closed = true;
             break;
           }
           nEdges++;
           halfEdge = halfEdge->getLastHalfEdge();
         }
       }
 
       if (!closed) {
         nOpenFaces++;
         openCountMap[nEdges]++;
       }
 
       edgeCountMap[nEdges]++;
     }
 
     std::cout << "==> Found " << nOpenFaces << " open faces from total of " << fFaceList.size()
               << std::endl;
     for (const auto& edgeCount : edgeCountMap)
       std::cout << "    -> all edges,  # edges: " << edgeCount.first
                 << ", count: " << edgeCount.second << std::endl;
     for (const auto& edgeCount : openCountMap)
       std::cout << "    -> open edges, # edges: " << edgeCount.first
                 << ", count: " << edgeCount.second << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
 
     // Clear internal containers that are no longer useful
     fSiteEventList.clear();
     fCircleEventList.clear();
     fCircleNodeList.clear();
 
     return;
   }

void voronoi2d::VoronoiDiagram::buildVoronoiDiagramBoost ( const dcel2d::PointList & pointList )

Definition at line 282 of file Voronoi.cxx.

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

   {
     // Insure all the local data structures have been cleared
     fHalfEdgeList.clear();
     fVertexList.clear();
     fFaceList.clear();
     fPointList.clear();
     fSiteEventList.clear();
     fCircleEventList.clear();
     fCircleNodeList.clear();
     fNumBadCircles = 0;
 
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "==> # input points: " << pointList.size() << std::endl;
 
     // Construct out voronoi diagram
     boost::polygon::voronoi_diagram<double> vd;
     boost::polygon::construct_voronoi(pointList.begin(), pointList.end(), &vd);
 
     // Make maps for translating from boost to me (or we can rewrite our code in boost... for now maps)
     using BoostEdgeToEdgeMap =
       std::map<const boost::polygon::voronoi_edge<double>*, dcel2d::HalfEdge*>;
     using BoostVertexToVertexMap =
       std::map<const boost::polygon::voronoi_vertex<double>*, dcel2d::Vertex*>;
     using BoostCellToFaceMap = std::map<const boost::polygon::voronoi_cell<double>*, dcel2d::Face*>;
 
     BoostEdgeToEdgeMap boostEdgeToEdgeMap;
     BoostVertexToVertexMap boostVertexToVertexMap;
     BoostCellToFaceMap boostCellToFaceMap;
 
     // Loop over the edges
     for (const auto& edge : vd.edges()) {
       const boost::polygon::voronoi_edge<double>* twin = edge.twin();
 
       boostTranslation(
         pointList, &edge, twin, boostEdgeToEdgeMap, boostVertexToVertexMap, boostCellToFaceMap);
       boostTranslation(
         pointList, twin, &edge, boostEdgeToEdgeMap, boostVertexToVertexMap, boostCellToFaceMap);
     }
 
     //std::cout << "==> Found " << nOpenFaces << " open faces from total of " << fFaceList.size() << std::endl;
     //for(const auto& edgeCount : edgeCountMap) std::cout << "    -> all edges,  # edges: " << edgeCount.first << ", count: " << edgeCount.second << std::endl;
     //for(const auto& edgeCount : openCountMap) std::cout << "    -> open edges, # edges: " << edgeCount.first << ", count: " << edgeCount.second << std::endl;
     std::cout << "==> Boost returns " << fHalfEdgeList.size() << " edges, " << fFaceList.size()
               << " faces, " << fVertexList.size() << " vertices " << std::endl;
     std::cout << "                  " << vd.edges().size() << " edges, " << vd.cells().size()
               << " faces, " << vd.vertices().size() << " vertices" << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
     std::cout << "*********************************************************************************"
                  "*********************************"
               << std::endl;
 
     // Clear internal containers that are no longer useful
     fSiteEventList.clear();
     fCircleEventList.clear();
     fCircleNodeList.clear();
 
     return;
   }

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 714 of file Voronoi.cxx.

References e, fNumBadCircles, x2, x3, y2, and y3.

Referenced by makeCircleEvent().

   {
     // The method is to translate the three points to a system where the first point is at the origin. Then we
     // are looking for a circle that passes through the origin and the two remaining (translated) points. In
     // this mode the circle radius is the distance from the origin to the circle center which simplifies the
     // calculation.
     double xCoord = p1[0];
     double yCoord = p1[1];
 
     double x2 = p2[0] - xCoord;
     double y2 = p2[1] - yCoord;
     double x3 = p3[0] - xCoord;
     double y3 = p3[1] - yCoord;
 
     double det = x2 * y3 - x3 * y2;
 
     // Points are colinear so cannot make a circle
     // Or, if negative then the midpoint is "left" of line between first and third points meaning the
     // circle curvature is the "wrong" way for making a circle event
     if (det <= 0.) // std::numeric_limits<double>::epsilon())
     //if (!(std::abs(det) > 0.)) // std::numeric_limits<double>::epsilon())
     {
       if (det > -std::numeric_limits<double>::epsilon())
         std::cout << "      --->Circle failure, det: " << det << ", mid x: " << p2[0]
                   << ", y: " << p2[1] << ",   l x: " << p1[0] << ", y: " << p1[1]
                   << ",   r x: " << p3[0] << ", y: " << p3[1] << std::endl;
 
       return false;
     }
 
     double p2sqr = x2 * x2 + y2 * y2;
     double p3sqr = x3 * x3 + y3 * y3;
 
     double cxpr = 0.5 * (y3 * p2sqr - y2 * p3sqr) / det;
     double cypr = 0.5 * (x2 * p3sqr - x3 * p2sqr) / det;
 
     radius = std::sqrt(cxpr * cxpr + cypr * cypr);
     center[0] = cxpr + xCoord;
     center[1] = cypr + yCoord;
     center[2] = 0.;
 
     // So... roundoff error can cause degenerate circles to get missed
     // We need to attack that problem by calculating the radius all possible
     // ways and then take the largest...
     dcel2d::Coords p1Rad = p1 - center;
     dcel2d::Coords p2Rad = p2 - center;
     dcel2d::Coords p3Rad = p3 - center;
 
     std::vector<float> radSqrVec;
 
     radSqrVec.push_back(p1Rad.norm());
     radSqrVec.push_back(p2Rad.norm());
     radSqrVec.push_back(p3Rad.norm());
 
     double maxRadius = *std::max_element(radSqrVec.begin(), radSqrVec.end());
 
     delta = std::max(5.e-7, maxRadius - radius);
 
     if (radius > 1000.) {
       //        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;
       fNumBadCircles++;
     }
 
     return true;
   }

bool voronoi2d::VoronoiDiagram::computeCircleCenter2	(	const dcel2d::Coords &	p1,
		const dcel2d::Coords &	p2,
		const dcel2d::Coords &	p3,
		dcel2d::Coords &	center,
		double &	radius
	)		const

private

Definition at line 785 of file Voronoi.cxx.

References util::abs().

   {
     // Compute the circle center as the intersection of the two perpendicular bisectors of rays between the points
     double slope12 = (p2[1] - p1[1]) / (p2[0] - p1[0]);
     double slope32 = (p3[1] - p2[1]) / (p3[0] - p2[0]);
 
     if (std::abs(slope12 - slope32) <= std::numeric_limits<double>::epsilon()) {
       std::cout << "       >>>> Matching slopes! points: (" << p1[0] << "," << p1[1] << "), ("
                 << p2[0] << "," << p2[1] << "), (" << p3[0] << "," << p3[1] << ")" << std::endl;
 
       return false;
     }
 
     center[0] = (slope12 * slope32 * (p3[1] - p1[1]) + slope12 * (p2[0] + p3[0]) -
                  slope32 * (p1[0] + p2[0])) /
                 (2. * (slope12 - slope32));
     center[1] = 0.5 * (p1[1] + p2[1]) - (center[0] - 0.5 * (p1[0] + p2[0])) / slope12;
     center[2] = 0.;
     radius = std::sqrt(std::pow((p2[0] - center[0]), 2) + std::pow((p2[1] - center[1]), 2));
 
     if (radius > 100.) {
       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;
     }
 
     return true;
   }

bool voronoi2d::VoronoiDiagram::computeCircleCenter3	(	const dcel2d::Coords &	p1,
		const dcel2d::Coords &	p2,
		const dcel2d::Coords &	p3,
		dcel2d::Coords &	center,
		double &	radius
	)		const

private

Definition at line 818 of file Voronoi.cxx.

References util::abs().

   {
     // Yet another bisector method to calculate the circle center...
     double temp = p2[0] * p2[0] + p2[1] * p2[1];
     double p1p2 = (p1[0] * p1[0] + p1[1] * p1[1] - temp) / 2.0;
     double p2p3 = (temp - p3[0] * p3[0] - p3[1] * p3[1]) / 2.0;
     double det = (p1[0] - p2[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p2[1]);
 
     if (std::abs(det) < 10. * std::numeric_limits<double>::epsilon()) {
       std::cout << "       >>>> Determinant zero! points: (" << p1[0] << "," << p1[1] << "), ("
                 << p2[0] << "," << p2[1] << "), (" << p3[0] << "," << p3[1] << ")" << std::endl;
 
       return false;
     }
 
     det = 1. / det;
 
     center[0] = (p1p2 * (p2[1] - p3[1]) - p2p3 * (p1[1] - p2[1])) * det;
     center[1] = (p2p3 * (p1[0] - p2[0]) - p1p2 * (p2[0] - p3[0])) * det;
     center[2] = 0.;
 
     radius = std::sqrt(std::pow((p1[0] - center[0]), 2) + std::pow((p1[1] - center[1]), 2));
 
     if (radius > 100.) {
       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;
     }
 
     return true;
   }

double voronoi2d::VoronoiDiagram::ComputeFaceArea ( )

private

Compute the area of the faces.

Definition at line 1270 of file Voronoi.cxx.

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

Referenced by terminateInfiniteEdges().

   {
     // Compute the area by taking advantage of
     // 1) the ability to decompose a convex hull into triangles,
     // 2) the ability to use the cross product to calculate the area
     // So the idea is to loop through all the faces and then follow the edges
     // around the face to compute the area of the face.
     // Note that a special case are the "infinite faces" which lie on the
     // convex hull of the event. Skip these for now...
 
     double totalArea(0.);
     int nNonInfiniteFaces(0);
     double smallestArea(std::numeric_limits<double>::max());
     double largestArea(0.);
 
     std::vector<std::pair<double, const dcel2d::Face*>> areaFaceVec;
 
     areaFaceVec.reserve(fFaceList.size());
 
     for (auto& face : fFaceList) {
       //        const dcel2d::Coords&   faceCoords = face.getCoords();
       const dcel2d::HalfEdge* halfEdge = face.getHalfEdge();
       double faceArea(0.);
       int numEdges(0);
       bool doNext = true;
 
       dcel2d::Coords faceCenter(0., 0., 0.);
 
       while (doNext) {
         if (halfEdge->getTargetVertex()) faceCenter += halfEdge->getTargetVertex()->getCoords();
 
         numEdges++;
 
         halfEdge = halfEdge->getNextHalfEdge();
 
         if (!halfEdge) {
           faceArea = std::numeric_limits<double>::max();
           doNext = false;
         }
 
         if (halfEdge == face.getHalfEdge()) doNext = false;
       }
 
       faceCenter /= numEdges;
 
       halfEdge = face.getHalfEdge();
       doNext = true;
 
       while (doNext) {
         const dcel2d::HalfEdge* twinEdge = halfEdge->getTwinHalfEdge();
 
         if (!halfEdge->getTargetVertex() || !twinEdge->getTargetVertex()) {
           faceArea = std::numeric_limits<double>::max();
           break;
         }
 
         // Recover the two vertex points
         const dcel2d::Coords& p1 = halfEdge->getTargetVertex()->getCoords();
         const dcel2d::Coords& p2 = twinEdge->getTargetVertex()->getCoords();
 
         // Define a quick 2D cross product here since it will used quite a bit!
         double dp1p0X = p1[0] - faceCenter[0];
         double dp1p0Y = p1[1] - faceCenter[1];
         double dp2p0X = p2[0] - faceCenter[0];
         double dp2p0Y = p2[1] - faceCenter[1];
 
         //faceArea += dp1p0X * dp2p0Y - dp1p0Y * dp2p0X;
         double crossProduct = dp1p0X * dp2p0Y - dp1p0Y * dp2p0X;
 
         faceArea += crossProduct;
 
         if (crossProduct < 0.) {
           dcel2d::Coords edgeVec = p1 - p2;
           std::cout << "--- negative cross: " << crossProduct << ", edgeLen: " << edgeVec.norm()
                     << ", x/y: " << edgeVec[0] << "/" << edgeVec[1] << std::endl;
         }
 
         //            numEdges++;
 
         halfEdge = halfEdge->getNextHalfEdge();
 
         if (!halfEdge) {
           faceArea = std::numeric_limits<double>::max();
           break;
         }
 
         if (halfEdge == face.getHalfEdge()) doNext = false;
       }
 
       areaFaceVec.emplace_back(faceArea, &face);
 
       if (faceArea < std::numeric_limits<double>::max() && faceArea > 0.) {
         nNonInfiniteFaces++;
         totalArea += faceArea;
         smallestArea = std::min(faceArea, smallestArea);
         largestArea = std::max(faceArea, largestArea);
       }
 
       if (faceArea < 1.e-4)
         std::cout << "---> face area <= 0: " << faceArea << ", with " << numEdges << " edges"
                   << std::endl;
 
       face.setFaceArea(faceArea);
     }
 
     // Calculate a truncated mean...
     std::sort(areaFaceVec.begin(), areaFaceVec.end(), [](const auto& left, const auto& right) {
       return left.first < right.first;
     });
 
     std::vector<std::pair<double, const dcel2d::Face*>>::iterator firstItr = std::find_if(
       areaFaceVec.begin(), areaFaceVec.end(), [](const auto& val) { return val.first > 0.; });
     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());
       });
 
     size_t nToKeep = 0.8 * std::distance(firstItr, lastItr);
 
     std::cout << ">>>>> nToKeep: " << nToKeep
               << ", last non infinite entry: " << std::distance(areaFaceVec.begin(), lastItr)
               << std::endl;
 
     double totalTruncArea = std::accumulate(
       firstItr, firstItr + nToKeep, 0., [](auto sum, const auto& val) { return sum + val.first; });
     double aveTruncArea = totalTruncArea / double(nToKeep);
 
     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;
     else
       std::cout << ">>>>> there are no non infinite faces" << std::endl;
 
     return totalArea;
   }

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 89 of file Voronoi.cxx.

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

   {
     // Define a quick 2D cross product here since it will used quite a bit!
     double deltaX = std::get<0>(p1) - std::get<0>(p0);
     double deltaY = std::get<1>(p1) - std::get<1>(p0);
     double dCheckX = std::get<0>(p2) - std::get<0>(p0);
     double dCheckY = std::get<1>(p2) - std::get<1>(p0);
 
     return ((deltaX * dCheckY) - (deltaY * dCheckX));
   }

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 1408 of file Voronoi.cxx.

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

Referenced by buildVoronoiDiagram().

   {
     // Find extremes in x to start
     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];
         });
 
     fXMin = minMaxItrX.first->getCoords()[0];
     fXMax = minMaxItrX.second->getCoords()[0];
 
     // To get the extremes in y we need to make a pass through the list
     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];
         });
 
     fYMin = minMaxItrY.first->getCoords()[1];
     fYMax = minMaxItrY.second->getCoords()[1];
 
     return;
   }

double voronoi2d::VoronoiDiagram::findNearestDistance ( const dcel2d::Point & point ) const

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

Definition at line 1494 of file Voronoi.cxx.

References findNearestEdge().

Referenced by getConvexHull().

   {
     double closestDistance;
 
     findNearestEdge(point, closestDistance);
 
     return closestDistance;
   }

VoronoiDiagram::PointPair voronoi2d::VoronoiDiagram::findNearestEdge	(	const dcel2d::Point &	point,
		double &	closestDistance
	)		const

Given an input Point, find the nearest edge.

Definition at line 1433 of file Voronoi.cxx.

References fPointList, and isLeft().

Referenced by findNearestDistance(), and getConvexHull().

   {
     // The idea is to find the nearest edge of the convex hull, defined by
     // two adjacent vertices of the hull, to the input point.
     // As near as I can tell, the best way to do this is to apply brute force...
     // Idea will be to iterate over pairs of points
     dcel2d::PointList::const_iterator curPointItr = fPointList.begin();
     dcel2d::Point prevPoint = *curPointItr++;
     dcel2d::Point curPoint = *curPointItr;
 
     // Set up the winner
     PointPair closestEdge(prevPoint, curPoint);
 
     closestDistance = std::numeric_limits<double>::max();
 
     // curPointItr is meant to point to the second point
     while (curPointItr != fPointList.end()) {
       if (curPoint != prevPoint) {
         // Dereference some stuff
         double xPrevToPoint = (std::get<0>(point) - std::get<0>(prevPoint));
         double yPrevToPoint = (std::get<1>(point) - std::get<1>(prevPoint));
         double xPrevToCur = (std::get<0>(curPoint) - std::get<0>(prevPoint));
         double yPrevToCur = (std::get<1>(curPoint) - std::get<1>(prevPoint));
         double edgeLength = std::sqrt(xPrevToCur * xPrevToCur + yPrevToCur * yPrevToCur);
 
         // Find projection onto convex hull edge
         double projection =
           ((xPrevToPoint * xPrevToCur) + (yPrevToPoint * yPrevToCur)) / edgeLength;
 
         // DOCA point
         dcel2d::Point docaPoint(std::get<0>(prevPoint) + projection * xPrevToCur / edgeLength,
                                 std::get<1>(prevPoint) + projection * yPrevToCur / edgeLength,
                                 0);
 
         if (projection > edgeLength) docaPoint = curPoint;
         if (projection < 0) docaPoint = prevPoint;
 
         double xDocaDist = std::get<0>(point) - std::get<0>(docaPoint);
         double yDocaDist = std::get<1>(point) - std::get<1>(docaPoint);
         double docaDist = xDocaDist * xDocaDist + yDocaDist * yDocaDist;
 
         if (docaDist < closestDistance) {
           closestEdge = PointPair(prevPoint, curPoint);
           closestDistance = docaDist;
         }
       }
 
       prevPoint = curPoint;
       curPoint = *curPointItr++;
     }
 
     closestDistance = std::sqrt(closestDistance);
 
     // Convention is convex hull vertices sorted in counter clockwise fashion so if the point
     // lays to the left of the nearest edge then it must be an interior point
     if (isLeft(closestEdge.first, closestEdge.second, point)) closestDistance = -closestDistance;
 
     return closestEdge;
   }

const dcel2d::PointList& voronoi2d::VoronoiDiagram::getConvexHull ( ) const

inline

recover the point list representing the convex hull

Definition at line 74 of file Voronoi.h.

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

Referenced by terminateInfiniteEdges().

74 { return fConvexHullList; }

voronoi2d::VoronoiDiagram::fConvexHullList

dcel2d::PointList fConvexHullList

Definition: Voronoi.h:214

void voronoi2d::VoronoiDiagram::getConvexHull ( const BSTNode * topNode )

private

this function recovers the convex hull

Definition at line 1009 of file Voronoi.cxx.

References util::abs(), 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().

   {
     // Assume the input node is the top of the binary search tree and represents
     // the beach line at the end of the sweep algorithm
     // The convex hull will then be the leaf elements along the beach line.
     // Note that this algorithm will give you the convex hull in a clockwise manner
     if (topNode) {
       const BSTNode* node = topNode;
 
       // Initialize the center
       fConvexHullCenter = dcel2d::Coords(0., 0., 0.);
 
       // Run down the left branches to find the right most leaf
       // We do it this way to make sure we get a CCW convex hull
       while (node->getRightChild())
         node = node->getRightChild();
 
       // Include a check on convexity...
       Eigen::Vector3f prevVec(0., 0., 0.);
       dcel2d::Coords lastPoint(0., 0., 0.);
 
       // Now we are going to traverse across the beach line and process the remaining leaves
       // This will identify the convex hull
       while (node) {
         if (!node->getLeftChild() && !node->getRightChild()) {
           // Add point to the convex hull list
           fConvexHullList.emplace_back(node->getEvent()->getPoint());
 
           // Mark the face as on the convex hull
           node->getFace()->setOnConvexHull();
 
           dcel2d::Coords nextPoint = node->getEvent()->getCoords();
           Eigen::Vector3f curVec = nextPoint - lastPoint;
 
           curVec.normalize();
 
           double dotProd = prevVec.dot(curVec);
 
           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;
 
           prevVec = curVec;
           lastPoint = nextPoint;
         }
 
         node = node->getPredecessor();
       }
 
       // Annoyingly, our algorithm does not contain only the convex hull points and so we need to skim out the renegades...
       lar_cluster3d::ConvexHull::PointList localList;
 
       for (const auto& edgePoint : fConvexHullList)
         localList.emplace_back(
           std::get<0>(edgePoint), std::get<1>(edgePoint), std::get<2>(edgePoint));
 
       // Sort the point vec by increasing x, then increase y
       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);
       });
 
       // Why do I need to do this?
       lar_cluster3d::ConvexHull convexHull(localList);
 
       // Clear the convex hull list...
       //        fConvexHullList.clear();
 
       std::cout << "~~~>> there are " << convexHull.getConvexHull().size()
                 << " convex hull points and " << fConvexHullList.size() << " infinite cells"
                 << std::endl;
 
       // Now rebuild it
       for (const auto& hullPoint : convexHull.getConvexHull()) {
         //            fConvexHullList.emplace_back(std::get<0>(hullPoint),std::get<1>(hullPoint),std::get<2>(hullPoint));
 
         std::cout << "~~~ Convex hull Point: " << std::get<0>(hullPoint) << ", "
                   << std::get<1>(hullPoint) << std::endl;
       }
     }
 
     return;
   }

VoronoiDiagram::PointPair voronoi2d::VoronoiDiagram::getExtremePoints ( ) const

Given an input Point, find the nearest edge.

Definition at line 1097 of file Voronoi.cxx.

References fConvexHullList.

Referenced by getConvexHull().

   {
     dcel2d::PointList::const_iterator nextPointItr = fConvexHullList.begin();
     dcel2d::PointList::const_iterator firstPointItr = nextPointItr++;
 
     float maxSeparation(0.);
 
     PointPair extremePoints(dcel2d::Point(0, 0, NULL), dcel2d::Point(0, 0, NULL));
 
     while (nextPointItr != fConvexHullList.end()) {
       // Get a vector representing the edge from the first to the current next point
       dcel2d::Point firstPoint = *firstPointItr++;
       dcel2d::Point nextPoint = *nextPointItr;
       Eigen::Vector2f firstEdge(std::get<0>(*firstPointItr) - std::get<0>(firstPoint),
                                 std::get<1>(*firstPointItr) - std::get<1>(firstPoint));
 
       // normalize it
       firstEdge.normalize();
 
       dcel2d::PointList::const_iterator endPointItr = nextPointItr;
 
       while (++endPointItr != fConvexHullList.end()) {
         dcel2d::Point endPoint = *endPointItr;
         Eigen::Vector2f nextEdge(std::get<0>(endPoint) - std::get<0>(nextPoint),
                                  std::get<1>(endPoint) - std::get<1>(nextPoint));
 
         // normalize it
         nextEdge.normalize();
 
         // Have we found the turnaround point?
         if (firstEdge.dot(nextEdge) < 0.) {
           Eigen::Vector2f separation(std::get<0>(nextPoint) - std::get<0>(firstPoint),
                                      std::get<1>(nextPoint) - std::get<1>(firstPoint));
           float separationDistance = separation.norm();
 
           if (separationDistance > maxSeparation) {
             extremePoints.first = firstPoint;
             extremePoints.second = nextPoint;
             maxSeparation = separationDistance;
           }
 
           // Regardless of thise being the maximum distance we have hit a turnaround point so
           // we need to break out of this loop
           break;
         }
 
         nextPointItr = endPointItr;
         nextPoint = endPoint;
       }
 
       // If we have hit the end of the convex hull without finding a turnaround point then we are not
       // going to find one so break out of the main loop
       if (endPointItr == fConvexHullList.end()) break;
 
       // Need to make sure we don't overrun the next point
       if (firstPointItr == nextPointItr) nextPointItr++;
     }
 
     return extremePoints;
   }

const dcel2d::FaceList& voronoi2d::VoronoiDiagram::getFaceList ( ) const

inline

Recover the list of faces.

Definition at line 59 of file Voronoi.h.

References fFaceList.

59 { return fFaceList; }

voronoi2d::VoronoiDiagram::fFaceList

dcel2d::FaceList & fFaceList

Definition: Voronoi.h:207

const dcel2d::VertexList& voronoi2d::VoronoiDiagram::getVertexList ( ) const

inline

Recover the list of vertices.

Definition at line 64 of file Voronoi.h.

References fVertexList.

64 { return fVertexList; }

voronoi2d::VoronoiDiagram::fVertexList

dcel2d::VertexList & fVertexList

Definition: Voronoi.h:206

double voronoi2d::VoronoiDiagram::getVoronoiDiagramArea ( ) const

inline

recover the area of the convex hull

Definition at line 69 of file Voronoi.h.

References fVoronoiDiagramArea.

69 { return fVoronoiDiagramArea; }

voronoi2d::VoronoiDiagram::fVoronoiDiagramArea

double fVoronoiDiagramArea

Definition: Voronoi.h:222

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 503 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().

   {
     BSTNode* circleNode = circleEvent->getBSTNode();
     BSTNode* arcNode = circleNode->getAssociated();
 
     // Recover the half edges for the breakpoints either side of the arc we're removing
     dcel2d::HalfEdge* leftHalfEdge = arcNode->getPredecessor()->getHalfEdge();
     dcel2d::HalfEdge* rightHalfEdge = arcNode->getSuccessor()->getHalfEdge();
 
     if (leftHalfEdge->getTwinHalfEdge()->getFace() != rightHalfEdge->getFace()) {
       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;
     }
 
     // Create the new vertex point
     // Don't forget we need to swap coordinates when returning to the outside world
     dcel2d::Coords vertexPos(circleEvent->circleCenter());
 
     fVertexList.push_back(dcel2d::Vertex(vertexPos, leftHalfEdge->getTwinHalfEdge()));
 
     dcel2d::Vertex* vertex = &fVertexList.back();
 
     // The edges we obtained "emanate" from the new vertex point,
     // their twins will point to it
     leftHalfEdge->getTwinHalfEdge()->setTargetVertex(vertex);
     rightHalfEdge->getTwinHalfEdge()->setTargetVertex(vertex);
 
     // Remove the arc which will return the new breakpoint
     // NOTE: invalidation of any possible circle events occurs in the call to this routine
     BSTNode* newBreakPoint = beachLine.removeLeaf(arcNode);
 
     // Now create edges associated with the new break point
     fHalfEdgeList.push_back(dcel2d::HalfEdge());
 
     dcel2d::HalfEdge* halfEdgeOne = &fHalfEdgeList.back();
 
     fHalfEdgeList.push_back(dcel2d::HalfEdge());
 
     dcel2d::HalfEdge* halfEdgeTwo = &fHalfEdgeList.back();
 
     // Associate the first edge with the new break point
     newBreakPoint->setHalfEdge(halfEdgeTwo);
 
     // These are twins
     halfEdgeOne->setTwinHalfEdge(halfEdgeTwo);
     halfEdgeTwo->setTwinHalfEdge(halfEdgeOne);
 
     // The second points to the vertex we just made
     halfEdgeTwo->setTargetVertex(vertex);
     vertex->setHalfEdge(halfEdgeOne);
 
     // halfEdgeTwo points to the vertex, face to left, so pairs with
     // the leftHalfEdge (leaving the vertex, face to left)
     halfEdgeTwo->setNextHalfEdge(leftHalfEdge);
     leftHalfEdge->setLastHalfEdge(halfEdgeTwo);
     halfEdgeTwo->setFace(leftHalfEdge->getFace());
 
     // halfEdgeOne emanates from the vertex, face to left, so pairs with
     // the twin of the rightHalfEdge (pointing to vertex, face to left)
     halfEdgeOne->setLastHalfEdge(rightHalfEdge->getTwinHalfEdge());
     rightHalfEdge->getTwinHalfEdge()->setNextHalfEdge(halfEdgeOne);
     halfEdgeOne->setFace(rightHalfEdge->getTwinHalfEdge()->getFace());
 
     // Finally, the twin of the leftHalfEdge points to the vertex (face to left)
     // and will pair with the rightHalfEdge emanating from the vertex
     // In this case they should already be sharing the same face
     rightHalfEdge->setLastHalfEdge(leftHalfEdge->getTwinHalfEdge());
     leftHalfEdge->getTwinHalfEdge()->setNextHalfEdge(rightHalfEdge);
 
     // We'll try to make circle events with the remnant arcs so get the pointers now
     // Note that we want the former left and right arcs to be the middle arcs in this case
     BSTNode* leftArc = newBreakPoint->getPredecessor();
     BSTNode* rightArc = newBreakPoint->getSuccessor();
 
     // Look for a new circle candidates
     // In this case, we'd like the right arc to be the middle of the right circle,
     // the left arc to be the middle of the left circle, hence the logic below
     makeRightCircleEvent(eventQueue, leftArc, circleEvent->xPos());
     makeLeftCircleEvent(eventQueue, rightArc, circleEvent->xPos());
 
     return;
   }

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 445 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().

   {
     // Insert the new site event into the beach line and recover the leaf for the
     // new arc in the beach line
     // NOTE: invalidation of any possible circle events occurs in the call to this routine
     BSTNode* newLeaf = beachLine.insertNewLeaf(siteEvent);
 
     // Create a new vertex for each site event
     fFaceList.emplace_back(
       dcel2d::Face(NULL, siteEvent->getCoords(), std::get<2>(siteEvent->getPoint())));
 
     dcel2d::Face* face = &fFaceList.back();
 
     newLeaf->setFace(face);
 
     // If we are the first site added to the BST then we are done
     if (!(newLeaf->getPredecessor() || newLeaf->getSuccessor())) return;
 
     // So, now we deal with creating the edges that will be mapped out by the two breakpoints as they
     // move with the beachline. Note that this is a single edge pair (edge and its twin) between the
     // two breakpoints that have been created.
     fHalfEdgeList.push_back(dcel2d::HalfEdge());
 
     dcel2d::HalfEdge* halfEdge = &fHalfEdgeList.back();
 
     fHalfEdgeList.push_back(dcel2d::HalfEdge());
 
     dcel2d::HalfEdge* halfTwin = &fHalfEdgeList.back();
 
     // Each half edge is the twin of the other
     halfEdge->setTwinHalfEdge(halfTwin);
     halfTwin->setTwinHalfEdge(halfEdge);
 
     // Point the breakpoint to the new edge
     newLeaf->getPredecessor()->setHalfEdge(halfEdge);
     newLeaf->getSuccessor()->setHalfEdge(halfTwin);
 
     // The second half edge corresponds to our face for the left breakpoint...
     face->setHalfEdge(halfTwin);
     halfTwin->setFace(face);
 
     // Update for the left leaf
     BSTNode* leftLeaf = newLeaf->getPredecessor()->getPredecessor();
 
     // This needs to be checked since the first face will not have an edge set to it
     if (!leftLeaf->getFace()->getHalfEdge()) leftLeaf->getFace()->setHalfEdge(halfEdge);
 
     halfEdge->setFace(leftLeaf->getFace());
 
     // Try to make circle events either side of this leaf
     makeLeftCircleEvent(eventQueue, newLeaf, siteEvent->xPos());
     makeRightCircleEvent(eventQueue, newLeaf, siteEvent->xPos());
 
     return;
   }

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 1158 of file Voronoi.cxx.

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

Referenced by terminateInfiniteEdges().

   {
     bool insideHull(true);
     dcel2d::Point vertexPos(vertex.getCoords()[0], vertex.getCoords()[1], NULL);
 
     // Now check to see if the vertex is inside the convex hull
     dcel2d::PointList::const_iterator hullItr = fConvexHullList.begin();
     dcel2d::Point firstPoint = *hullItr++;
 
     // We assume here the convex hull is stored in a CCW order
     while (hullItr != fConvexHullList.end()) {
       dcel2d::Point secondPoint = *hullItr++;
 
       // CCW order means we check to see if this point lies to left of line from first to second point
       // in order for the point to be "inside" the convex hull
       // Note that we are looking to reject points that are outside...
       if (!isLeft(firstPoint, secondPoint, vertexPos)) {
         insideHull = false;
         break;
       }
 
       firstPoint = secondPoint;
     }
 
     return insideHull;
   }

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 78 of file Voronoi.cxx.

References crossProduct().

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

   {
     // Use the cross product to determine if the check point lies to the left, on or right
     // of the line defined by points p0 and p1
     return crossProduct(p0, p1, pCheck) > 0;
   }

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 1185 of file Voronoi.cxx.

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

   {
     bool outsideHull(false);
     dcel2d::Point vertexPos(vertex.getCoords()[0], vertex.getCoords()[1], NULL);
 
     distToConvexHull = std::numeric_limits<double>::max();
     firstHullPointItr = fConvexHullList.begin();
 
     // Now check to see if the vertex is inside the convex hull
     dcel2d::PointList::const_iterator hullItr = firstHullPointItr;
     dcel2d::PointList::const_iterator firstPointItr = hullItr++;
 
     while (hullItr != fConvexHullList.end()) {
       dcel2d::PointList::const_iterator secondPointItr = hullItr++;
 
       // Dereference some stuff
       double xPrevToPoint = (std::get<0>(vertexPos) - std::get<0>(*firstPointItr));
       double yPrevToPoint = (std::get<1>(vertexPos) - std::get<1>(*firstPointItr));
       double xPrevToCur = (std::get<0>(*secondPointItr) - std::get<0>(*firstPointItr));
       double yPrevToCur = (std::get<1>(*secondPointItr) - std::get<1>(*firstPointItr));
       double edgeLength = std::sqrt(xPrevToCur * xPrevToCur + yPrevToCur * yPrevToCur);
 
       // Find projection onto convex hull edge
       double projection = ((xPrevToPoint * xPrevToCur) + (yPrevToPoint * yPrevToCur)) / edgeLength;
 
       // DOCA point
       dcel2d::Point docaPoint(std::get<0>(*firstPointItr) + projection * xPrevToCur / edgeLength,
                               std::get<1>(*firstPointItr) + projection * yPrevToCur / edgeLength,
                               0);
 
       if (projection > edgeLength) docaPoint = *secondPointItr;
       if (projection < 0) docaPoint = *firstPointItr;
 
       double xDocaDist = std::get<0>(vertexPos) - std::get<0>(docaPoint);
       double yDocaDist = std::get<1>(vertexPos) - std::get<1>(docaPoint);
       double docaDist = xDocaDist * xDocaDist + yDocaDist * yDocaDist;
 
       if (docaDist < distToConvexHull) {
         firstHullPointItr = firstPointItr;
         intersection = dcel2d::Coords(std::get<0>(docaPoint), std::get<1>(docaPoint), 0.);
         distToConvexHull = docaDist;
       }
 
       // Check to see if this point is outside the convex hull
       if (isLeft(*firstPointItr, *secondPointItr, vertexPos)) outsideHull = true;
 
       firstPointItr = secondPointItr;
     }
 
     return outsideHull;
   }

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 673 of file Voronoi.cxx.

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

Referenced by makeLeftCircleEvent(), and makeRightCircleEvent().

   {
     // It might be that we don't create a new circle
     IEvent* circle = 0;
 
     // First step is to calculate the center and radius of the circle determined by the three input site events
     dcel2d::Coords center;
     double radius;
     double deltaR;
 
     IEvent* p1 = arc1->getEvent();
     IEvent* p2 = arc2->getEvent();
     IEvent* p3 = arc3->getEvent();
 
     // Compute the center of the circle. Note that this will also automagically check that breakpoints are
     // converging so that this is a circle we want
     if (computeCircleCenter(
           p1->getCoords(), p2->getCoords(), p3->getCoords(), center, radius, deltaR)) {
       double circleBottomX = center[0] - radius;
 
       // Now check if the bottom of this circle lies below the beach line
       if (beachLinePos >= circleBottomX - 10. * deltaR) {
         // Making a circle event!
         dcel2d::Point circleBottom(circleBottomX, center[1], NULL);
 
         fCircleEventList.emplace_back(circleBottom, center);
 
         circle = &fCircleEventList.back();
       }
       else if (circleBottomX - beachLinePos < 1.e-4)
         std::cout << "==> Circle close, beachLine: " << beachLinePos
                   << ", circleBottomX: " << circleBottomX << ", deltaR: " << deltaR
                   << ", d: " << circleBottomX - beachLinePos << std::endl;
     }
 
     return circle;
   }

void voronoi2d::VoronoiDiagram::makeLeftCircleEvent	(	EventQueue &	eventQueue,
		BSTNode *	leaf,
		double	beachLine
	)

private

Definition at line 591 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().

   {
 
     // Check status of triplet of site events to the left of this new leaf
     if (leaf->getPredecessor()) {
       BSTNode* midLeaf = leaf->getPredecessor()->getPredecessor();
 
       if (midLeaf && midLeaf->getPredecessor()) {
         BSTNode* edgeLeaf = midLeaf->getPredecessor()->getPredecessor();
 
         // edge leaves must be different
         if (leaf->getEvent()->getPoint() == edgeLeaf->getEvent()->getPoint()) return;
 
         if (edgeLeaf) {
           IEvent* circleEvent = makeCircleEvent(edgeLeaf, midLeaf, leaf, beachLine);
 
           // Did we succeed in making a circle event?
           if (circleEvent) {
             // Add to the circle node list
             fCircleNodeList.emplace_back(circleEvent);
 
             BSTNode* circleNode = &fCircleNodeList.back();
 
             // If there was an associated circle event to this node, invalidate it
             if (midLeaf->getAssociated()) {
               midLeaf->getAssociated()->setAssociated(NULL);
               midLeaf->getAssociated()->getEvent()->setInvalid();
             }
 
             // Now reset to point at the new one
             midLeaf->setAssociated(circleNode);
             circleNode->setAssociated(midLeaf);
 
             eventQueue.push(circleEvent);
           }
         }
       }
     }
 
     return;
   }

void voronoi2d::VoronoiDiagram::makeRightCircleEvent	(	EventQueue &	eventQueue,
		BSTNode *	leaf,
		double	beachLine
	)

private

Definition at line 633 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().

   {
     // Check status of triplet of site events to the left of this new leaf
     if (leaf->getSuccessor()) {
       BSTNode* midLeaf = leaf->getSuccessor()->getSuccessor();
 
       if (midLeaf && midLeaf->getSuccessor()) {
         BSTNode* edgeLeaf = midLeaf->getSuccessor()->getSuccessor();
 
         if (leaf->getEvent()->getPoint() == edgeLeaf->getEvent()->getPoint()) return;
 
         if (edgeLeaf) {
           IEvent* circleEvent = makeCircleEvent(leaf, midLeaf, edgeLeaf, beachLine);
 
           // Did we succeed in making a circle event?
           if (circleEvent) {
             // Add to the circle node list
             fCircleNodeList.emplace_back(circleEvent);
 
             BSTNode* circleNode = &fCircleNodeList.back();
 
             // If there was an associated circle event to this node, invalidate it
             if (midLeaf->getAssociated()) {
               midLeaf->getAssociated()->setAssociated(NULL);
               midLeaf->getAssociated()->getEvent()->setInvalid();
             }
 
             // Now reset to point at the new one
             midLeaf->setAssociated(circleNode);
             circleNode->setAssociated(midLeaf);
 
             eventQueue.push(circleEvent);
           }
         }
       }
     }
 
     return;
   }

void voronoi2d::VoronoiDiagram::mergeDegenerateVertices ( )

private

merge degenerate vertices (found by zero length edges)

Definition at line 1240 of file Voronoi.cxx.

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

Referenced by terminateInfiniteEdges().

   {
     dcel2d::HalfEdgeList::iterator edgeItr = fHalfEdgeList.begin();
 
     while (edgeItr != fHalfEdgeList.end()) {
       dcel2d::HalfEdge* halfEdge = &(*edgeItr);
       dcel2d::HalfEdge* twinEdge = halfEdge->getTwinHalfEdge();
 
       // Make sure we are not looking at an infinite edge
       if (halfEdge->getTargetVertex() && twinEdge->getTargetVertex()) {
         dcel2d::Coords vtxPosDiff =
           halfEdge->getTargetVertex()->getCoords() - twinEdge->getTargetVertex()->getCoords();
 
         if (vtxPosDiff.norm() < 1.e-3) {
           std::cout << "***>> found a degenerate vertex! "
                     << halfEdge->getTargetVertex()->getCoords()[0] << "/"
                     << halfEdge->getTargetVertex()->getCoords()[1] << ", d: " << vtxPosDiff.norm()
                     << std::endl;
           edgeItr++;
         }
         else
           edgeItr++;
       }
       else
         edgeItr++;
     }
 
     return;
   }

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 854 of file Voronoi.cxx.

Referenced by buildVoronoiDiagram().

   {
     // Need to complete processing of the beachline, the remaning leaves represent the site points with "infinite"
     // edges which we need to terminate at our bounding box.
     // For now let's just do step one which is to find the convex hull
     const BSTNode* node = beachLine.getTopNode();
 
     // Do the convex hull independently but before terminating edges
     getConvexHull(node);
 
     if (node) {
       // Run down the left branches to find the left most leaf
       while (node->getLeftChild())
         node = node->getLeftChild();
 
       // Things to keep track of...
       int nodeCount(0);
       int nBadBreaks(0);
       double lastBreakPoint = std::numeric_limits<double>::lowest();
 
       // Now we are going to traverse across the beach line and process the remaining leaves
       // This will identify the infinite edges and find the convex hull
       while (node) {
         std::cout << "-----------------------------------------------------------------------------"
                      "---------------"
                   << std::endl;
 
         // Do we have a leaf?
         if (!node->getLeftChild() && !node->getRightChild()) {
           std::cout << "node " << nodeCount << " has leaf: " << node
                     << ", face: " << node->getFace()
                     << ", pos: " << node->getEvent()->getCoords()[0] << "/"
                     << node->getEvent()->getCoords()[1] << std::endl;
         }
         // Otherwise it is a break point
         else {
           RootsPair roots;
           double breakPoint = fUtilities.computeBreak(beachLinePos,
                                                       node->getPredecessor()->getEvent(),
                                                       node->getSuccessor()->getEvent(),
                                                       roots);
           double leftArcVal =
             fUtilities.computeArcVal(beachLinePos, breakPoint, node->getPredecessor()->getEvent());
           double rightArcVal =
             fUtilities.computeArcVal(beachLinePos, breakPoint, node->getSuccessor()->getEvent());
 
           dcel2d::HalfEdge* halfEdge = node->getHalfEdge();
           dcel2d::HalfEdge* halfTwin = halfEdge->getTwinHalfEdge();
           dcel2d::Vertex* vertex = halfEdge->getTargetVertex();
 
           std::cout << "node " << nodeCount << " has break: " << breakPoint
                     << ", beachPos: " << beachLinePos << " (last: " << lastBreakPoint
                     << "), leftArcVal: " << leftArcVal << ", rightArcVal: " << rightArcVal
                     << std::endl;
           if (lastBreakPoint > breakPoint)
             std::cout << "     ***** lastBreakPoint larger than current break point *****"
                       << std::endl;
           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;
           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;
 
           const dcel2d::Coords& leftLeafPos = node->getPredecessor()->getEvent()->getCoords();
           const dcel2d::Coords& rightLeafPos = node->getSuccessor()->getEvent()->getCoords();
           Eigen::Vector3f lrPosDiff = rightLeafPos - leftLeafPos;
           Eigen::Vector3f lrPosSum = rightLeafPos + leftLeafPos;
 
           lrPosDiff.normalize();
           lrPosSum *= 0.5;
 
           dcel2d::Coords leafPos = node->getSuccessor()->getEvent()->getCoords();
           dcel2d::Coords leafDir = lrPosDiff;
 
           if (vertex) {
             dcel2d::Coords breakDir(leafDir[1], -leafDir[0], 0.);
             dcel2d::Coords breakPos = lrPosSum;
             dcel2d::Coords vertexPos = vertex->getCoords();
 
             // Get the arclength from the break position to the intersection of the two lines
             dcel2d::Coords breakToLeafPos = leafPos - vertexPos;
             double arcLenToLine = breakToLeafPos.dot(breakDir);
 
             // Now get position
             dcel2d::Coords breakVertexPos = vertexPos + arcLenToLine * breakDir;
 
             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;
 
             // At this point we need to create a vertex and then terminate the half edges here...
             fVertexList.push_back(dcel2d::Vertex(breakVertexPos, halfEdge));
 
             dcel2d::Vertex* breakVertex = &fVertexList.back();
 
             halfTwin->setTargetVertex(breakVertex);
           }
           else {
             std::cout << "****** null vertex!!! Skipping to next node... *********" << std::endl;
           }
 
           if (lastBreakPoint > breakPoint) nBadBreaks++;
 
           lastBreakPoint = breakPoint;
         }
 
         if (node->getAssociated()) {
           BSTNode* associated = node->getAssociated();
           IEvent* iEvent = associated->getEvent();
 
           std::cout << "     -> associated circle: " << iEvent->isCircle()
                     << ", is valid: " << iEvent->isValid() << std::endl;
         }
 
         node = node->getSuccessor();
         nodeCount++;
       }
 
       // Now that we have the convex hull, loop over vertices to see if they are inside or outside of the convex hull
       dcel2d::VertexList::iterator curVertexItr = fVertexList.begin();
 
       size_t nVerticesInitial = fVertexList.size();
 
       // Loop over all vertices to begin with
       while (curVertexItr != fVertexList.end()) {
         // Dereference vertex
         dcel2d::Vertex& vertex = *curVertexItr;
 
         bool outsideHull = !isInsideConvexHull(vertex);
 
         // Do we need to drop this vertex?
         if (outsideHull) { curVertexItr = fVertexList.erase(curVertexItr); }
         else
           curVertexItr++;
       }
 
       mergeDegenerateVertices();
       ComputeFaceArea();
 
       std::cout << "Loop over beachline done, saved " << fConvexHullList.size()
                 << " arcs, encountered " << nBadBreaks << " bad break points" << std::endl;
       std::cout << "-- started with " << nVerticesInitial << " vertices, found "
                 << fVertexList.size() << " inside convex hull" << std::endl;
     }
 
     return;
   }