I find that in the boost library, for the voronoi diagram, some edge data is infinite. According to the instruction, it has to be clipped. But I cant find how to do it. Could anybody can give me a example code, please?
Thx
I find that in the boost library, for the voronoi diagram, some edge data is infinite. According to the instruction, it has to be clipped. But I cant find how to do it. Could anybody can give me a example code, please?
Thx
This is a rather stale question, but I found it when trying to solve the exact same problem. It's only fair to share the solution I hit upon so the next poor sap doesn't have to figure it out themselves.
I have no idea if this is a particularly good way to do things, and I suspect it would struggle if you started using curved cells, but it worked OK for my purposes.
The basic problem is that you only have one vertex for the infinite edges, so you have to calculate the direction vector yourself. The direction to use is perpendicular to the vector between the two points separated by the edge.
#include <vector>
#include <boost/polygon/voronoi.hpp>
using boost::polygon::voronoi_builder;
using boost::polygon::voronoi_diagram;
typedef boost::polygon::point_data<int> point;
typedef voronoi_diagram<double>::cell_type cell_type;
typedef voronoi_diagram<double>::edge_type edge_type;
typedef voronoi_diagram<double>::vertex_type vertex_type;
int main(int argc, char *argv[])
{
std::vector<point> points;
// Populate with random points
for (int i = 0; i < 50; i++)
{
points.push_back(point(60 + rand() % 500, 60 + rand() % 500));
}
voronoi_diagram<double> vd;
construct_voronoi(points.begin(), points.end(), &vd);
// vd now contains the voronoi diagram. Let's visualise it
// pseudocode 'draw_line(x1, y1, x2, y2)'
for (voronoi_diagram<double>::const_cell_iterator it = vd.cells().begin();
it != vd.cells().end(); ++it)
{
const cell_type& cell = *it;
const edge_type* edge = cell.incident_edge();
do
{
if (edge->is_primary())
{
// Easy finite edge case
if (edge->is_finite())
{
// Without this check each edge would be drawn twice
// as they are really half-edges
if (edge->cell()->source_index() <
edge->twin()->cell()->source_index())
{
draw_line(edge->vertex0()->x(), edge->vertex0()->y(),
edge->vertex1()->x(), edge->vertex1()->y());
}
}
else
{
// This covers the infinite edge case in question.
const vertex_type* v0 = edge->vertex0();
// Again, only consider half the half-edges, ignore edge->vertex1()
// to avoid overdrawing the lines
if (v0)
{
// Direction of infinite edge if perpendicular to direction
// between the points owning the two half edges.
// Take the rotated right vector and multiply by a large
// enough number to reach your bounding box
point p1 = points[edge->cell()->source_index()];
point p2 = points[edge->twin()->cell()->source_index()];
int end_x = (p1.y() - p2.y()) * 640;
int end_y = (p1.x() - p2.x()) * -640;
draw_line(v0->x(), v0->y(),
end_x, end_y);
}
}
}
edge = edge->next();
} while (edge != cell.incident_edge());
}
}
I found this code segment here: http://www.boost.org/doc/libs/1_55_0/libs/polygon/example/voronoi_visualizer.cpp
void clip_infinite_edge(
const edge_type& edge, std::vector<point_type>* clipped_edge) {
const cell_type& cell1 = *edge.cell();
const cell_type& cell2 = *edge.twin()->cell();
point_type origin, direction;
// Infinite edges could not be created by two segment sites.
if (cell1.contains_point() && cell2.contains_point()) {
point_type p1 = retrieve_point(cell1);
point_type p2 = retrieve_point(cell2);
origin.x((p1.x() + p2.x()) * 0.5);
origin.y((p1.y() + p2.y()) * 0.5);
direction.x(p1.y() - p2.y());
direction.y(p2.x() - p1.x());
} else {
origin = cell1.contains_segment() ?
retrieve_point(cell2) :
retrieve_point(cell1);
segment_type segment = cell1.contains_segment() ?
retrieve_segment(cell1) :
retrieve_segment(cell2);
coordinate_type dx = high(segment).x() - low(segment).x();
coordinate_type dy = high(segment).y() - low(segment).y();
if ((low(segment) == origin) ^ cell1.contains_point()) {
direction.x(dy);
direction.y(-dx);
} else {
direction.x(-dy);
direction.y(dx);
}
}
coordinate_type side = xh(brect_) - xl(brect_);
coordinate_type koef =
side / (std::max)(fabs(direction.x()), fabs(direction.y()));
if (edge.vertex0() == NULL) {
clipped_edge->push_back(point_type(
origin.x() - direction.x() * koef,
origin.y() - direction.y() * koef));
} else {
clipped_edge->push_back(
point_type(edge.vertex0()->x(), edge.vertex0()->y()));
}
if (edge.vertex1() == NULL) {
clipped_edge->push_back(point_type(
origin.x() + direction.x() * koef,
origin.y() + direction.y() * koef));
} else {
clipped_edge->push_back(
point_type(edge.vertex1()->x(), edge.vertex1()->y()));
}
}
It might be missing some class variables or methods, but the logic is what is important here.
I wrote a program to clip the infinite edges in Voronoi cells so the edges are confined inside a rectangular bounding box. The code is in https://github.com/datajaguar/jaguardb/blob/master/src/JagCGAL.cc
getIntersectionPointWithBox() method basically checks each side ( total 4 sides) of the bounding box and gets the intersection point if there is one. fillInCorners() method checks missing points around the bounding box and add the corner coordinates to the edges. To make a polygon for each cell, it adds the first point at the end to make the polygon closed.