Project Scipy Voronoi diagram from 3d to 2d

Question

I am trying to find a way to calculate a 2d Power Diagram in Python. For this I want to make use of the fact that a 2d power diagram can be interpreted as the intersection of a regular 3d voronoi diagram with a plane.

With the SciPy Voronoi module I can calculate a 3d Voronoi diagram - is there a possibility to intersect it with a plane and convert it to a 2d diagram?

Answer 1

There isn't yet a power diagram capability inside of SciPy.

Converting a 3D Voronoi diagram into a 2D power diagram is likely to be difficult and, at least within Python, slow.

To work around this, I've developed a C++ program based around CGAL that can then be wrapped by Python.

With the Python function and C++ code in hand, a power diagram can be generated using:

polys, bbox = GetPowerDiagram(pts)

And the result looks like this if the power diagram is cropped to the bounding box of the point cloud:

and this if it is uncropped (but converted to a finite representation):

The code below worked as of 02019-01-12 using GCC 7.3.0, Python 3.6.7, CGAL 4.11 (1041101000). The code can also be acquired on Github here.

Python code

#!/usr/bin/env python3
#Power Diagramer
#Author: Richard Barnes (rbarnes.org)
#!/usr/bin/env python3
import plumbum
import random
import shapely
import shapely.wkt
from matplotlib import pyplot as plt
from descartes import PolygonPatch

def GetPowerDiagram(points, ray_length=1000, crop=True):
  """Generates a power diagram of a set of points.

  Arguments:

    points - A list of points of the form `[(x,y,weight), (x,y,weight), ...]`

    ray_length - The power diagram contains infinite rays. The direction vector
                 of those rays will be multiplied by `ray_length` and the ends
                 of the rays connected in order to form a finite representation
                 of the polygon

    crop       - If `True`, then the bounded representation above is cropped to
                 the bounding box of the point cloud
  """
  powerd = plumbum.local["./power_diagramer.exe"]

  #Format output for reading by power_diagramer.exe
  points = [map(str,x) for x in points]
  points = [' '.join(x) for x in points]
  points = '\n'.join(points)
  points = '{raylen}\n{crop}\n{points}'.format(
    raylen = ray_length,
    crop   = 'CROP' if crop else 'NOCROP',
    points = points
  )

  #Run the command
  polygons = (powerd["-"] << points)()

  #Get the output of `power_diagramer.exe`. It is in WKT format, one polygon per
  #line.
  polygons = polygons.split("\n")
  polygons = [x.strip() for x in polygons]
  polygons = [x for x in polygons if len(x)>0]
  polygons = [shapely.wkt.loads(x) for x in polygons]

  #Generate bounding box for ease in plotting
  bbox = [x.bounds for x in polygons]
  minx = min([x[0] for x in bbox])
  miny = min([x[1] for x in bbox])
  maxx = max([x[2] for x in bbox])
  maxy = max([x[3] for x in bbox])

  return polygons, (minx,miny,maxx,maxy)

POINT_COUNT = 100
pts         = []
for i in range(POINT_COUNT):
  x      = random.uniform(0,100) 
  y      = random.uniform(0,100)
  weight = random.uniform(0,10)
  pts.append((x,y,weight))

polys, (minx, miny, maxx, maxy) = GetPowerDiagram(pts, ray_length=1000, crop=True)

fig = plt.figure(1, figsize=(5,5), dpi=90)
ax = fig.add_subplot(111)
ax.set_xlim(minx,maxx)
ax.set_ylim(miny,maxy)
for poly in polys:
  ax.add_patch(PolygonPatch(poly))
plt.show()

Makefile

all:
    #-frounding-math is GCC specific, but required for any CGAL code compiled with
    #GCC
    g++ -O3 power_diagram_lib.cpp -o power_diagramer.exe -Wall -lCGAL -lgmp -lgmpxx -lmpfr -frounding-math

And the resulting executable named power_diagramer.exe and placed in the same directory as the Python script.

C++ code

//Finds the cropped Voronoi diagram of a set of points and saves it as WKT
//Compile with: g++ -O3 main.cpp -o power_diagramer.exe -Wall -lCGAL -lgmp
//Author: Richard Barnes (rbarnes.org)
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_filtered_traits_2.h>
#include <CGAL/Regular_triangulation_adaptation_traits_2.h>
#include <CGAL/Regular_triangulation_adaptation_policies_2.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/Voronoi_diagram_2.h>
#include <CGAL/Boolean_set_operations_2.h>
#include <CGAL/bounding_box.h>
#include <CGAL/Polygon_2.h>
#include <iostream>
#include <cstdint>
#include <string>
#include <memory>

typedef CGAL::Exact_predicates_exact_constructions_kernel K;
typedef CGAL::Regular_triangulation_2<K> RT;

typedef CGAL::Regular_triangulation_adaptation_traits_2<RT>         AT;
typedef CGAL::Regular_triangulation_degeneracy_removal_policy_2<RT> DRP;
typedef CGAL::Voronoi_diagram_2<RT, AT, DRP> VD;

int main(int argc, char **argv){
  if(argc!=2){
    std::cerr<<"Synax: "<<argv[0]<<" <FILENAME>"<<std::endl;
    std::cerr<<"<FILENAME> may be a file or '-'. The latter reads from stdin."<<std::endl;
    std::cerr<<"File has the format:"<<std::endl;
    std::cerr<<"<RAY_LENGTH>"        <<std::endl;
    std::cerr<<"<CROP/NOCROP>"       <<std::endl;
    std::cerr<<"<X> <Y> <WEIGHT>"    <<std::endl;
    std::cerr<<"<X> <Y> <WEIGHT>"    <<std::endl;
    std::cerr<<"..."                 <<std::endl;
    std::cerr                        <<std::endl;
    std::cerr<<"'RAY_LENGTH' is a multiplier that extends infinite rays"<<std::endl;
    std::cerr<<"'CROP' will crop the power diagram to the bounding box of the input points"<<std::endl;
    return -1;
  }

  std::string filename = argv[1];

  //Output formatting
  std::cout.precision(4);          //Number of digits of decimal precision
  std::cout.setf(std::ios::fixed); //Don't use scientific notation

  //Used to convert otherwise infinite rays into looooong line segments by
  //multiplying the components of the direction vector of a ray by this value.
  int RAY_LENGTH = 1000;

  //Create a pointer to the correct input stream
  std::istream *fin;
  if(filename=="-")
    fin = &std::cin;
  else
    fin = new std::ifstream(filename);

  std::string crop_string;
  bool do_crop = false;

  (*fin)>>RAY_LENGTH;
  (*fin)>>crop_string;
  if(crop_string=="CROP")
    do_crop = true;
  else if(crop_string=="NOCROP")
    do_crop = false;
  else {
    std::cerr<<"Crop value must be 'CROP' or 'NOCROP'!"<<std::endl;
    return -1;
  }

  //Read in points from the command line
  RT::Weighted_point wp;
  std::vector<RT::Weighted_point> wpoints;
  while((*fin)>>wp)
    wpoints.push_back(wp);

  //Clean up the stream pointer
  if(filename!="-")
    delete fin;

  //Create a Regular Triangulation from the points
  RT rt(wpoints.begin(), wpoints.end());
  rt.is_valid();

  //Wrap the triangulation with a Voronoi diagram adaptor. This is necessary to
  //get the Voronoi faces.
  VD vd(rt);

  //CGAL often returns objects that are either segments or rays. This converts
  //these objects into segments. If the object would have resolved into a ray,
  //that ray is intersected with the bounding box defined above and returned as
  //a segment.
  const auto ConvertToSeg = [&](const CGAL::Object seg_obj, bool outgoing) -> K::Segment_2 {
    //One of these will succeed and one will have a NULL pointer
    const K::Segment_2 *dseg = CGAL::object_cast<K::Segment_2>(&seg_obj);
    const K::Ray_2     *dray = CGAL::object_cast<K::Ray_2>(&seg_obj);
    if (dseg) { //Okay, we have a segment
      return *dseg;
    } else {    //Must be a ray
      const auto &source = dray->source();
      const auto dsx     = source.x();
      const auto dsy     = source.y();
      const auto &dir    = dray->direction();
      const auto tpoint  = K::Point_2(dsx+RAY_LENGTH*dir.dx(),dsy+RAY_LENGTH*dir.dy());
      if(outgoing)
        return K::Segment_2(
          dray->source(),
          tpoint
        );
      else
        return K::Segment_2(
          tpoint,
          dray->source()
        );
    }
  };

  //Loop over the faces of the Voronoi diagram in some arbitrary order
  for(VD::Face_iterator fit = vd.faces_begin(); fit!=vd.faces_end();++fit){
    CGAL::Polygon_2<K> pgon;

    //Edge circulators traverse endlessly around a face. Make a note of the
    //starting point so we know when to quit.
    VD::Face::Ccb_halfedge_circulator ec_start = fit->ccb();

    //Find a bounded edge to start on
    for(;ec_start->is_unbounded();ec_start++){}

    //Current location of the edge circulator
    VD::Face::Ccb_halfedge_circulator ec = ec_start;

    //In WKT format each polygon must begin and end with the same point
    K::Point_2 first_point;

    do {
      //A half edge circulator representing a ray doesn't carry direction
      //information. To get it, we take the dual of the dual of the half-edge.
      //The dual of a half-edge circulator is the edge of a Delaunay triangle.
      //The dual of the edge of Delaunay triangle is either a segment or a ray.
      // const CGAL::Object seg_dual = rt.dual(ec->dual());
      const CGAL::Object seg_dual = vd.dual().dual(ec->dual());

      //Convert the segment/ray into a segment
      const auto this_seg = ConvertToSeg(seg_dual, ec->has_target());

      pgon.push_back(this_seg.source());
      if(ec==ec_start)
        first_point = this_seg.source();

      //If the segment has no target, it's a ray. This means that the next
      //segment will also be a ray. We need to connect those two rays with a
      //segment. The following accomplishes this.
      if(!ec->has_target()){
        const CGAL::Object nseg_dual = vd.dual().dual(ec->next()->dual());
        const auto next_seg = ConvertToSeg(nseg_dual, ec->next()->has_target());
        pgon.push_back(next_seg.target());
      }
    } while ( ++ec != ec_start ); //Loop until we get back to the beginning

    if(do_crop){
      //Find the bounding box of the points. This will be used to crop the Voronoi
      //diagram later.
      const CGAL::Bbox_2 bbox = CGAL::bbox_2(wpoints.begin(), wpoints.end());

      //In order to crop the Voronoi diagram, we need to convert the bounding box
      //into a polygon. You'd think there'd be an easy way to do this. But there
      //isn't (or I haven't found it).
      CGAL::Polygon_2<K> bpoly;
      bpoly.push_back(K::Point_2(bbox.xmin(),bbox.ymin()));
      bpoly.push_back(K::Point_2(bbox.xmax(),bbox.ymin()));
      bpoly.push_back(K::Point_2(bbox.xmax(),bbox.ymax()));
      bpoly.push_back(K::Point_2(bbox.xmin(),bbox.ymax()));

      //Perform the intersection. Since CGAL is very general, it believes the
      //result might be multiple polygons with holes.
      std::list<CGAL::Polygon_with_holes_2<K>> isect;
      CGAL::intersection(pgon, bpoly, std::back_inserter(isect));

      //But we know better. The intersection of a convex polygon and a box is
      //always a single polygon without holes. Let's assert this.
      assert(isect.size()==1);

      //And recover the polygon of interest
      auto &poly_w_holes = isect.front();
      pgon               = poly_w_holes.outer_boundary();
    }

    //Print the polygon as a WKT polygon
    std::cout<<"POLYGON ((";
    for(auto v=pgon.vertices_begin();v!=pgon.vertices_end();v++)
      std::cout<<v->x()<<" "<<v->y()<<", ";
    std::cout<<pgon.vertices_begin()->x()<<" "<<pgon.vertices_begin()->y()<<"))\n";
  }

  return 0;
}