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::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	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 &, 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::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 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 148 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 146 of file Voronoi.h.

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

private

Definition at line 147 of file Voronoi.h.

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

private

Definition at line 92 of file Voronoi.h.

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

Definition at line 37 of file Voronoi.h.

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

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.

                                                                                                                        :
     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 ( )

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

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

 {
     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 267 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 692 of file Voronoi.cxx.

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

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,
		double &	delta
	)		const

private

Definition at line 764 of file Voronoi.cxx.

 {
     // 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,
		double &	delta
	)		const

private

Definition at line 795 of file Voronoi.cxx.

 {
     // 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 1155 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().

 {
     // 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)
         {
             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();
             
             // 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();
                 doNext   = false;
             }
             
             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 93 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 1279 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 1357 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 1296 of file Voronoi.cxx.

References fPointList, isLeft(), and max.

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

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

Referenced by terminateInfiniteEdges().

78 {return fConvexHullList;}

voronoi2d::VoronoiDiagram::fConvexHullList

dcel2d::PointList fConvexHullList

Definition: Voronoi.h:191

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

 {
     // 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;
 }

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;}

voronoi2d::VoronoiDiagram::fFaceList

dcel2d::FaceList & fFaceList

Definition: Voronoi.h:184

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;}

voronoi2d::VoronoiDiagram::fVertexList

dcel2d::VertexList & fVertexList

Definition: Voronoi.h:183

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;}

voronoi2d::VoronoiDiagram::fVoronoiDiagramArea

double fVoronoiDiagramArea

Definition: Voronoi.h:199

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

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

 {
     // 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 1042 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 84 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 1071 of file Voronoi.cxx.

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

 {
     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 655 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 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().

 {
     
     // 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 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().

 {
     // 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 1127 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 831 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;
 }