从3d到2d的Project Scipy Voronoi图 [英] Project Scipy Voronoi diagram from 3d to 2d

问题描述

我正在尝试找到一种在Python中计算2d Power Diagram的方法。为此，我想利用以下事实：

，如果未裁剪（但转换为有限表示），则为：

以下代码的工作原理是使用GCC 7.3.0，Python 3.6.7，CGAL 4.11（1041101000）的02019-01-12版本。也可以在Github的此处上获取代码。

Python代码

 ＃！ / usr / bin / env python3 
 #Power Diagramer 
＃作者：Richard Barnes（rbarnes.org）
＃！/ usr / bin / env python3 
 import plumum 
随机进口
匀称进口
匀称进口.wkt 
从matplotlib进口pyplot as plt 
从笛卡尔进口进口PolygonPatch 
 
 def GetPowerDiagram（points，ray_length = 1000，crop = True）：
生成一组点的功率图。
 
参数：
 
 points-点的列表形式`[[（x，y，weight），（x，y，weight），...]`
 
 ray_length-功率图包含无限射线，其方向矢量
射线将被乘以ray_length和射线的末端
为了形成多边形
 
的有限表示形式
裁剪-如果为 True，则将上面的有界表示形式裁剪为
点云$ b $的边界框b 
供电= plumbum.local [ ./ power_diagramer.exe] 
 
＃通过power_diagramer.exe读取的格式化输出
点= [map（str， x）表示x的积分] 
积分= [''.join（x）表示x的积分] 
积分='\n'.join（points）
积分='{ raylen} \n {crop} \n {points}'。format（
 raylen = ray_length，
作物='CROP'如果作物为'NOCROP'，
积分=点
）
 
＃运行命令
polygons =（powerd [-]<< points）（）
 
＃获取`power_diagramer.exe`的输出。它为WKT格式，每
 #line个多边形。 
的多边形= polygons.split（ \n）
的多边形= [x.strip（）表示x的多边形] 
多边形= [x表示x的多边形，如果len（x） > 0] 
多边形= [多边形中x的shape.wkt.loads（x）] 
 
＃生成边界框以便于绘制
 bbox = [x.bounds对于多边形中的x] 
 minx = min（[对于bbox中x的x [0 [0]]）
 miny = min（[对于bbox中x的x [1]]]）
 maxx = max （[xbbox中x的x [2]]）
 maxy = max（[bbox中x中的x的[x [3]]）
 
返回多边形，（minx，miny，maxx，maxy ）
 
 POINT_COUNT = 100 
 pts = [] 
对于范围内的i（POINT_COUNT）：
x = random.uniform（0,100）
y = random。均匀（0,100）
重量= random.uniform（0,10）
 pts.append（（x，y，weight））
 
多边形，（minx，miny，maxx ，maxy）= GetPowerDiagram（pts，ray_length = 1000，crop = True）
 
图= plt.figure（1，图大小=（5,5），dpi = 90）
轴= fig.add_subplot（111）
 ax.set_xlim（minx，maxx）
 ax.set_ylim（miny，maxy）
 for pol中的poly ys：
 ax.add_patch（PolygonPatch（poly））
 plt.show（）
  

Makefile

  all：
＃-rounding-math特定于GCC，但对于任何编译的CGAL代码都是必需的与
 #GCC 
 g ++ -O3 power_diagram_lib.cpp -o power_diagramer.exe -Wall -lCGAL -lgmp -lgmpxx -lmpfr -frounding-math 
  

生成的可执行文件名为 power_diagramer.exe ，并与Python脚本放在同一目录中。

C ++代码

  //点的Voronoi图并将其保存为WKT 
 //编译为：g ++ -O3 main.cpp -o power_diagramer.exe -Wall -lCGAL -lgmp 
 //作者：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>在; 
 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>可以是文件或'-'。后者从stdin中读取。<< std :: endl; 
 std :: cerr<<文件具有以下格式：<< 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'是扩展无限射线的倍增器。 
 std :: cerr<<" CROP会将功率图裁剪到输入点的边界框中。 
返回-1; 
} 
 
 std :: string filename = argv [1]; 
 
 //输出格式
 std :: cout.precision（4）; //十进制精度的位数
 std :: cout.setf（std :: ios :: fixed）; //请勿使用科学记数法
 
 //用于通过
 //将光线的方向向量的分量乘以该值，将否则的无限光线转换为线段。 
 int RAY_LENGTH = 1000; 
 
 //创建一个指向正确输入流的指针
 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'或'NOCROP'！ 
返回-1; 
} 
 
 //从命令行中读取点
 RT :: Weighted_point wp; 
 std :: vector< RT :: Weighted_point>积分
 while（（* fin）>> wp）
 wpoints.push_back（wp）; 
 
 //清理流指针
 if（filename！=-）
 delete fin; 
 
 //从点
中创建规则三角剖分RT rt（wpoints.begin（），wpoints.end（））; 
 rt.is_valid（）; 
 
 //用Voronoi图适配器包装三角剖分。这对于
 //获得Voronoi的脸是必需的。 
 VD vd（rt）; 
 
 // CGAL通常返回段或射线的对象。这会将
 //这些对象转换成段。如果对象将分解为射线，则
 //将该射线与上面定义的边界框相交，并以
 //一个段返回。 
 const auto ConvertToSeg = [&]（const CGAL :: Object seg_obj，bool支出）-> K :: Segment_2 {
 //其中之一将成功，并且其中一个将具有NULL指针
 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）{//好，我们有一个段
 return * dseg; 
} else {//必须是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 
返回K :: Segment_2（
 tpoint，
 dray-> source（）
）; 
} 
}; 
 
 //以任意顺序环绕Voronoi图的表面
 for（VD :: Face_iterator fit = vd.faces_begin（）; fit！= vd.faces_end（）; + + fit）{
 CGAL :: Polygon_2< K> gon 
 
 //边缘循环器无休止地绕着脸移动。记下
 //起点，以便我们知道何时退出。 
 VD :: Face :: Ccb_halfedge_circulator ec_start = fit-> ccb（）; 
 
 //找到有界边以
开始，用于（; ec_start-> is_unbounded（）; ec_start ++）{} 
 
 //边缘循环器
 VD :: Face :: Ccb_halfedge_circulator ec = ec_start; 
 
 //以WKT格式，每个多边形必须以相同的点开始和结束
 K :: Point_2 first_point; 
 
 do {
 //代表光线的半边环行器不会携带方向
 //信息。为了得到它，我们采用半边对边的对边。 
 //半边循环器的对偶是Delaunay三角形的边。 
 // Delaunay三角形的边缘对角线是线段或射线。 
 // const CGAL :: Object seg_dual = rt.dual（ec-> dual（））; 
 const CGAL :: Object seg_dual = vd.dual（）。dual（ec-> dual（））; 
 
 //将线段/射线转换为线段
 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（！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）; //循环直到我们回到开头
 
 if（do_crop）{
 //找到点的边界框。这将在以后裁剪Voronoi 
 //图表时使用。 
 const CGAL :: Bbox_2 bbox = CGAL :: bbox_2（wpoints.begin（），wpoints.end（））; 
 
 //为了裁剪Voronoi图，我们需要将边界框
 //转换成多边形。您会认为有一种简单的方法可以做到这一点。但是没有
 //（或我没有找到它）。 
 CGAL :: Polygon_2< K>保利
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（）））; 
 
 //执行交点。由于CGAL非常笼统，因此它认为
 //结果可能是带有孔的多个多边形。 
 std :: list< CGAL :: Polygon_with_holes_2< K>岛
 CGAL :: intersection（pgon，bpoly，std :: back_inserter（isect））; 
 
 //但是我们知道的更好。凸多边形和盒子的交集是
 //始终是没有孔的单个多边形。断言一下。 
 assert（isect.size（）== 1）; 
 
 //然后恢复感兴趣的多边形
 auto& poly_w_holes = isect.front（）; 
 pgon = poly_w_holes.outer_boundary（）; 
} 
 
 //将多边形打印为WKT多边形
 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; 
} 
 
返回0； 
} 
  

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?

解决方案

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

这篇关于从3d到2d的Project Scipy Voronoi图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！