voronoi2d::VoronoiDiagram Class Reference

VoronoiDiagram class definiton. More...

#include <Voronoi.h>

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

Public Member Functions
	VoronoiDiagram (dcel2d::HalfEdgeList &, dcel2d::VertexList &, dcel2d::FaceList &)
	Constructor. More...
	~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 &, 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 31 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 151 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 149 of file Voronoi.h.

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

private

Definition at line 150 of file Voronoi.h.

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

private

Definition at line 95 of file Voronoi.h.

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

Definition at line 35 of file Voronoi.h.

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

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

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

{
}

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.

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

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

{
    // 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 265 of file Voronoi.cxx.

{
    // 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 690 of file Voronoi.cxx.

{
    // 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 762 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 793 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 1215 of file Voronoi.cxx.

{
    // 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 91 of file Voronoi.cxx.

{
    // 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 1346 of file Voronoi.cxx.

{
    // 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 1424 of file Voronoi.cxx.

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

{
    // 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 76 of file Voronoi.h.

{return fConvexHullList;}

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

private

this function recovers the convex hull

Definition at line 963 of file Voronoi.cxx.

{
    // 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 1040 of file Voronoi.cxx.

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

{return fFaceList;}

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

inline

Recover the list of vertices.

Definition at line 66 of file Voronoi.h.

{return fVertexList;}

double voronoi2d::VoronoiDiagram::getVoronoiDiagramArea ( ) const

inline

recover the area of the convex hull

Definition at line 71 of file Voronoi.h.

{return fVoronoiDiagramArea;}

void voronoi2d::VoronoiDiagram::handleCircleEvents ( BeachLine &  beachLine, EventQueue &  eventQueue, IEvent *  circleEvent  )

private

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

Definition at line 472 of file Voronoi.cxx.

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

{
    // 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 1102 of file Voronoi.cxx.

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

{
    // 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 1131 of file Voronoi.cxx.

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

{
    // 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 557 of file Voronoi.cxx.

{

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

{
    // 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 1187 of file Voronoi.cxx.

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

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

Member Data Documentation

CircleEventList voronoi2d::VoronoiDiagram::fCircleEventList

private

Definition at line 191 of file Voronoi.h.

BSTNodeList voronoi2d::VoronoiDiagram::fCircleNodeList

private

Definition at line 192 of file Voronoi.h.

dcel2d::Coords voronoi2d::VoronoiDiagram::fConvexHullCenter

private

Definition at line 195 of file Voronoi.h.

dcel2d::PointList voronoi2d::VoronoiDiagram::fConvexHullList

private

Definition at line 194 of file Voronoi.h.

dcel2d::FaceList& voronoi2d::VoronoiDiagram::fFaceList

private

Definition at line 187 of file Voronoi.h.

dcel2d::HalfEdgeList& voronoi2d::VoronoiDiagram::fHalfEdgeList

private

Definition at line 185 of file Voronoi.h.

int voronoi2d::VoronoiDiagram::fNumBadCircles

mutableprivate

Definition at line 201 of file Voronoi.h.

dcel2d::PointList voronoi2d::VoronoiDiagram::fPointList

private

Definition at line 189 of file Voronoi.h.

SiteEventList voronoi2d::VoronoiDiagram::fSiteEventList

private

Definition at line 190 of file Voronoi.h.

EventUtilities voronoi2d::VoronoiDiagram::fUtilities

private

Definition at line 204 of file Voronoi.h.

dcel2d::VertexList& voronoi2d::VoronoiDiagram::fVertexList

private

Definition at line 186 of file Voronoi.h.

double voronoi2d::VoronoiDiagram::fVoronoiDiagramArea

private

Definition at line 202 of file Voronoi.h.

double voronoi2d::VoronoiDiagram::fXMax

private

Definition at line 198 of file Voronoi.h.

double voronoi2d::VoronoiDiagram::fXMin

private

Definition at line 197 of file Voronoi.h.

double voronoi2d::VoronoiDiagram::fYMax

private

Definition at line 200 of file Voronoi.h.

double voronoi2d::VoronoiDiagram::fYMin

private

Definition at line 199 of file Voronoi.h.

---

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

• larreco/larreco/RecoAlg/Cluster3DAlgs/Voronoi/Voronoi.h
• larreco/larreco/RecoAlg/Cluster3DAlgs/Voronoi/Voronoi.cxx