当前位置：首页 > news >正文

Mesh地形曲面提取等高线

news 2025/6/17 20:14:42

Mesh地形曲面提取等高线

本文演示如何提取等值高度来构建地形图。

示例

第一步是从三角剖分的所有面中提取经过该面的等值面，以线段的形式。下面的函数可以测试一个isovalue是否穿过一个三角面，并提取它:

利用这些函数，我们可以创建一个线段图，稍后再处理成一组折线。为此，我们使用boost::adjacency_list结构，并记录从端点位置到图顶点的映射关系。

下面的代码计算了50个均匀分布在点云最小和最大高度之间的等值线，并创建了一个包含所有等值线的图

#include<iostream>
#include<cmath>#include<CGAL/Surface_mesh.h>
#include<CGAL/Surface_mesh/IO/PLY.h>#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Projection_traits_xy_3.h>
#include <CGAL/Delaunay_triangulation_2.h>#include <boost/graph/adjacency_list.hpp>
#include <CGAL/boost/graph/split_graph_into_polylines.h>
#include <CGAL/IO/WKT.h>#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
#include <CGAL/Polyline_simplification_2/simplify.h>
#include <CGAL/Polyline_simplification_2/Squared_distance_cost.h>#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/compute_average_spacing.h>using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;
using Point_2 = Kernel::Point_2;
using Point_3 = Kernel::Point_3;
using Segment_3 = Kernel::Segment_3;
// Triangulated Irregular Network
// 不规则三角网
using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;
using Mesh = CGAL::Surface_mesh<Point_3>;//
namespace PS = CGAL::Polyline_simplification_2;
using CDT_vertex_base = PS::Vertex_base_2<Projection_traits>;
using CDT_face_base = CGAL::Constrained_triangulation_face_base_2<Projection_traits>;
using CDT_TDS = CGAL::Triangulation_data_structure_2<CDT_vertex_base, CDT_face_base>;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<Projection_traits, CDT_TDS>;
using CTP = CGAL::Constrained_triangulation_plus_2<CDT>;#ifdef CGAL_LINKED_WITH_TBB
using Concurrency_tag = CGAL::Parallel_tag;
#else
using Concurrency_tag = CGAL::Sequential_tag;
#endif//判断三角面是否与某个等值交叉
bool face_has_isovalue(TIN::Face_handle fh, double isovalue)
{bool above = false, below = false;for (int i = 0; i < 3; ++i){// Face has isovalue if one of its vertices is above and another one below// 如果面的其中一个顶点位于上方，另一个顶点位于下方，则面具有等值if (fh->vertex(i)->point().z() > isovalue)above = true;if (fh->vertex(i)->point().z() < isovalue)below = true;}return (above && below);
}
//在三角面上抽取某个等值线段
Segment_3 isocontour_in_face(TIN::Face_handle fh, double isovalue)
{Point_3 source;Point_3 target;bool source_found = false;for (int i = 0; i < 3; ++i){Point_3 p0 = fh->vertex((i + 1) % 3)->point();Point_3 p1 = fh->vertex((i + 2) % 3)->point();// Check if the isovalue crosses segment (p0,p1)//检查等值是否跨越线段 （p0，p1）即，三角形的某一条边if ((p0.z() - isovalue) * (p1.z() - isovalue) > 0)continue;double zbottom = p0.z();double ztop = p1.z();if (zbottom > ztop){std::swap(zbottom, ztop);std::swap(p0, p1);}// Compute position of segment vertex// 计算线段顶点的位置double ratio = (isovalue - zbottom) / (ztop - zbottom);Point_3 p = CGAL::barycenter(p0, (1 - ratio), p1, ratio);if (source_found)target = p;else{source = p;source_found = true;}}return Segment_3(source, target);
}template <typename Graph>
class Polylines_visitor
{
private:std::vector<std::vector<Point_3> >& polylines;Graph& graph;
public:Polylines_visitor(Graph& graph, std::vector<std::vector<Point_3> >& polylines): polylines(polylines), graph(graph) { }void start_new_polyline(){polylines.push_back(std::vector<Point_3>());}void add_node(typename Graph::vertex_descriptor vd){polylines.back().push_back(graph[vd]);}void end_polyline(){// filter small polylines// 过滤小的折线if (polylines.back().size() < 50)polylines.pop_back();}
};
int main() {Mesh mesh;CGAL::IO::read_PLY("dtm.ply", mesh);CGAL::Bbox_3 bbox = CGAL::bbox_3(mesh.points().begin(), mesh.points().end());std::cout << "with:" << bbox.x_span() << "\nheight:" << bbox.y_span() << std::endl;TIN tin(mesh.points().begin(), mesh.points().end());std::array<double, 50> isovalues; // Contour 50 isovaluesfor (std::size_t i = 0; i < isovalues.size(); ++i)isovalues[i] = bbox.zmin() + ((i + 1) * (bbox.zmax() - bbox.zmin()) / (isovalues.size() - 2));// First find on each face if they are crossed by some isovalues and// extract segments in a graph//首先在每个面上查找它们是否与某些等值交叉，并在图形中提取线段// 定义一个无向图，顶点属性值存储Point3using Segment_graph = boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, Point_3>;Segment_graph graph;//定义点到图中顶点索引的映射using Map_p2v = std::map<Point_3, Segment_graph::vertex_descriptor>;Map_p2v map_p2v;for (TIN::Face_handle vh : tin.finite_face_handles())for (double iv : isovalues)if (face_has_isovalue(vh, iv)){Segment_3 segment = isocontour_in_face(vh, iv);for (const Point_3& p : { segment.source(), segment.target() }){// Only insert end points of segments once to get a well connected graph// 只需插入一次线段的端点即可获得连接良好的图形Map_p2v::iterator iter;bool inserted;std::tie(iter, inserted) = map_p2v.insert(std::make_pair(p, Segment_graph::vertex_descriptor()));if (inserted){//向图中插入点，并设置值iter->second = boost::add_vertex(graph);graph[iter->second] = p;}}//向图中插入边boost::add_edge(map_p2v[segment.source()], map_p2v[segment.target()], graph);}// Split segments into polylines// 将线段分割为多段线std::vector<std::vector<Point_3> > polylines;Polylines_visitor<Segment_graph> visitor(graph, polylines);CGAL::split_graph_into_polylines(graph, visitor);std::cerr << polylines.size() << " polylines computed, with "<< map_p2v.size() << " vertices in total" << std::endl;// Output to WKT filestd::ofstream contour_ofile("contour.wkt");contour_ofile.precision(18);CGAL::IO::write_multi_linestring_WKT(contour_ofile, polylines);contour_ofile.close();// Construct constrained Delaunay triangulation with polylines as constraints// 以折线作为约束构造约束型 Delaunay 三角剖分// 计算平均间距double spacing = CGAL::compute_average_spacing<Concurrency_tag>(mesh.points(), 6);spacing *= 2;CTP ctp;for (const std::vector<Point_3>& poly : polylines)ctp.insert_constraint(poly.begin(), poly.end());// Simplification algorithm with limit on distance// 具有距离限制的简化算法PS::simplify(ctp, PS::Squared_distance_cost(), PS::Stop_above_cost_threshold(16 * spacing * spacing));polylines.clear();for (CTP::Constraint_id cid : ctp.constraints()){polylines.push_back(std::vector<Point_3>());polylines.back().reserve(ctp.vertices_in_constraint(cid).size());for (CTP::Vertex_handle vh : ctp.vertices_in_constraint(cid))polylines.back().push_back(vh->point());}std::size_t nb_vertices= std::accumulate(polylines.begin(), polylines.end(), std::size_t(0),[](std::size_t size, const std::vector<Point_3>& poly) -> std::size_t{ return size + poly.size(); });std::cerr << nb_vertices<< " vertices remaining after simplification ("<< 100. * (nb_vertices / double(map_p2v.size())) << "%)" << std::endl;// Output to WKT file// 输出到 WKT 文件std::ofstream simplified_ofile("simplified.wkt");simplified_ofile.precision(18);CGAL::IO::write_multi_linestring_WKT(simplified_ofile, polylines);simplified_ofile.close();return 0;
}