如何使用 CGAL 获取多边形的中轴？
How do you get the medial axis of a multipolygon using CGAL?
我想使用 CGAL 提取多边形的中轴。
查看 CGAL 文档我只看到了制作 straight skeleton 的函数。
如何使用 CGAL 获取中轴？
概览
如果我们有一个多边形（具有直边的平面域，可能是凸边），那么形成多边形中轴的边是多边形 Voronoi diagram.
的子集
那么，我们的策略是找到多边形的 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边都是中轴的一部分。事实证明，那些部分由 凹面 顶点定义的边（又名“reflex vertices") of the polygon are not part of the medial axis. The concave vertices of the polygon are those that "point inwards". Computationally, we can find them by using cross-products. (Incidentally, the area of a polygon can be calculated by summing the cross-products of each of its vertices。）
移除上面讨论的边缘给我们最终的中轴：
代码
以下 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
// 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()
我想使用 CGAL 提取多边形的中轴。
查看 CGAL 文档我只看到了制作 straight skeleton 的函数。
如何使用 CGAL 获取中轴？
概览
如果我们有一个多边形（具有直边的平面域，可能是凸边），那么形成多边形中轴的边是多边形 Voronoi diagram.
的子集那么，我们的策略是找到多边形的 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边都是中轴的一部分。事实证明，那些部分由 凹面 顶点定义的边（又名“reflex vertices") of the polygon are not part of the medial axis. The concave vertices of the polygon are those that "point inwards". Computationally, we can find them by using cross-products. (Incidentally, the area of a polygon can be calculated by summing the cross-products of each of its vertices。）
移除上面讨论的边缘给我们最终的中轴：
代码
以下 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
// 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()