CGAL Mesh地形曲面提取等高线

CGAL 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 isovalues
    for (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
    //首先在每个面上查找它们是否与某些等值交叉，并在图形中提取线段
    // 定义一个无向图，顶点属性值存储Point3
    using 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 file
    std::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;
}

参考及相关链接

1. https://blog.csdn.net/mrbaolong/article/details/141679342?spm=1001.2014.3001.5501
2. https://doc.cgal.org/latest/Manual/tuto_gis.html