c++ - 如何使用 CGAL 获得多面体的中轴？

Question

我想使用 CGAL 提取 MultiPolygon 的中轴。

浏览 CGAL 文档，我只看到制作直骨架的函数。

如何使用 CGAL 获得中轴？

标签： c++computational-geometrycgal

Answer 1

概述

如果我们有一个多边形（一个带有直边的平面域，可能是凸的），那么形成多边形中轴的边是多多边形Voronoi 图的子集。

因此，我们的策略是找到多边形的 Voronoi 图，然后从图中删除边，直到剩下中轴。如果我们有一个 MultiPolygon，我们只需对其每个组成多边形重复该过程。

我将在下面直观地演示该过程，然后提供代码来完成它。我将使用的测试多边形的数据是这个WKT：

MULTIPOLYGON(((223.83073496659313 459.9703924976702,256.1247216035629 304.08821449578033,1.1135857461033538 187.63424160514955,468.8195991091309 189.21374898389445,310.69042316258424 318.7630937196408,432.07126948775 451.93407043677666,453.2293986636971 612.1086753947116,32.29398663697134 616.325100949687,223.83073496659313 459.9703924976702)))

多边形看起来像这样：

多边形的完整 Voronoi 图如下所示：

在这一点上，我们可以考虑 Voronoi 图中的三种类型的顶点：正确位于多边形内的顶点、位于其边界上的顶点以及正确位于多边形外的顶点。很明显，中轴不包含 Voronoi 图的任何边缘，该边缘部分由位于多边形外部的点定义。因此，我们删除了这些外部边缘。留给我们这个：

并非上图中的所有 Voronoi 边都是中轴的一部分。事实证明，部分由多边形的凹顶点（又称“反射顶点”）定义的那些边不是中轴的一部分。多边形的凹顶点是“指向内”的顶点。在计算上，我们可以通过使用叉积来找到它们。（顺便说一下，多边形的面积可以通过对每个顶点的叉积求和来计算。）

去除上面讨论的边缘给我们最终的中轴：

编码

以下 C++ 代码计算多面体的中轴。

// Compile with: clang++ -DBOOST_ALL_NO_LIB -DCGAL_USE_GMPXX=1 -O3 -g -Wall -Wextra -pedantic -march=native -frounding-math main.cpp -lgmpxx -lmpfr -lgmp
#include <CGAL/Boolean_set_operations_2.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/IO/WKT.h>
#include <CGAL/Polygon_with_holes_2.h>
#include <CGAL/Segment_Delaunay_graph_2.h>
#include <CGAL/Segment_Delaunay_graph_adaptation_policies_2.h>
#include <CGAL/Segment_Delaunay_graph_adaptation_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_traits_2.h>
#include <CGAL/squared_distance_2.h>
#include <CGAL/Voronoi_diagram_2.h>

#include <algorithm>
#include <deque>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <set>
#include <stdexcept>
#include <unordered_set>

typedef CGAL::Exact_predicates_exact_constructions_kernel              K;
typedef CGAL::Segment_Delaunay_graph_traits_2<K>                       Gt;
typedef CGAL::Segment_Delaunay_graph_2<Gt>                             SDG2;
typedef CGAL::Segment_Delaunay_graph_adaptation_traits_2<SDG2>         AT;
typedef CGAL::Segment_Delaunay_graph_degeneracy_removal_policy_2<SDG2> AP;
typedef CGAL::Voronoi_diagram_2<SDG2, AT, AP>      VoronoiDiagram;
typedef AT::Site_2                                 Site_2;
typedef AT::Point_2                                Point_2;
typedef VoronoiDiagram::Locate_result              Locate_result;
typedef VoronoiDiagram::Vertex_handle              Vertex_handle;
typedef VoronoiDiagram::Face_handle                Face_handle;
typedef VoronoiDiagram::Halfedge_handle            Halfedge_handle;
typedef VoronoiDiagram::Ccb_halfedge_circulator    Ccb_halfedge_circulator;
typedef VoronoiDiagram::Bounded_halfedges_iterator BHE_Iter;
typedef VoronoiDiagram::Halfedge                   Halfedge;
typedef VoronoiDiagram::Vertex                     Vertex;
typedef CGAL::Polygon_with_holes_2<K> Polygon;
typedef std::deque<Polygon> MultiPolygon;

/// Creates a hash of a Point_2, used for making O(1) point lookups
// struct Point2Hash {
//   size_t operator()(const Point_2 &pt) const {
//     std::hash<double> hasher;
//     auto seed = hasher(pt.x());
//     // boost::hash_combine from https://stackoverflow.com/q/35985960/752843
//     seed ^= hasher(pt.y()) + 0x9e3779b9 + (seed<<6) + (seed>>2);
//     return seed;
//   }
// };

typedef std::set<Point_2> Point2_Set;
typedef std::map<Vertex_handle, int> VH_Int_Map;


/// Holds a more accessible description of the Voronoi diagram
struct MedialData {
  /// Map of vertices comprising the Voronoi diagram
  VH_Int_Map vertex_handles;
  /// List of edges in the diagram (pairs of the vertices above)
  std::vector<std::pair<int, int>> edges;
  /// Medial axis up governor. 1:1 correspondance with edges above.
  std::vector<VoronoiDiagram::Delaunay_graph::Vertex_handle> ups;
  /// Medial axis down governor. 1:1 correspondance with edges above.
  std::vector<VoronoiDiagram::Delaunay_graph::Vertex_handle> downs;
};


/// Read well-known text from @p filename to obtain shape boundary
MultiPolygon get_wkt_from_file(const std::string& filename){
  std::ifstream fin(filename);
  MultiPolygon mp;
  CGAL::read_multi_polygon_WKT(fin, mp);

  if(mp.empty()){
    throw std::runtime_error("WKT file '" + filename + "' was empty!");
  }
  for(const auto &poly: mp){
    if(poly.outer_boundary().size()==0){
      throw std::runtime_error("WKT file '" + filename + "' contained a polygon without an outer boundary!");
    }
  }

  return mp;
}


/// Converts a MultiPolygon into its corresponding Voronoi diagram
VoronoiDiagram convert_mp_to_voronoi_diagram(const MultiPolygon &mp){
  VoronoiDiagram vd;

  const auto add_segments_to_vd = [&](const auto &poly){
    for(std::size_t i=0;i<poly.size();i++){
      std::cerr<<i<<" "<<std::fixed<<std::setprecision(10)<<poly[i]<<std::endl;
      // Modulus to close the loop
      vd.insert(
        Site_2::construct_site_2(poly[i], poly[(i+1)%poly.size()])
      );
    }
  };

  for(const auto &poly: mp){                    // For each polygon in MultiPolygon
    std::cout<<poly<<std::endl;                 // Print polygon to screen for debugging
    add_segments_to_vd(poly.outer_boundary());  // Add the outer boundary
    for(const auto &hole : poly.holes()){       // And any holes
      add_segments_to_vd(hole);
    }
  }

  if(!vd.is_valid()){
    throw std::runtime_error("Voronoi Diagram was not valid!");
  }

  return vd;
}


/// Find @p item in collection @p c or add it if not present.
/// Returns the index of `item`'s location
int find_or_add(VH_Int_Map &c, const Vertex_handle &item){
  // Map means we can do this in log(N) time
  if(c.count(item) == 0){
    c.emplace(item, c.size());
    return c.size() - 1;
  }

  return c.at(item);
}


/// Convert a map of <T, int> pairs to a vector of `T` ordered by increasing int
std::vector<Vertex_handle> map_to_ordered_vector(const VH_Int_Map &m){
  std::vector<std::pair<Vertex_handle, int>> to_sort(m.begin(), m.end());
  to_sort.reserve(m.size());
  std::sort(to_sort.begin(), to_sort.end(), [](const auto &a, const auto &b){
    return a.second < b.second;
  });

  std::vector<Vertex_handle> ret;
  ret.reserve(to_sort.size());
  std::transform(begin(to_sort), end(to_sort), std::back_inserter(ret),
    [](auto const& pair){ return pair.first; }
  );

  return ret;
}


/// Find vertex handles which are in the interior of the MultiPolygon
std::set<Vertex_handle> identify_vertex_handles_inside_mp(
  const VoronoiDiagram &vd,
  const MultiPolygon &mp
){
  // Used to accelerate interior lookups by avoiding Point-in-Polygon checks for
  // vertices we've already considered
  std::set<Vertex_handle> considered;
  // The set of interior vertices we are building
  std::set<Vertex_handle> interior;

  for (
      auto edge_iter = vd.bounded_halfedges_begin();
      edge_iter != vd.bounded_halfedges_end();
      edge_iter++
  ) {
    // Determine if an orientation implies an interior vertex
    const auto inside = [](const auto &orientation){
      return orientation == CGAL::ON_ORIENTED_BOUNDARY || orientation == CGAL::POSITIVE;
    };

    // Determine if a vertex is in the interior of the multipolygon and, if so,
    // add it to `interior`
    const auto vertex_in_mp_interior = [&](const Vertex_handle& vh){
      // Skip vertices which have already been considered, since a vertex may
      // be connected to multiple halfedges
      if(considered.count(vh)!=0){
        return;
      }
      // Ensure we don't look at a vertex twice
      considered.insert(vh);
      // Determine if the vertex is inside of any polygon of the MultiPolygon
      const auto inside_of_a_poly = std::any_of(
        mp.begin(), mp.end(), [&](const auto &poly) {
          return inside(CGAL::oriented_side(vh->point(), poly));
        }
      );
      // If the vertex was inside the MultiPolygon make a note of it
      if(inside_of_a_poly){
        interior.insert(vh);
      }
    };

    // Check both vertices of the current halfedge of the Voronoi diagram
    vertex_in_mp_interior(edge_iter->source());
    vertex_in_mp_interior(edge_iter->target());
  }

  return interior;
}


/// The medial axis is formed by building a Voronoi diagram and then removing
/// the edges of the diagram which connect to the concave points of the
/// MultiPolygon. Here, we identify those concave points
Point2_Set identify_concave_points_of_mp(const MultiPolygon &mp){
  Point2_Set concave_points;

  // Determine cross-product, given three points. The sign of the cross-product
  // determines whether the point is concave or convex.
  const auto z_cross_product = [](const Point_2 &pt1, const Point_2 &pt2, const Point_2 &pt3){
    const auto dx1 = pt2.x() - pt1.x();
    const auto dy1 = pt2.y() - pt1.y();
    const auto dx2 = pt3.x() - pt2.x();
    const auto dy2 = pt3.y() - pt2.y();
    return dx1 * dy2 - dy1 * dx2;
  };

  // Loop through all the points in a polygon, get their cross products, and
  // add any concave points to the set we're building.
  // `sense` should be `1` for outer boundaries and `-1` for holes, since holes
  // will have points facing outward.
  const auto consider_polygon = [&](const auto &poly, const double sense){
    for(size_t i=0;i<poly.size()+1;i++){
      const auto zcp = z_cross_product(
        poly[(i + 0) % poly.size()],
        poly[(i + 1) % poly.size()],
        poly[(i + 2) % poly.size()]
      );
      if(sense*zcp < 0){
        concave_points.insert(poly[(i + 1) % poly.size()]);
      }
    }
  };

  // Loop over the polygons of the MultiPolygon, as well as their holes
  for(const auto &poly: mp){
    // Outer boundary has positive sense
    consider_polygon(poly.outer_boundary(), 1);
    for(const auto &hole: poly.holes()){
      // Inner boundaries (holes) have negative (opposite) sense
      consider_polygon(hole, -1);
    }
  }

  return concave_points;
}


/// Print the points which collectively comprise the medial axis
void print_medial_axis_points(const MedialData &md, const std::string &filename){
  std::ofstream fout(filename);
  fout<<"x,y"<<std::endl;
  for (const auto &vh : map_to_ordered_vector(md.vertex_handles)) {
    fout << vh->point().x() << "," << vh->point().y() << std::endl;
  }
}


/// Prints the edges of the medial diagram
void print_medial_axis_edges(const MedialData &md, const std::string &filename){
  std::ofstream fout(filename);
  fout<<"SourceIdx,TargetIdx,UpGovernorIsPoint,DownGovernorIsPoint"<<std::endl;
  for (std::size_t i = 0; i < md.edges.size(); i++) {
    fout << md.edges[i].first        << ","
          << md.edges[i].second      << ","
          << md.ups[i]->is_point()   << "," // Is up-governor a point?
          << md.downs[i]->is_point()        // Is down-governor a point?
          << std::endl;
  }
}


MedialData filter_voronoi_diagram_to_medial_axis(
  const VoronoiDiagram &vd,
  const MultiPolygon &mp
){
  MedialData ret;

  const auto interior = identify_vertex_handles_inside_mp(vd, mp);
  const auto concave_points = identify_concave_points_of_mp(mp);

  // Returns true if a point is a concave point of the MultiPolygon
  const auto pconcave = [&](const Point_2 &pt){
    return concave_points.count(pt) != 0;
  };

  // The Voronoi diagram is comprised of a number of vertices connected by lines
  // Here, we go through each edge of the Voronoi diagram and determine which
  // vertices it's incident on. We add these vertices to `ret.vertex_handles`
  // so that they will have unique ids.

  // The `up` and `down` refer to the medial axis governors - that which
  // constrains each edge of the Voronoi diagram
  for (
      auto edge_iter = vd.bounded_halfedges_begin();
      edge_iter != vd.bounded_halfedges_end();
      edge_iter++
  ) {
    const Halfedge& halfedge = *edge_iter;
    const Vertex_handle& v1p = halfedge.source();
    const Vertex_handle& v2p = halfedge.target();

    // Filter Voronoi diagram to only the part in the interior of the
    // MultiPolygon
    if(interior.count(v1p)==0 || interior.count(v2p)==0){
      continue;
    }

    // Drop those edges of the diagram which are not part of the medial axis
    if(pconcave(v1p->point()) || pconcave(v2p->point())){
      continue;
    }

    // Get unique ids for edge vertex handle that's part of the medial axis
    const auto id1 = find_or_add(ret.vertex_handles, v1p);
    const auto id2 = find_or_add(ret.vertex_handles, v2p);
    ret.edges.emplace_back(id1, id2);

    // Keep track of the medial axis governors
    ret.ups.push_back(halfedge.up());
    ret.downs.push_back(halfedge.down());
  }

  return ret;
}


int main(int argc, char** argv) {
  if(argc!=2){
      std::cerr<<"Syntax: "<<argv[0]<<" <Shape Boundary WKT>"<<std::endl;
      return -1;
  }

  CGAL::set_pretty_mode(std::cout);

  const auto mp = get_wkt_from_file(argv[1]);
  const auto voronoi = convert_mp_to_voronoi_diagram(mp);
  const auto ma_data = filter_voronoi_diagram_to_medial_axis(voronoi, mp);

  print_medial_axis_points(ma_data, "voronoi_points.csv");
  print_medial_axis_edges(ma_data, "voronoi_edges.csv");

  return 0;
}

您可以使用此 Python 脚本绘制生成的中轴：

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import numpy as np
from shapely import wkt

fig, ax = plt.subplots()

# Output file from C++ medial axis code
pts = np.loadtxt("build/voronoi_points.csv", skiprows=1, delimiter=',')

fig4 = wkt.loads(open("fig4_data.wkt").read())
for geom in fig4.geoms:
  xs, ys = geom.exterior.xy
  ax.plot(xs, ys, '-ok', lw=4)

# Output file from C++ medial axis code
ma = np.loadtxt("build/voronoi_edges.csv", dtype=int, skiprows=1, delimiter=',')
for mal in ma:
  print(mal)
  ax.plot((pts[mal[0]][0], pts[mal[1]][0]), (pts[mal[0]][1], pts[mal[1]][1]), '-o')

plt.show()