Mesh 网格曲面栅格化

使用CGAL将地形曲面网格体栅格化。

Code

#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 <CGAL/Polygon_mesh_processing/locate.h>

#include "include/Color_ramp.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;

// Triangulated Irregular Network
using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;
using Mesh = CGAL::Surface_mesh<Point_3>;

int main(){
   //读取网格曲面
    Mesh mesh;
    CGAL::IO::read_PLY("test.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;
    int width = std::ceil(bbox.x_span());
    int height = std::ceil(bbox.y_span());
    width = std::max(width, height);
    height = width;
    // 构建Delaunay_triangulation 三角网
    TIN tin (mesh.points().begin(), mesh.points().end());
    TIN::Face_handle location;

    // Use PPM format (Portable PixMap) for simplicity
    std::ofstream raster_ofile("raster.ppm", std::ios_base::binary);
    // PPM header
    raster_ofile << "P6" << std::endl // magic number
        << width << " " << height << std::endl // dimensions of the image
        << 255 << std::endl; // maximum color value
    // Use rainbow color ramp output
    Color_ramp color_ramp;
    //查询并插值，使用三角形重心插值
    for (std::size_t y = 0; y < height; ++y)
        for (std::size_t x = 0; x < width; ++x)
        {
            Point_3 query(bbox.xmin() + x * (bbox.xmax() - bbox.xmin()) / double(width),
                bbox.ymin() + (height - y) * (bbox.ymax() - bbox.ymin()) / double(height),
                0); // not relevant for location in 2D
            location = tin.locate(query, location);
            // Points outside the convex hull will be colored black
            std::array<unsigned char, 3> colors { 0, 0, 0 };
            // 查询到网格点所在三角形面
            if (!tin.is_infinite(location))
            {
                std::array<double, 3> barycentric_coordinates
                    = CGAL::Polygon_mesh_processing::barycentric_coordinates
                    (Point_2(location->vertex(0)->point().x(), location->vertex(0)->point().y()),
                        Point_2(location->vertex(1)->point().x(), location->vertex(1)->point().y()),
                        Point_2(location->vertex(2)->point().x(), location->vertex(2)->point().y()),
                        Point_2(query.x(), query.y()),
                        Kernel());
                  // 插值
                double height_at_query
                    = (barycentric_coordinates[0] * location->vertex(0)->point().z()
                        + barycentric_coordinates[1] * location->vertex(1)->point().z()
                        + barycentric_coordinates[2] * location->vertex(2)->point().z());
                // Color ramp generates a color depending on a value from 0 to 1
                // 高程值映射（拉伸或缩放）
                double height_ratio = (height_at_query - bbox.zmin()) / (bbox.zmax() - bbox.zmin());
                height_ratio = height_ratio > 1.0 ? 1.0 : height_ratio;
                colors = color_ramp.get(height_ratio);
            }
            //写入映射后的RGB
            raster_ofile.write(reinterpret_cast<char*>(&colors), 3);
        }
    raster_ofile.close();
    return 0;
}

构建

cmake_minimum_required(VERSION 3.1...3.23)
project(main)
find_package(CGAL REQUIRED)
create_single_source_cgal_program("main.cpp")

cmake -B build -S . -DCMAKE_TOOLCHAIN_FILE=D:\vcpkg\scripts\buildsystems\vcpkg.cmake
cmake --build build --config Debug

Result

网格体

栅格化符号化结果：红色->黄色->绿色->蓝色->白色（高程值递减）

反过程

参考及相关链接

1. https://blog.csdn.net/mrbaolong/article/details/137610742
2. https://doc.cgal.org/latest/Manual/tuto_gis.html
3. stl2png-Mesh转栅格
4. 高度栅格图转网格体Mesh