21.1 实际应用案例研究
本章将探讨CGAL在五个关键领域的实际应用:建筑信息模型(BIM)、3D打印、GIS地理信息系统、分子建模以及计算机视觉。每个案例都包含领域背景介绍、核心算法讲解以及可直接编译运行的代码示例。
相关章节
- 20.1 性能优化与调试 - 实际应用中的性能优化
- 16.1 AABB树 - 碰撞检测与空间查询
- 15.1 包围盒与PCA - 几何分析基础
- 12.1 点云处理 - 点云配准与重建
- 10.1 表面网格处理 - 网格修复与简化
21.1.1 案例一:建筑信息模型(BIM)
21.1.1.1 领域背景
建筑信息模型(Building Information Modeling, BIM)是建筑行业数字化转型的核心技术。BIM不仅包含建筑的几何信息,还集成了材料、成本、进度等多维度数据。CGAL在BIM中的应用主要集中在几何处理层面:
- 空间分析:房间体积计算、净高分析、视线分析
- 碰撞检测:管道与结构、设备与空间的干涉检查
- 施工模拟:4D进度模拟中的几何冲突预判
- 能源分析:建筑外表面面积、朝向分析用于能耗计算
21.1.1.2 核心算法
空间规划与碰撞检测
碰撞检测是BIM应用中最常见的几何问题。CGAL提供了多种空间查询结构:
// bim_collision_detection.cpp
// 建筑构件碰撞检测示例
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <iostream>
#include <vector>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef K::Triangle_3 Triangle_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
typedef CGAL::AABB_face_graph_triangle_primitive<Surface_mesh> Primitive;
typedef CGAL::AABB_tree<K::FT, Primitive> Tree;
typedef Tree::Point_and_primitive_id Point_and_primitive_id;
// 建筑构件类
struct BuildingElement {
std::string name;
std::string category; // "wall", "beam", "pipe", "duct"
Surface_mesh mesh;
CGAL::Bbox_3 bbox;
};
// 碰撞检测结果
struct Collision {
std::string element_a;
std::string element_b;
double penetration_depth;
Point_3 contact_point;
};
// 构建AABB树加速结构
Tree build_aabb_tree(const Surface_mesh& mesh) {
Tree tree(faces(mesh).first, faces(mesh).second, mesh);
tree.accelerate_distance_queries();
return tree;
}
// 检测两个构件是否碰撞
bool detect_collision(const BuildingElement& elem_a,
const BuildingElement& elem_b,
Collision& result) {
// 首先进行包围盒快速剔除
if (!CGAL::do_overlap(elem_a.bbox, elem_b.bbox)) {
return false;
}
// 使用CGAL的网格求交检测
bool intersecting = PMP::do_intersect(elem_a.mesh, elem_b.mesh);
if (intersecting) {
result.element_a = elem_a.name;
result.element_b = elem_b.name;
// 计算近似接触点(取两个包围盒中心的中点)
Point_3 center_a(
(elem_a.bbox.xmin() + elem_a.bbox.xmax()) / 2,
(elem_a.bbox.ymin() + elem_a.bbox.ymax()) / 2,
(elem_a.bbox.zmin() + elem_a.bbox.zmax()) / 2
);
Point_3 center_b(
(elem_b.bbox.xmin() + elem_b.bbox.xmax()) / 2,
(elem_b.bbox.ymin() + elem_b.bbox.ymax()) / 2,
(elem_b.bbox.zmin() + elem_b.bbox.zmax()) / 2
);
result.contact_point = Point_3(
(center_a.x() + center_b.x()) / 2,
(center_a.y() + center_b.y()) / 2,
(center_a.z() + center_b.z()) / 2
);
// 估算穿透深度
double dx = std::max(0.0,
std::min(elem_a.bbox.xmax() - elem_b.bbox.xmin(),
elem_b.bbox.xmax() - elem_a.bbox.xmin()));
double dy = std::max(0.0,
std::min(elem_a.bbox.ymax() - elem_b.bbox.ymin(),
elem_b.bbox.ymax() - elem_a.bbox.ymin()));
double dz = std::max(0.0,
std::min(elem_a.bbox.zmax() - elem_b.bbox.zmin(),
elem_b.bbox.zmax() - elem_a.bbox.zmin()));
result.penetration_depth = std::min({dx, dy, dz});
}
return intersecting;
}
// 批量碰撞检测
std::vector<Collision> batch_collision_detection(
const std::vector<BuildingElement>& elements) {
std::vector<Collision> collisions;
// 使用空间分割策略:按类别分组减少检测对数
std::map<std::string, std::vector<size_t>> category_groups;
for (size_t i = 0; i < elements.size(); ++i) {
category_groups[elements[i].category].push_back(i);
}
// 定义需要检测的类别对(某些类别之间不需要检测)
std::vector<std::pair<std::string, std::string>> check_pairs = {
{"pipe", "beam"},
{"pipe", "wall"},
{"duct", "beam"},
{"duct", "wall"},
{"equipment", "wall"}
};
for (const auto& pair : check_pairs) {
const auto& group_a = category_groups[pair.first];
const auto& group_b = category_groups[pair.second];
for (size_t idx_a : group_a) {
for (size_t idx_b : group_b) {
Collision col;
if (detect_collision(elements[idx_a], elements[idx_b], col)) {
collisions.push_back(col);
}
}
}
}
return collisions;
}
// 创建示例建筑构件(简化表示)
BuildingElement create_wall(const std::string& name,
double x, double y, double z,
double width, double height, double thickness) {
BuildingElement elem;
elem.name = name;
elem.category = "wall";
// 创建立方体网格表示墙体
std::vector<Point_3> points = {
Point_3(x, y, z),
Point_3(x + width, y, z),
Point_3(x + width, y + thickness, z),
Point_3(x, y + thickness, z),
Point_3(x, y, z + height),
Point_3(x + width, y, z + height),
Point_3(x + width, y + thickness, z + height),
Point_3(x, y + thickness, z + height)
};
for (const auto& p : points) {
elem.mesh.add_vertex(p);
}
// 添加面(简化,实际应用需要完整闭合网格)
elem.bbox = elem.mesh.bbox();
return elem;
}
int main() {
std::cout << "=== BIM碰撞检测系统 ===" << std::endl;
// 创建示例建筑模型
std::vector<BuildingElement> elements;
// 创建墙体
elements.push_back(create_wall("Wall_01", 0, 0, 0, 10, 3, 0.3));
elements.push_back(create_wall("Wall_02", 0, 5, 0, 10, 3, 0.3));
std::cout << "建筑模型包含 " << elements.size() << " 个构件" << std::endl;
std::cout << "碰撞检测功能演示完成" << std::endl;
return 0;
}房间体积与净高分析
// bim_space_analysis.cpp
// 建筑空间分析:体积计算与净高检查
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/Convex_hull_3/dual/halfspace_intersection_3.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <iostream>
#include <vector>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef K::Plane_3 Plane_3;
typedef CGAL::Polyhedron_3<K> Polyhedron;
typedef CGAL::Nef_polyhedron_3<K> Nef_polyhedron;
// 房间定义
struct Room {
std::string name;
std::string room_id;
std::vector<Point_3> floor_boundary; // 地面边界点
double floor_elevation;
double ceiling_elevation;
double required_min_height;
};
// 空间分析结果
struct SpaceAnalysisResult {
std::string room_name;
double volume;
double floor_area;
double min_ceiling_height;
double max_ceiling_height;
bool meets_height_requirement;
std::vector<Point_3> low_clearance_zones;
};
// 计算房间体积(使用棱柱法)
double calculate_room_volume(const Room& room) {
if (room.floor_boundary.size() < 3) {
return 0.0;
}
// 计算地面多边形面积
double floor_area = 0.0;
const auto& pts = room.floor_boundary;
size_t n = pts.size();
// 使用鞋带公式计算面积(假设点按顺序排列)
for (size_t i = 0; i < n; ++i) {
size_t j = (i + 1) % n;
floor_area += pts[i].x() * pts[j].y();
floor_area -= pts[j].x() * pts[i].y();
}
floor_area = std::abs(floor_area) / 2.0;
// 体积 = 底面积 × 高度
double height = room.ceiling_elevation - room.floor_elevation;
return floor_area * height;
}
// 净高分析
SpaceAnalysisResult analyze_room_clearance(const Room& room,
const std::vector<Plane_3>& obstacles) {
SpaceAnalysisResult result;
result.room_name = room.name;
// 基础几何计算
result.volume = calculate_room_volume(room);
// 计算地面面积
double floor_area = 0.0;
const auto& pts = room.floor_boundary;
for (size_t i = 0; i < pts.size(); ++i) {
size_t j = (i + 1) % pts.size();
floor_area += pts[i].x() * pts[j].y();
floor_area -= pts[j].x() * pts[i].y();
}
result.floor_area = std::abs(floor_area) / 2.0;
// 基础净高
double base_height = room.ceiling_elevation - room.floor_elevation;
result.max_ceiling_height = base_height;
result.min_ceiling_height = base_height;
// 检测障碍物导致的局部净高不足
for (const auto& plane : obstacles) {
// 简化处理:检查障碍物平面是否位于房间内
// 实际应用需要更复杂的几何计算
for (const auto& pt : room.floor_boundary) {
Point_3 test_point(pt.x(), pt.y(), room.floor_elevation + base_height / 2);
double dist = CGAL::to_double(plane.distance(test_point));
if (std::abs(dist) < base_height / 2) {
// 发现潜在干涉
double obstacle_height = base_height - std::abs(dist);
if (obstacle_height < result.min_ceiling_height) {
result.min_ceiling_height = obstacle_height;
}
if (obstacle_height < room.required_min_height) {
result.low_clearance_zones.push_back(test_point);
}
}
}
}
result.meets_height_requirement =
result.min_ceiling_height >= room.required_min_height;
return result;
}
// 生成空间分析报告
void generate_space_report(const std::vector<SpaceAnalysisResult>& results) {
std::cout << "\n========== 空间分析报告 ==========" << std::endl;
std::cout << std::left << std::setw(15) << "房间名称"
<< std::setw(12) << "体积(m³)"
<< std::setw(12) << "面积(m²)"
<< std::setw(12) << "最小净高"
<< std::setw(12) << "最大净高"
<< std::setw(10) << "合规" << std::endl;
std::cout << std::string(75, '-') << std::endl;
for (const auto& r : results) {
std::cout << std::left << std::setw(15) << r.room_name
<< std::setw(12) << std::fixed << std::setprecision(2) << r.volume
<< std::setw(12) << r.floor_area
<< std::setw(12) << r.min_ceiling_height
<< std::setw(12) << r.max_ceiling_height
<< std::setw(10) << (r.meets_height_requirement ? "是" : "否")
<< std::endl;
if (!r.low_clearance_zones.empty()) {
std::cout << " [警告] 发现 " << r.low_clearance_zones.size()
<< " 处净高不足区域" << std::endl;
}
}
}
int main() {
std::cout << "=== BIM空间分析系统 ===" << std::endl;
// 创建示例房间
std::vector<Room> rooms;
Room room1;
room1.name = "会议室";
room1.room_id = "R001";
room1.floor_elevation = 0.0;
room1.ceiling_elevation = 3.0;
room1.required_min_height = 2.4;
room1.floor_boundary = {
Point_3(0, 0, 0),
Point_3(8, 0, 0),
Point_3(8, 6, 0),
Point_3(0, 6, 0)
};
rooms.push_back(room1);
Room room2;
room2.name = "走廊";
room2.room_id = "R002";
room2.floor_elevation = 0.0;
room2.ceiling_elevation = 2.8;
room2.required_min_height = 2.2;
room2.floor_boundary = {
Point_3(8, 0, 0),
Point_3(12, 0, 0),
Point_3(12, 2, 0),
Point_3(8, 2, 0)
};
rooms.push_back(room2);
// 执行分析
std::vector<SpaceAnalysisResult> results;
std::vector<Plane_3> obstacles; // 示例:无障碍物
for (const auto& room : rooms) {
results.push_back(analyze_room_clearance(room, obstacles));
}
// 生成报告
generate_space_report(results);
return 0;
}21.1.1.3 扩展资源
- IFC标准: Industry Foundation Classes, BIM数据交换的国际标准
- CGAL与IFC: 使用IfcOpenShell库将IFC模型转换为CGAL可处理的网格
- 相关CGAL包: Polygon Mesh Processing, AABB Tree, Nef Polyhedra
21.1.2 案例二:3D打印流程
21.1.2.1 领域背景
3D打印(增材制造)是将数字模型转化为实体对象的技术。CGAL在3D打印流程中的应用贯穿始终:
- 模型修复:修复非流形几何、补洞、确保水密性
- 支撑生成:自动生成可移除的支撑结构
- 切片处理:将3D模型转换为2D轮廓层
- 路径规划:生成最优的打印头运动路径
21.1.2.2 核心算法
模型修复与验证
// print_model_repair.cpp
// 3D打印模型修复与验证
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <iostream>
#include <vector>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
// 模型质量评估结果
struct ModelQualityReport {
bool is_manifold;
bool is_watertight;
bool has_self_intersections;
int num_connected_components;
int num_boundary_edges;
int num_degenerate_faces;
double volume;
double surface_area;
std::vector<std::string> issues;
};
// 检查模型是否适合3D打印
ModelQualityReport validate_print_model(const Surface_mesh& mesh) {
ModelQualityReport report;
// 1. 检查是否为流形
report.is_manifold = mesh.is_valid() && PMP::is_valid_polygon_mesh(mesh);
if (!report.is_manifold) {
report.issues.push_back("模型不是有效的流形网格");
}
// 2. 检查是否水密(无边界边)
std::vector<Surface_mesh::Halfedge_index> border_edges;
PMP::border_halfedges(faces(mesh), mesh, std::back_inserter(border_edges));
report.num_boundary_edges = border_edges.size();
report.is_watertight = (report.num_boundary_edges == 0);
if (!report.is_watertight) {
report.issues.push_back("模型存在 " +
std::to_string(report.num_boundary_edges) + " 条边界边,不是水密模型");
}
// 3. 检查自相交
std::vector<std::pair<Surface_mesh::Face_index, Surface_mesh::Face_index>>
intersecting_faces;
PMP::self_intersections(mesh, std::back_inserter(intersecting_faces));
report.has_self_intersections = !intersecting_faces.empty();
if (report.has_self_intersections) {
report.issues.push_back("模型存在 " +
std::to_string(intersecting_faces.size()) + " 对自相交面");
}
// 4. 检查连通分量
std::vector<Surface_mesh::Face_index> face_cc;
report.num_connected_components =
PMP::connected_components(mesh, boost::make_function_output_iterator(
[&face_cc](const std::pair<Surface_mesh::Face_index, int>& p) {
face_cc.push_back(p.first);
}));
if (report.num_connected_components > 1) {
report.issues.push_back("模型包含 " +
std::to_string(report.num_connected_components) + " 个不连通部分");
}
// 5. 计算几何属性
if (report.is_watertight && !report.has_self_intersections) {
report.volume = CGAL::to_double(PMP::volume(mesh));
report.surface_area = CGAL::to_double(PMP::area(mesh));
}
return report;
}
// 修复常见问题
bool repair_model(Surface_mesh& mesh, ModelQualityReport& report) {
bool repaired = false;
// 1. 修复退化面
int degenerate_removed = PMP::remove_degenerate_faces(mesh);
if (degenerate_removed > 0) {
report.issues.push_back("移除了 " + std::to_string(degenerate_removed) + " 个退化面");
repaired = true;
}
// 2. 缝合接近的边界
double tolerance = 0.001; // 1mm容差
int stitched = PMP::stitch_borders(mesh, tolerance);
if (stitched > 0) {
report.issues.push_back("缝合了 " + std::to_string(stitched) + " 对边界边");
repaired = true;
}
// 3. 修复孔洞(简化版,实际应用需要更复杂的算法)
// PMP::triangulate_hole() 可用于填充小孔
// 4. 统一面朝向
if (PMP::is_outward_oriented(mesh)) {
PMP::reverse_face_orientations(mesh);
report.issues.push_back("统一了面片朝向");
repaired = true;
}
return repaired;
}
// 生成打印可行性报告
void print_validation_report(const ModelQualityReport& report) {
std::cout << "\n========== 3D打印模型验证报告 ==========" << std::endl;
std::cout << "流形网格: " << (report.is_manifold ? "是" : "否") << std::endl;
std::cout << "水密性: " << (report.is_watertight ? "是" : "否") << std::endl;
std::cout << "自相交: " << (report.has_self_intersections ? "有" : "无") << std::endl;
std::cout << "连通分量数: " << report.num_connected_components << std::endl;
std::cout << "边界边数: " << report.num_boundary_edges << std::endl;
std::cout << "体积: " << std::fixed << std::setprecision(4)
<< report.volume << " mm³" << std::endl;
std::cout << "表面积: " << report.surface_area << " mm²" << std::endl;
if (!report.issues.empty()) {
std::cout << "\n发现问题:" << std::endl;
for (const auto& issue : report.issues) {
std::cout << " - " << issue << std::endl;
}
}
// 判断是否可打印
bool printable = report.is_manifold && report.is_watertight &&
!report.has_self_intersections;
std::cout << "\n打印可行性: " << (printable ? "可打印" : "需修复") << std::endl;
}
// 创建测试用简单立方体网格
Surface_mesh create_test_cube(double size) {
Surface_mesh mesh;
// 创建立方体顶点
auto v0 = mesh.add_vertex(Point_3(0, 0, 0));
auto v1 = mesh.add_vertex(Point_3(size, 0, 0));
auto v2 = mesh.add_vertex(Point_3(size, size, 0));
auto v3 = mesh.add_vertex(Point_3(0, size, 0));
auto v4 = mesh.add_vertex(Point_3(0, 0, size));
auto v5 = mesh.add_vertex(Point_3(size, 0, size));
auto v6 = mesh.add_vertex(Point_3(size, size, size));
auto v7 = mesh.add_vertex(Point_3(0, size, size));
// 添加6个面(每个面2个三角形)
mesh.add_face(v0, v1, v2); mesh.add_face(v0, v2, v3); // 底面
mesh.add_face(v4, v6, v5); mesh.add_face(v4, v7, v6); // 顶面
mesh.add_face(v0, v4, v5); mesh.add_face(v0, v5, v1); // 前面
mesh.add_face(v2, v6, v7); mesh.add_face(v2, v7, v3); // 后面
mesh.add_face(v0, v3, v7); mesh.add_face(v0, v7, v4); // 左面
mesh.add_face(v1, v5, v6); mesh.add_face(v1, v6, v2); // 右面
return mesh;
}
int main() {
std::cout << "=== 3D打印模型验证工具 ===" << std::endl;
// 创建测试模型
Surface_mesh mesh = create_test_cube(10.0);
std::cout << "模型顶点数: " << mesh.number_of_vertices() << std::endl;
std::cout << "模型面数: " << mesh.number_of_faces() << std::endl;
// 验证模型
ModelQualityReport report = validate_print_model(mesh);
print_validation_report(report);
return 0;
}切片与路径规划
// print_slicing.cpp
// 3D打印切片处理
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/clip.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/create_straight_skeleton_2.h>
#include <CGAL/Straight_skeleton_2.h>
#include <iostream>
#include <vector>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef K::Point_3 Point_3;
typedef K::Plane_3 Plane_3;
typedef K::Segment_3 Segment_3;
typedef CGAL::Polygon_2<K> Polygon_2;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
// 切片层定义
struct SliceLayer {
double z_height;
std::vector<Polygon_2> contours; // 外轮廓
std::vector<Polygon_2> holes; // 内孔洞
std::vector<Point_2> infill_paths; // 填充路径
};
// 切片参数
struct SlicingParameters {
double layer_height; // 层厚(mm)
double first_layer_height; // 首层高度
double infill_density; // 填充密度(0-1)
int wall_count; // 壁厚(圈数)
double nozzle_diameter; // 喷嘴直径
};
// 计算切片平面
std::vector<Plane_3> compute_slicing_planes(double min_z, double max_z,
const SlicingParameters& params) {
std::vector<Plane_3> planes;
double current_z = min_z + params.first_layer_height;
while (current_z <= max_z) {
// 创建水平切片平面
planes.push_back(Plane_3(Point_3(0, 0, current_z), K::Vector_3(0, 0, 1)));
current_z += params.layer_height;
}
return planes;
}
// 执行切片(简化版)
SliceLayer slice_at_height(const Surface_mesh& mesh, double z_height) {
SliceLayer layer;
layer.z_height = z_height;
// 创建切片平面
Plane_3 plane(Point_3(0, 0, z_height), K::Vector_3(0, 0, 1));
// 收集与平面相交的边
std::vector<Segment_3> intersections;
for (auto edge : mesh.edges()) {
auto h = mesh.halfedge(edge);
auto v1 = mesh.source(h);
auto v2 = mesh.target(h);
Point_3 p1 = mesh.point(v1);
Point_3 p2 = mesh.point(v2);
// 检查边是否跨越平面
double d1 = p1.z() - z_height;
double d2 = p2.z() - z_height;
if (d1 * d2 < 0) { // 符号不同,说明跨越平面
// 线性插值计算交点
double t = std::abs(d1) / (std::abs(d1) + std::abs(d2));
Point_3 intersection(
p1.x() + t * (p2.x() - p1.x()),
p1.y() + t * (p2.y() - p1.y()),
z_height
);
// 简化为2D点存储
// 实际应用需要构建完整的多边形轮廓
}
}
// 注意:完整实现需要:
// 1. 将交点连接成闭合轮廓
// 2. 处理多个轮廓和孔洞
// 3. 使用CGAL的2D多边形操作
return layer;
}
// 生成填充路径(直线填充)
std::vector<Point_2> generate_infill(const Polygon_2& contour,
const std::vector<Polygon_2>& holes,
double density,
double line_spacing) {
std::vector<Point_2> paths;
// 计算边界框
double min_x = contour.bbox().xmin();
double max_x = contour.bbox().xmax();
double min_y = contour.bbox().ymin();
double max_y = contour.bbox().ymax();
// 生成平行线填充
double spacing = line_spacing / density;
for (double y = min_y; y <= max_y; y += spacing) {
// 创建水平扫描线
// 计算与轮廓的交点
// 排序交点并连接成路径段
// 简化的伪代码:
// intersections = compute_intersections(contour, y)
// for i = 0 to intersections.size()-1 step 2:
// paths.push_back(intersections[i])
// paths.push_back(intersections[i+1])
}
return paths;
}
// 生成G-code(简化版)
void generate_gcode(const std::vector<SliceLayer>& layers,
const SlicingParameters& params,
std::ostream& output) {
output << "; Generated by CGAL Slicer" << std::endl;
output << "; Layer height: " << params.layer_height << std::endl;
output << std::fixed << std::setprecision(3);
// 初始化G-code
output << "G28 ; Home all axes" << std::endl;
output << "G1 Z" << params.first_layer_height << " F300" << std::endl;
output << "M104 S200 ; Set hotend temp" << std::endl;
output << "M140 S60 ; Set bed temp" << std::endl;
for (const auto& layer : layers) {
output << "\n; Layer at Z=" << layer.z_height << std::endl;
output << "G1 Z" << layer.z_height << " F300" << std::endl;
// 打印外轮廓
for (const auto& contour : layer.contours) {
output << "; Outer contour" << std::endl;
for (auto it = contour.vertices_begin();
it != contour.vertices_end(); ++it) {
output << "G1 X" << it->x() << " Y" << it->y() << " F1500" << std::endl;
}
}
// 打印填充
output << "; Infill" << std::endl;
for (size_t i = 0; i < layer.infill_paths.size(); i += 2) {
if (i + 1 < layer.infill_paths.size()) {
output << "G1 X" << layer.infill_paths[i].x()
<< " Y" << layer.infill_paths[i].y() << std::endl;
output << "G1 X" << layer.infill_paths[i+1].x()
<< " Y" << layer.infill_paths[i+1].y() << std::endl;
}
}
}
// 结束G-code
output << "\nG28 ; Home" << std::endl;
output << "M104 S0 ; Turn off hotend" << std::endl;
output << "M140 S0 ; Turn off bed" << std::endl;
}
int main() {
std::cout << "=== 3D打印切片系统 ===" << std::endl;
SlicingParameters params;
params.layer_height = 0.2;
params.first_layer_height = 0.2;
params.infill_density = 0.2;
params.wall_count = 2;
params.nozzle_diameter = 0.4;
std::cout << "切片参数:" << std::endl;
std::cout << " 层厚: " << params.layer_height << "mm" << std::endl;
std::cout << " 填充密度: " << params.infill_density * 100 << "%" << std::endl;
std::cout << " 壁厚: " << params.wall_count << "圈" << std::endl;
// 计算切片平面
auto planes = compute_slicing_planes(0, 10, params);
std::cout << "\n生成 " << planes.size() << " 个切片层" << std::endl;
return 0;
}21.1.2.3 扩展资源
- STL文件格式: 3D打印的标准文件格式,CGAL提供STL读写支持
- 支撑生成算法: 基于悬垂角检测的自动支撑生成
- 相关CGAL包: Polygon Mesh Processing, 2D Polygon, Straight Skeleton
21.1.3 案例三:GIS地理信息系统
21.1.3.1 领域背景
地理信息系统(GIS)用于采集、存储、分析和展示地理空间数据。CGAL在GIS中的应用包括:
- 地形分析:数字高程模型(DEM)处理、坡度坡向计算
- 可视性分析:视域分析、视线通视判断
- 水文分析:流域提取、河网生成
- 空间索引:大规模地形数据的快速查询
21.1.3.2 核心算法
地形分析与TIN生成
// gis_terrain_analysis.cpp
// GIS地形分析:DEM处理和TIN生成
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Triangulation_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/interpolation_functions.h>
#include <iostream>
#include <vector>
#include <fstream>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef K::Point_3 Point_3;
typedef K::Iso_rectangle_2 Iso_rectangle_2;
// 顶点信息:高程值
struct VertexInfo {
double elevation;
double slope; // 坡度(度)
double aspect; // 坡向(度,北为0)
};
// 面片信息
struct FaceInfo {
double avg_elevation;
bool is_valid;
};
typedef CGAL::Triangulation_vertex_base_with_info_2<VertexInfo, K> Vb;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo, K> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds> Delaunay_triangulation;
// 地形点云数据
struct TerrainPoint {
double x, y, z;
};
// 地形分析结果
struct TerrainAnalysis {
double min_elevation;
double max_elevation;
double mean_elevation;
double total_area;
std::vector<double> slope_distribution; // 坡度分布直方图
};
// 从不规则点云构建TIN
Delaunay_triangulation build_tin(const std::vector<TerrainPoint>& points) {
Delaunay_triangulation dt;
for (const auto& p : points) {
VertexInfo info;
info.elevation = p.z;
info.slope = 0.0;
info.aspect = 0.0;
auto vh = dt.insert(Point_2(p.x, p.y));
vh->info() = info;
}
return dt;
}
// 计算每个顶点的坡度和坡向
void compute_slope_aspect(Delaunay_triangulation& dt) {
for (auto vit = dt.finite_vertices_begin();
vit != dt.finite_vertices_end(); ++vit) {
// 收集相邻三角形
std::vector<Point_3> neighbors;
auto circulator = dt.incident_vertices(vit);
auto done = circulator;
do {
if (!dt.is_infinite(circulator)) {
Point_2 p = circulator->point();
double z = circulator->info().elevation;
neighbors.push_back(Point_3(p.x(), p.y(), z));
}
} while (++circulator != done);
if (neighbors.size() < 3) continue;
// 使用最小二乘法拟合局部平面
// 简化的坡度计算:使用相邻点的高差
double max_slope = 0.0;
double sum_aspect_x = 0.0;
double sum_aspect_y = 0.0;
int count = 0;
Point_2 center = vit->point();
double center_z = vit->info().elevation;
for (const auto& n : neighbors) {
double dx = n.x() - center.x();
double dy = n.y() - center.y();
double dz = n.z() - center_z;
double dist = std::sqrt(dx*dx + dy*dy);
if (dist > 0) {
double slope = std::atan2(std::abs(dz), dist) * 180.0 / M_PI;
max_slope = std::max(max_slope, slope);
// 坡向计算
if (dz != 0) {
sum_aspect_x += dx / dist * (dz > 0 ? 1 : -1);
sum_aspect_y += dy / dist * (dz > 0 ? 1 : -1);
count++;
}
}
}
vit->info().slope = max_slope;
if (count > 0) {
vit->info().aspect = std::atan2(sum_aspect_y / count,
sum_aspect_x / count) * 180.0 / M_PI;
if (vit->info().aspect < 0) vit->info().aspect += 360.0;
}
}
}
// 自然邻域插值:估算任意点的高程
double interpolate_elevation(const Delaunay_triangulation& dt,
const Point_2& query_point) {
// 检查点是否在凸包内
if (dt.dimension() != 2) return 0.0;
typedef std::vector<std::pair<Point_2, K::FT> > Point_coordinate_vector;
Point_coordinate_vector coords;
K::FT norm = CGAL::natural_neighbor_coordinates_2(dt, query_point,
std::back_inserter(coords)).second;
if (norm == 0) return 0.0;
// 加权插值
double elevation = 0.0;
for (const auto& coord : coords) {
auto vh = dt.nearest_vertex(coord.first);
elevation += CGAL::to_double(coord.second) * vh->info().elevation / CGAL::to_double(norm);
}
return elevation;
}
// 生成等高线(简化版)
std::vector<std::vector<Point_2>> generate_contours(const Delaunay_triangulation& dt,
double contour_interval) {
std::vector<std::vector<Point_2>> contours;
// 确定高程范围
double min_z = std::numeric_limits<double>::max();
double max_z = std::numeric_limits<double>::lowest();
for (auto vit = dt.finite_vertices_begin();
vit != dt.finite_vertices_end(); ++vit) {
min_z = std::min(min_z, vit->info().elevation);
max_z = std::max(max_z, vit->info().elevation);
}
// 对每个等高线高程值
for (double z = std::ceil(min_z / contour_interval) * contour_interval;
z <= max_z; z += contour_interval) {
std::vector<Point_2> contour_line;
// 遍历所有边,寻找与等高线相交的边
for (auto eit = dt.finite_edges_begin();
eit != dt.finite_edges_end(); ++eit) {
auto v1 = eit->first->vertex((eit->second + 1) % 3);
auto v2 = eit->first->vertex((eit->second + 2) % 3);
double z1 = v1->info().elevation;
double z2 = v2->info().elevation;
// 检查边是否跨越等高线
if ((z1 <= z && z2 >= z) || (z2 <= z && z1 >= z)) {
if (std::abs(z2 - z1) > 1e-6) {
double t = (z - z1) / (z2 - z1);
Point_2 p1 = v1->point();
Point_2 p2 = v2->point();
Point_2 intersection(
p1.x() + t * (p2.x() - p1.x()),
p1.y() + t * (p2.y() - p1.y())
);
contour_line.push_back(intersection);
}
}
}
if (!contour_line.empty()) {
contours.push_back(contour_line);
}
}
return contours;
}
// 地形统计分析
TerrainAnalysis analyze_terrain(const Delaunay_triangulation& dt) {
TerrainAnalysis analysis;
analysis.min_elevation = std::numeric_limits<double>::max();
analysis.max_elevation = std::numeric_limits<double>::lowest();
double sum_elevation = 0.0;
int vertex_count = 0;
// 高程统计
for (auto vit = dt.finite_vertices_begin();
vit != dt.finite_vertices_end(); ++vit) {
double z = vit->info().elevation;
analysis.min_elevation = std::min(analysis.min_elevation, z);
analysis.max_elevation = std::max(analysis.max_elevation, z);
sum_elevation += z;
vertex_count++;
}
analysis.mean_elevation = sum_elevation / vertex_count;
// 坡度分布统计
analysis.slope_distribution.resize(9, 0); // 0-5, 5-10, ..., 35-40, >40
for (auto vit = dt.finite_vertices_begin();
vit != dt.finite_vertices_end(); ++vit) {
double slope = vit->info().slope;
int bin = std::min(static_cast<int>(slope / 5), 8);
analysis.slope_distribution[bin]++;
}
// 计算总面积(所有三角形面积之和)
analysis.total_area = 0.0;
for (auto fit = dt.finite_faces_begin();
fit != dt.finite_faces_end(); ++fit) {
auto v1 = fit->vertex(0)->point();
auto v2 = fit->vertex(1)->point();
auto v3 = fit->vertex(2)->point();
// 2D三角形面积
double area = std::abs(
(v2.x() - v1.x()) * (v3.y() - v1.y()) -
(v3.x() - v1.x()) * (v2.y() - v1.y())
) / 2.0;
analysis.total_area += area;
}
return analysis;
}
int main() {
std::cout << "=== GIS地形分析系统 ===" << std::endl;
// 创建示例地形点云
std::vector<TerrainPoint> points;
// 生成简单的山丘地形
for (int i = 0; i <= 20; ++i) {
for (int j = 0; j <= 20; ++j) {
double x = i * 10.0;
double y = j * 10.0;
// 高斯山丘
double dx = x - 100.0;
double dy = y - 100.0;
double z = 50.0 * std::exp(-(dx*dx + dy*dy) / 2000.0);
points.push_back({x, y, z});
}
}
std::cout << "输入点数: " << points.size() << std::endl;
// 构建TIN
Delaunay_triangulation dt = build_tin(points);
std::cout << "TIN构建完成" << std::endl;
std::cout << "三角形数: " << dt.number_of_faces() << std::endl;
// 计算坡度和坡向
compute_slope_aspect(dt);
std::cout << "坡度/坡向计算完成" << std::endl;
// 地形分析
TerrainAnalysis analysis = analyze_terrain(dt);
std::cout << "\n地形统计:" << std::endl;
std::cout << " 最小高程: " << analysis.min_elevation << "m" << std::endl;
std::cout << " 最大高程: " << analysis.max_elevation << "m" << std::endl;
std::cout << " 平均高程: " << analysis.mean_elevation << "m" << std::endl;
std::cout << " 总面积: " << analysis.total_area << "m²" << std::endl;
// 插值测试
Point_2 query(100, 100);
double interpolated_z = interpolate_elevation(dt, query);
std::cout << "\n插值测试 (100, 100): " << interpolated_z << "m" << std::endl;
return 0;
}可见性分析
// gis_visibility.cpp
// GIS可见性分析:视域和视线分析
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Arrangement_2.h>
#include <CGAL/Arr_segment_traits_2.h>
#include <CGAL/Triangular_expansion_visibility_2.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>
#include <iostream>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef K::Point_3 Point_3;
typedef K::Segment_2 Segment_2;
typedef K::Triangle_3 Triangle_3;
typedef K::Ray_3 Ray_3;
typedef CGAL::Arr_segment_traits_2<K> Traits_2;
typedef CGAL::Arrangement_2<Traits_2> Arrangement_2;
// 视域分析参数
struct ViewshedParams {
Point_3 observer_position;
double observer_height; // 观测者眼睛高度
double max_distance; // 最大可视距离
double view_angle; // 水平视场角(度)
double vertical_angle; // 垂直视场角(度)
};
// 视线分析结果
struct LineOfSightResult {
bool is_visible;
double distance_to_target;
Point_3 intersection_point;
double elevation_difference;
};
// 基于TIN的视线分析
class TerrainVisibilityAnalyzer {
public:
// 从TIN构建AABB树加速结构
void build_from_tin(const std::vector<Point_3>& tin_vertices,
const std::vector<std::array<int, 3>>& tin_faces) {
triangles_.clear();
for (const auto& f : tin_faces) {
triangles_.push_back(Triangle_3(
tin_vertices[f[0]],
tin_vertices[f[1]],
tin_vertices[f[2]]
));
}
// 构建AABB树
// 实际应用需要完整的AABB树实现
}
// 两点间视线分析
LineOfSightResult analyze_line_of_sight(const Point_3& observer,
const Point_3& target) const {
LineOfSightResult result;
result.distance_to_target = std::sqrt(
CGAL::to_double((target - observer).squared_length())
);
// 创建视线射线
Ray_3 line_of_sight(observer, target);
// 简化的视线分析
// 实际应用需要遍历地形三角形进行求交测试
// 检查视线是否被地形遮挡
result.is_visible = true;
result.elevation_difference = target.z() - observer.z();
return result;
}
// 视域分析(计算从观测点可见的所有区域)
std::vector<Point_3> compute_viewshed(const ViewshedParams& params,
int num_rays) const {
std::vector<Point_3> visible_points;
// 在水平面上均匀发射射线
for (int i = 0; i < num_rays; ++i) {
double angle = 2.0 * M_PI * i / num_rays;
double dx = std::cos(angle) * params.max_distance;
double dy = std::sin(angle) * params.max_distance;
Point_3 ray_end(
params.observer_position.x() + dx,
params.observer_position.y() + dy,
params.observer_position.z()
);
// 沿射线采样检查可见性
// 实际应用需要更精细的采样和地形求交
LineOfSightResult los = analyze_line_of_sight(
params.observer_position, ray_end
);
if (los.is_visible) {
visible_points.push_back(ray_end);
}
}
return visible_points;
}
// 计算可视域面积
double compute_visible_area(const ViewshedParams& params,
int num_samples) const {
double visible_area = 0.0;
double sample_spacing = params.max_distance / std::sqrt(num_samples);
for (int i = 0; i < num_samples; ++i) {
// 随机或均匀采样
double r = params.max_distance * std::sqrt(static_cast<double>(i) / num_samples);
double theta = 2.0 * M_PI * i / num_samples * 10; // 黄金角分布
Point_3 sample_point(
params.observer_position.x() + r * std::cos(theta),
params.observer_position.y() + r * std::sin(theta),
0 // 需要查询地形高程
);
// 检查可见性并累加面积
LineOfSightResult los = analyze_line_of_sight(
params.observer_position, sample_point
);
if (los.is_visible) {
visible_area += sample_spacing * sample_spacing;
}
}
return visible_area;
}
private:
std::vector<Triangle_3> triangles_;
};
// 2D多边形可见性(用于建筑或室内空间)
std::vector<Point_2> compute_polygon_visibility(const std::vector<Point_2>& polygon,
const Point_2& observer) {
std::vector<Point_2> visible_region;
// 使用CGAL的可见性算法
// 1. 构建多边形的Arrangement
// 2. 使用Triangular_expansion_visibility_2计算可见区域
// 简化的实现:返回多边形顶点作为可见区域近似
Arrangement_2 arr;
// 添加多边形边到arrangement
for (size_t i = 0; i < polygon.size(); ++i) {
size_t j = (i + 1) % polygon.size();
CGAL::insert(arr, Segment_2(polygon[i], polygon[j]));
}
// 找到观测点所在的面
auto face = arr.point_location(observer);
// 计算可见性(需要完整的可见性算法实现)
// CGAL::Triangular_expansion_visibility_2<Arrangement_2> visibility(arr);
// visibility.compute_visibility(observer, face, visible_region);
return visible_region;
}
// 生成视域分析报告
void generate_viewshed_report(const ViewshedParams& params,
const std::vector<Point_3>& visible_points) {
std::cout << "\n========== 视域分析报告 ==========" << std::endl;
std::cout << "观测点位置: (" << params.observer_position.x() << ", "
<< params.observer_position.y() << ", "
<< params.observer_position.z() << ")" << std::endl;
std::cout << "观测高度: " << params.observer_height << "m" << std::endl;
std::cout << "最大可视距离: " << params.max_distance << "m" << std::endl;
std::cout << "可见点数: " << visible_points.size() << std::endl;
// 计算覆盖角度
if (visible_points.size() >= 2) {
double min_angle = 2 * M_PI;
double max_angle = 0;
for (const auto& p : visible_points) {
double dx = p.x() - params.observer_position.x();
double dy = p.y() - params.observer_position.y();
double angle = std::atan2(dy, dx);
if (angle < 0) angle += 2 * M_PI;
min_angle = std::min(min_angle, angle);
max_angle = std::max(max_angle, angle);
}
double coverage = (max_angle - min_angle) * 180.0 / M_PI;
std::cout << "水平覆盖角度: " << coverage << "°" << std::endl;
}
}
int main() {
std::cout << "=== GIS可见性分析系统 ===" << std::endl;
// 设置视域分析参数
ViewshedParams params;
params.observer_position = Point_3(100, 100, 50);
params.observer_height = 1.7; // 人眼高度
params.max_distance = 500.0;
params.view_angle = 360.0;
params.vertical_angle = 30.0;
TerrainVisibilityAnalyzer analyzer;
// 执行视域分析
auto visible_points = analyzer.compute_viewshed(params, 360);
// 生成报告
generate_viewshed_report(params, visible_points);
// 视线分析示例
Point_3 target(200, 150, 45);
auto los_result = analyzer.analyze_line_of_sight(params.observer_position, target);
std::cout << "\n视线分析 (到目标点):" << std::endl;
std::cout << " 距离: " << los_result.distance_to_target << "m" << std::endl;
std::cout << " 可见性: " << (los_result.is_visible ? "可见" : "遮挡") << std::endl;
std::cout << " 高差: " << los_result.elevation_difference << "m" << std::endl;
return 0;
}21.1.3.3 扩展资源
- GDAL: 地理空间数据抽象库,与CGAL配合使用处理栅格数据
- PROJ: 坐标参考系统转换库
- 相关CGAL包: 2D Triangulation, Interpolation, Visibility
21.1.4 案例四:分子建模与生物信息学
21.1.4.1 领域背景
分子建模是计算化学和生物信息学的核心领域。CGAL在分子建模中的应用包括:
- 分子表面计算:溶剂可及表面(SAS)、分子表面(SES/VDW)
- 口袋检测:蛋白质活性位点识别
- 分子对接:配体-受体相互作用几何分析
- 分子动力学:邻居列表构建、空间索引
21.1.4.2 核心算法
分子表面计算
// molecular_surface.cpp
// 分子表面计算:VDW表面和溶剂可及表面
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Regular_triangulation_cell_base_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/alpha_wrap_3.h>
#include <iostream>
#include <vector>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef K::Weighted_point_3 Weighted_point_3;
typedef K::Sphere_3 Sphere_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
// 原子定义
struct Atom {
std::string element; // 元素符号
Point_3 position; // 原子坐标(埃)
double radius; // VDW半径(埃)
double charge; // 部分电荷
int residue_id; // 残基编号
std::string residue_name;
};
// 分子表面类型
enum SurfaceType {
VDW_SURFACE, // 范德华表面
SAS_SURFACE, // 溶剂可及表面
SES_SURFACE // 溶剂排除表面(Connolly表面)
};
// 表面属性
struct SurfaceProperties {
double total_area;
double volume;
double avg_curvature;
int num_pockets; // 口袋数量
std::vector<double> pocket_volumes;
};
// 分子类
class Molecule {
public:
void add_atom(const Atom& atom) {
atoms_.push_back(atom);
}
const std::vector<Atom>& get_atoms() const { return atoms_; }
// 计算VDW表面(使用Alpha Shape或Union of Spheres)
Surface_mesh compute_vdw_surface(double probe_radius = 0.0) const {
Surface_mesh surface;
// 方法1:使用加权Alpha Shape
// 将原子表示为加权点(球体)
std::vector<Weighted_point_3> weighted_points;
for (const auto& atom : atoms_) {
double weight = -(atom.radius + probe_radius) * (atom.radius + probe_radius);
weighted_points.push_back(Weighted_point_3(atom.position, weight));
}
// 使用CGAL的alpha_wrap_3生成包裹网格
// 这是一个近似方法,实际分子表面需要更精确的算法
if (!weighted_points.empty()) {
// 参数:相对Alpha值和偏移量
double relative_alpha = 10.0; // 需要根据分子大小调整
double offset = 0.5;
CGAL::alpha_wrap_3(weighted_points, relative_alpha, offset, surface);
}
return surface;
}
// 计算溶剂可及表面
Surface_mesh compute_sas_surface(double solvent_radius = 1.4) const {
// SAS = VDW表面 + 溶剂半径
return compute_vdw_surface(solvent_radius);
}
// 计算分子表面积
double compute_surface_area(const Surface_mesh& surface) const {
return CGAL::to_double(PMP::area(surface));
}
// 计算分子体积
double compute_volume(const Surface_mesh& surface) const {
if (PMP::is_closed(surface)) {
return CGAL::to_double(PMP::volume(surface));
}
return 0.0;
}
// 计算质心
Point_3 compute_centroid() const {
double cx = 0, cy = 0, cz = 0;
for (const auto& atom : atoms_) {
cx += atom.position.x();
cy += atom.position.y();
cz += atom.position.z();
}
double n = atoms_.size();
return Point_3(cx/n, cy/n, cz/n);
}
// 获取包围盒
CGAL::Bbox_3 get_bounding_box() const {
if (atoms_.empty()) return CGAL::Bbox_3();
double min_x = atoms_[0].position.x() - atoms_[0].radius;
double max_x = atoms_[0].position.x() + atoms_[0].radius;
double min_y = atoms_[0].position.y() - atoms_[0].radius;
double max_y = atoms_[0].position.y() + atoms_[0].radius;
double min_z = atoms_[0].position.z() - atoms_[0].radius;
double max_z = atoms_[0].position.z() + atoms_[0].radius;
for (const auto& atom : atoms_) {
min_x = std::min(min_x, atom.position.x() - atom.radius);
max_x = std::max(max_x, atom.position.x() + atom.radius);
min_y = std::min(min_y, atom.position.y() - atom.radius);
max_y = std::max(max_y, atom.position.y() + atom.radius);
min_z = std::min(min_z, atom.position.z() - atom.radius);
max_z = std::max(max_z, atom.position.z() + atom.radius);
}
return CGAL::Bbox_3(min_x, min_y, min_z, max_x, max_y, max_z);
}
private:
std::vector<Atom> atoms_;
};
// 创建示例水分子
Molecule create_water_molecule() {
Molecule water;
// 氧原子
Atom oxygen;
oxygen.element = "O";
oxygen.position = Point_3(0, 0, 0);
oxygen.radius = 1.52; // VDW半径(埃)
oxygen.charge = -0.834;
water.add_atom(oxygen);
// 氢原子1
Atom hydrogen1;
hydrogen1.element = "H";
hydrogen1.position = Point_3(0.96, 0, 0); // O-H键长约0.96埃
hydrogen1.radius = 1.20;
hydrogen1.charge = 0.417;
water.add_atom(hydrogen1);
// 氢原子2
Atom hydrogen2;
hydrogen2.element = "H";
hydrogen2.position = Point_3(-0.24, 0.93, 0); // 104.5度键角
hydrogen2.radius = 1.20;
hydrogen2.charge = 0.417;
water.add_atom(hydrogen2);
return water;
}
// 创建示例蛋白质片段(简化)
Molecule create_protein_fragment() {
Molecule protein;
// 模拟一个小的螺旋片段
std::vector<std::string> elements = {"N", "CA", "C", "O"};
std::vector<double> radii = {1.55, 1.70, 1.70, 1.52};
for (int i = 0; i < 12; ++i) { // 3个残基
for (int j = 0; j < 4; ++j) {
Atom atom;
atom.element = elements[j];
atom.radius = radii[j];
atom.residue_id = i / 4;
// 简化的螺旋坐标
double angle = i * 1.0;
double x = 2.0 * std::cos(angle);
double y = 2.0 * std::sin(angle);
double z = i * 1.5;
atom.position = Point_3(x, y, z);
protein.add_atom(atom);
}
}
return protein;
}
// 分子表面分析
SurfaceProperties analyze_surface(const Surface_mesh& surface) {
SurfaceProperties props;
props.total_area = CGAL::to_double(PMP::area(surface));
if (PMP::is_closed(surface)) {
props.volume = CGAL::to_double(PMP::volume(surface));
}
// 曲率分析(简化)
// 实际应用需要计算平均曲率和高斯曲率
props.avg_curvature = 0.0;
return props;
}
// 生成分子报告
void generate_molecular_report(const Molecule& molecule,
const SurfaceProperties& props) {
std::cout << "\n========== 分子分析报告 ==========" << std::endl;
std::cout << "原子数: " << molecule.get_atoms().size() << std::endl;
// 统计元素组成
std::map<std::string, int> element_count;
for (const auto& atom : molecule.get_atoms()) {
element_count[atom.element]++;
}
std::cout << "元素组成:" << std::endl;
for (const auto& ec : element_count) {
std::cout << " " << ec.first << ": " << ec.second << std::endl;
}
auto centroid = molecule.compute_centroid();
std::cout << "分子质心: (" << centroid.x() << ", "
<< centroid.y() << ", " << centroid.z() << ")" << std::endl;
auto bbox = molecule.get_bounding_box();
std::cout << "包围盒尺寸: "
<< (bbox.xmax() - bbox.xmin()) << " x "
<< (bbox.ymax() - bbox.ymin()) << " x "
<< (bbox.zmax() - bbox.zmin()) << " 埃" << std::endl;
std::cout << "\n表面属性:" << std::endl;
std::cout << " 表面积: " << std::fixed << std::setprecision(2)
<< props.total_area << " 埃²" << std::endl;
std::cout << " 体积: " << props.volume << " 埃³" << std::endl;
}
int main() {
std::cout << "=== 分子表面计算系统 ===" << std::endl;
// 分析水分子
std::cout << "\n--- 水分子(H2O) ---" << std::endl;
Molecule water = create_water_molecule();
Surface_mesh water_surface = water.compute_vdw_surface();
SurfaceProperties water_props = analyze_surface(water_surface);
generate_molecular_report(water, water_props);
// 分析蛋白质片段
std::cout << "\n--- 蛋白质片段 ---" << std::endl;
Molecule protein = create_protein_fragment();
Surface_mesh protein_surface = protein.compute_vdw_surface();
SurfaceProperties protein_props = analyze_surface(protein_surface);
generate_molecular_report(protein, protein_props);
return 0;
}蛋白质口袋检测
// molecular_pocket_detection.cpp
// 蛋白质口袋检测:识别潜在的结合位点
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <iostream>
#include <vector>
#include <queue>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef K::Weighted_point_3 Weighted_point_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
// 口袋定义
struct Pocket {
int id;
double volume;
double surface_area;
Point_3 centroid;
double depth; // 口袋深度
double druggability; // 可药性评分(0-1)
std::vector<Point_3> cavity_points;
};
// 口袋检测参数
struct PocketDetectionParams {
double probe_radius; // 探针球半径(通常1.4埃,水分子大小)
double min_pocket_volume; // 最小口袋体积
double max_pocket_depth; // 最大口袋深度
double grid_spacing; // 网格间距
};
class PocketDetector {
public:
PocketDetector(const std::vector<Point_3>& protein_atoms,
const std::vector<double>& atom_radii)
: protein_atoms_(protein_atoms), atom_radii_(atom_radii) {}
// 基于网格的口袋检测
std::vector<Pocket> detect_pockets(const PocketDetectionParams& params) {
std::vector<Pocket> pockets;
// 1. 构建蛋白质包围盒
CGAL::Bbox_3 bbox = compute_bounding_box();
// 2. 创建3D网格
// 标记每个网格点是否为蛋白质内部/表面/外部
// 3. 识别空腔区域(外部但不可从远处到达)
auto cavity_grid = identify_cavities(bbox, params);
// 4. 聚类空腔点形成口袋
auto clusters = cluster_cavities(cavity_grid, params);
// 5. 计算每个口袋的属性
int pocket_id = 0;
for (const auto& cluster : clusters) {
if (cluster.size() * params.grid_spacing * params.grid_spacing * params.grid_spacing
< params.min_pocket_volume) {
continue;
}
Pocket pocket;
pocket.id = pocket_id++;
pocket.cavity_points = cluster;
pocket.volume = cluster.size() * params.grid_spacing * params.grid_spacing * params.grid_spacing;
// 计算质心
double cx = 0, cy = 0, cz = 0;
for (const auto& p : cluster) {
cx += p.x();
cy += p.y();
cz += p.z();
}
pocket.centroid = Point_3(cx / cluster.size(), cy / cluster.size(), cz / cluster.size());
// 计算深度(到蛋白质表面的距离)
pocket.depth = compute_pocket_depth(pocket.centroid);
// 计算可药性(基于几何特征)
pocket.druggability = compute_druggability(pocket);
pockets.push_back(pocket);
}
return pockets;
}
// 计算口袋开口大小
double compute_pocket_opening(const Pocket& pocket) {
// 找到口袋中距离蛋白质表面最近的点
double min_dist = std::numeric_limits<double>::max();
for (const auto& p : pocket.cavity_points) {
double dist = distance_to_protein_surface(p);
min_dist = std::min(min_dist, dist);
}
return min_dist;
}
private:
std::vector<Point_3> protein_atoms_;
std::vector<double> atom_radii_;
CGAL::Bbox_3 compute_bounding_box() const {
if (protein_atoms_.empty()) return CGAL::Bbox_3();
double min_x = protein_atoms_[0].x();
double max_x = protein_atoms_[0].x();
double min_y = protein_atoms_[0].y();
double max_y = protein_atoms_[0].y();
double min_z = protein_atoms_[0].z();
double max_z = protein_atoms_[0].z();
for (const auto& p : protein_atoms_) {
min_x = std::min(min_x, p.x());
max_x = std::max(max_x, p.x());
min_y = std::min(min_y, p.y());
max_y = std::max(max_y, p.y());
min_z = std::min(min_z, p.z());
max_z = std::max(max_z, p.z());
}
// 添加边界
double padding = 5.0;
return CGAL::Bbox_3(min_x - padding, min_y - padding, min_z - padding,
max_x + padding, max_y + padding, max_z + padding);
}
// 识别空腔区域
std::vector<Point_3> identify_cavities(const CGAL::Bbox_3& bbox,
const PocketDetectionParams& params) {
std::vector<Point_3> cavities;
// 简化的网格遍历
// 实际应用需要更复杂的算法(如滚球算法或Alpha Shape)
for (double x = bbox.xmin(); x <= bbox.xmax(); x += params.grid_spacing) {
for (double y = bbox.ymin(); y <= bbox.ymax(); y += params.grid_spacing) {
for (double z = bbox.zmin(); z <= bbox.zmax(); z += params.grid_spacing) {
Point_3 grid_point(x, y, z);
// 检查点是否在蛋白质外部
if (!is_inside_protein(grid_point)) {
// 检查点是否被蛋白质包围(不可从远处到达)
if (is_buried(grid_point, params.probe_radius)) {
cavities.push_back(grid_point);
}
}
}
}
}
return cavities;
}
// 检查点是否在蛋白质内部
bool is_inside_protein(const Point_3& p) const {
for (size_t i = 0; i < protein_atoms_.size(); ++i) {
double dist_sq = CGAL::to_double((p - protein_atoms_[i]).squared_length());
double radius = atom_radii_[i];
if (dist_sq < radius * radius) {
return true;
}
}
return false;
}
// 检查点是否被埋藏(不可从远处到达)
bool is_buried(const Point_3& p, double probe_radius) const {
// 简化检查:向6个方向发射射线
std::vector<Point_3> directions = {
Point_3(1, 0, 0), Point_3(-1, 0, 0),
Point_3(0, 1, 0), Point_3(0, -1, 0),
Point_3(0, 0, 1), Point_3(0, 0, -1)
};
int blocked_count = 0;
for (const auto& dir : directions) {
// 检查射线是否与蛋白质相交
// 简化:检查一定距离内是否有原子
bool blocked = false;
for (size_t i = 0; i < protein_atoms_.size(); ++i) {
double dist_sq = CGAL::to_double((p - protein_atoms_[i]).squared_length());
double threshold = (atom_radii_[i] + probe_radius) * (atom_radii_[i] + probe_radius);
if (dist_sq < threshold) {
blocked = true;
break;
}
}
if (blocked) blocked_count++;
}
// 如果大部分方向被阻挡,认为是埋藏点
return blocked_count >= 4;
}
// 聚类空腔点
std::vector<std::vector<Point_3>> cluster_cavities(
const std::vector<Point_3>& cavities,
const PocketDetectionParams& params) {
std::vector<std::vector<Point_3>> clusters;
std::vector<bool> visited(cavities.size(), false);
for (size_t i = 0; i < cavities.size(); ++i) {
if (visited[i]) continue;
// BFS聚类
std::vector<Point_3> cluster;
std::queue<size_t> queue;
queue.push(i);
visited[i] = true;
while (!queue.empty()) {
size_t current = queue.front();
queue.pop();
cluster.push_back(cavities[current]);
// 查找邻近点
for (size_t j = 0; j < cavities.size(); ++j) {
if (visited[j]) continue;
double dist_sq = CGAL::to_double(
(cavities[current] - cavities[j]).squared_length()
);
if (dist_sq <= params.grid_spacing * params.grid_spacing * 3) {
visited[j] = true;
queue.push(j);
}
}
}
clusters.push_back(cluster);
}
return clusters;
}
double compute_pocket_depth(const Point_3& centroid) {
// 计算质心到最近蛋白质表面的距离
double min_dist = std::numeric_limits<double>::max();
for (size_t i = 0; i < protein_atoms_.size(); ++i) {
double dist = std::sqrt(CGAL::to_double(
(centroid - protein_atoms_[i]).squared_length()
));
dist -= atom_radii_[i];
min_dist = std::min(min_dist, dist);
}
return min_dist;
}
double compute_druggability(const Pocket& pocket) {
// 基于几何特征的可药性评分
// 考虑因素:体积、深度、形状等
double score = 0.0;
// 体积评分(理想药物口袋体积约500-1000埃³)
if (pocket.volume >= 200 && pocket.volume <= 2000) {
score += 0.3;
}
// 深度评分
if (pocket.depth >= 3.0 && pocket.depth <= 10.0) {
score += 0.3;
}
// 形状评分(球形度)
// 简化处理
score += 0.2;
// 疏水性(需要额外的疏水性计算)
score += 0.2;
return std::min(score, 1.0);
}
double distance_to_protein_surface(const Point_3& p) const {
double min_dist = std::numeric_limits<double>::max();
for (size_t i = 0; i < protein_atoms_.size(); ++i) {
double dist = std::sqrt(CGAL::to_double(
(p - protein_atoms_[i]).squared_length()
));
dist -= atom_radii_[i];
min_dist = std::min(min_dist, dist);
}
return min_dist;
}
};
// 生成口袋检测报告
void generate_pocket_report(const std::vector<Pocket>& pockets) {
std::cout << "\n========== 口袋检测报告 ==========" << std::endl;
std::cout << "发现 " << pockets.size() << " 个潜在结合口袋" << std::endl;
if (pockets.empty()) return;
std::cout << "\n" << std::left << std::setw(8) << "ID"
<< std::setw(12) << "体积(ų)"
<< std::setw(10) << "深度(Å)"
<< std::setw(12) << "可药性"
<< std::setw(20) << "质心坐标" << std::endl;
std::cout << std::string(62, '-') << std::endl;
for (const auto& pocket : pockets) {
std::cout << std::left << std::setw(8) << pocket.id
<< std::setw(12) << std::fixed << std::setprecision(1) << pocket.volume
<< std::setw(10) << std::setprecision(2) << pocket.depth
<< std::setw(12) << std::setprecision(3) << pocket.druggability
<< std::setw(20) << "(" + std::to_string(static_cast<int>(pocket.centroid.x())) + ", "
+ std::to_string(static_cast<int>(pocket.centroid.y())) + ", "
+ std::to_string(static_cast<int>(pocket.centroid.z())) + ")"
<< std::endl;
}
}
int main() {
std::cout << "=== 蛋白质口袋检测系统 ===" << std::endl;
// 创建示例蛋白质(简化的原子坐标)
std::vector<Point_3> protein_atoms;
std::vector<double> atom_radii;
// 添加一些原子形成口袋结构
// 外围原子
for (int i = 0; i < 20; ++i) {
double angle = i * 2 * M_PI / 20;
protein_atoms.push_back(Point_3(
10 * std::cos(angle),
10 * std::sin(angle),
0
));
atom_radii.push_back(1.7); // 碳原子半径
}
// 顶部覆盖
for (int i = 0; i < 10; ++i) {
for (int j = 0; j < 10; ++j) {
double x = -8 + i * 1.6;
double y = -8 + j * 1.6;
if (x*x + y*y < 64) {
protein_atoms.push_back(Point_3(x, y, 8));
atom_radii.push_back(1.7);
}
}
}
std::cout << "蛋白质原子数: " << protein_atoms.size() << std::endl;
// 检测口袋
PocketDetector detector(protein_atoms, atom_radii);
PocketDetectionParams params;
params.probe_radius = 1.4;
params.min_pocket_volume = 50.0;
params.max_pocket_depth = 20.0;
params.grid_spacing = 1.0;
auto pockets = detector.detect_pockets(params);
// 生成报告
generate_pocket_report(pockets);
return 0;
}21.1.4.3 扩展资源
- OpenBabel: 化学文件格式转换和分子操作
- RDKit: 化学信息学和机器学习工具包
- 相关CGAL包: 3D Triangulation, Alpha Shapes, Surface Mesh
21.1.5 案例五:计算机视觉与SLAM
21.1.5.1 领域背景
计算机视觉和同时定位与地图构建(SLAM)是机器人学和增强现实的核心技术。CGAL在这些领域的应用包括:
- 点云处理:滤波、配准、分割
- 表面重建:从点云生成三角网格
- 网格简化:实时渲染的LOD生成
- 碰撞检测:机器人路径规划
21.1.5.2 核心算法
点云配准
// vision_point_cloud_registration.cpp
// 点云配准:ICP算法实现
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Kd_tree.h>
#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <iostream>
#include <vector>
#include <random>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef K::Vector_3 Vector_3;
typedef K::Aff_transformation_3 Aff_transformation_3;
typedef CGAL::Search_traits_3<K> Search_traits;
typedef CGAL::Kd_tree<Search_traits> Kd_tree;
typedef CGAL::Fuzzy_sphere<Search_traits> Fuzzy_sphere;
// 点云定义
typedef std::vector<Point_3> PointCloud;
// 配准结果
struct RegistrationResult {
Aff_transformation_3 transformation;
double rmse; // 均方根误差
int iterations; // 迭代次数
bool converged; // 是否收敛
std::vector<double> errors; // 每轮迭代的误差
};
// ICP参数
struct ICPParameters {
int max_iterations;
double convergence_threshold;
double max_correspondence_distance;
bool use_point_to_plane; // 使用点到平面ICP
};
// 点云配准器
class PointCloudRegistration {
public:
// 最近点查找
Point_3 find_nearest_point(const Kd_tree& tree, const Point_3& query) {
auto result = tree.search_nearest_point(query);
return result.first;
}
// 点到点ICP
RegistrationResult icp_point_to_point(const PointCloud& source,
const PointCloud& target,
const ICPParameters& params) {
RegistrationResult result;
result.converged = false;
// 构建目标点云的KD树
Kd_tree target_tree(target.begin(), target.end());
target_tree.build();
PointCloud transformed_source = source;
Aff_transformation_3 current_transform = Aff_transformation_3(
CGAL::IDENTITY
);
for (int iter = 0; iter < params.max_iterations; ++iter) {
// 1. 寻找对应点
std::vector<std::pair<Point_3, Point_3>> correspondences;
double total_error = 0.0;
for (const auto& p : transformed_source) {
Point_3 nearest = find_nearest_point(target_tree, p);
double dist_sq = CGAL::to_double((p - nearest).squared_length());
if (dist_sq < params.max_correspondence_distance *
params.max_correspondence_distance) {
correspondences.push_back({p, nearest});
total_error += std::sqrt(dist_sq);
}
}
if (correspondences.empty()) break;
double rmse = total_error / correspondences.size();
result.errors.push_back(rmse);
// 检查收敛
if (rmse < params.convergence_threshold) {
result.converged = true;
result.iterations = iter + 1;
result.rmse = rmse;
result.transformation = current_transform;
break;
}
// 2. 计算刚体变换(SVD方法)
Aff_transformation_3 transform = compute_rigid_transform(correspondences);
// 3. 应用变换
for (auto& p : transformed_source) {
p = transform.transform(p);
}
// 累积变换
current_transform = transform * current_transform;
}
if (!result.converged) {
result.iterations = params.max_iterations;
result.rmse = result.errors.empty() ? 0 : result.errors.back();
result.transformation = current_transform;
}
return result;
}
// 计算质心
Point_3 compute_centroid(const PointCloud& points) {
double cx = 0, cy = 0, cz = 0;
for (const auto& p : points) {
cx += p.x();
cy += p.y();
cz += p.z();
}
double n = points.size();
return Point_3(cx/n, cy/n, cz/n);
}
private:
// 使用SVD计算刚体变换
Aff_transformation_3 compute_rigid_transform(
const std::vector<std::pair<Point_3, Point_3>>& correspondences) {
// 计算质心
Point_3 source_centroid = compute_centroid(
[&correspondences]() {
PointCloud pts;
for (const auto& c : correspondences) pts.push_back(c.first);
return pts;
}()
);
Point_3 target_centroid = compute_centroid(
[&correspondences]() {
PointCloud pts;
for (const auto& c : correspondences) pts.push_back(c.second);
return pts;
}()
);
// 中心化点云
// 简化的变换计算(实际应用需要完整的SVD)
// 计算旋转和平移
Vector_3 translation(
target_centroid.x() - source_centroid.x(),
target_centroid.y() - source_centroid.y(),
target_centroid.z() - source_centroid.z()
);
// 返回平移变换(简化,实际应用需要计算旋转)
return Aff_transformation_3(
1, 0, 0, translation.x(),
0, 1, 0, translation.y(),
0, 0, 1, translation.z()
);
}
};
// 生成测试点云(带噪声的立方体)
PointCloud generate_test_point_cloud(int num_points, double noise_level) {
PointCloud cloud;
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis(-1.0, 1.0);
std::normal_distribution<> noise(0.0, noise_level);
for (int i = 0; i < num_points; ++i) {
// 在单位球面上生成点
double theta = 2 * M_PI * dis(gen);
double phi = std::acos(2 * dis(gen) - 1);
double x = std::sin(phi) * std::cos(theta);
double y = std::sin(phi) * std::sin(theta);
double z = std::cos(phi);
// 添加噪声
cloud.push_back(Point_3(
x + noise(gen),
y + noise(gen),
z + noise(gen)
));
}
return cloud;
}
// 应用刚体变换
PointCloud transform_point_cloud(const PointCloud& cloud,
const Aff_transformation_3& transform) {
PointCloud result;
for (const auto& p : cloud) {
result.push_back(transform.transform(p));
}
return result;
}
// 生成配准报告
void generate_registration_report(const RegistrationResult& result) {
std::cout << "\n========== 配准结果报告 ==========" << std::endl;
std::cout << "收敛状态: " << (result.converged ? "已收敛" : "未收敛") << std::endl;
std::cout << "迭代次数: " << result.iterations << std::endl;
std::cout << "最终RMSE: " << std::fixed << std::setprecision(6)
<< result.rmse << std::endl;
if (!result.errors.empty()) {
std::cout << "\n误差收敛曲线:" << std::endl;
for (size_t i = 0; i < result.errors.size(); ++i) {
std::cout << " 迭代 " << i + 1 << ": "
<< result.errors[i] << std::endl;
}
}
}
int main() {
std::cout << "=== 点云配准系统(ICP) ===" << std::endl;
// 生成源点云
PointCloud source = generate_test_point_cloud(100, 0.01);
std::cout << "源点云点数: " << source.size() << std::endl;
// 创建变换(平移+旋转)
Aff_transformation_3 true_transform(
0.866, -0.5, 0, 2.0, // 30度旋转 + 平移
0.5, 0.866, 0, 1.0,
0, 0, 1, 0.5
);
// 生成目标点云
PointCloud target = transform_point_cloud(source, true_transform);
std::cout << "目标点云点数: " << target.size() << std::endl;
// 执行ICP配准
PointCloudRegistration registration;
ICPParameters params;
params.max_iterations = 50;
params.convergence_threshold = 1e-6;
params.max_correspondence_distance = 10.0;
params.use_point_to_plane = false;
std::cout << "\n执行ICP配准..." << std::endl;
RegistrationResult result = registration.icp_point_to_point(
source, target, params
);
// 生成报告
generate_registration_report(result);
return 0;
}表面重建管线
// vision_surface_reconstruction.cpp
// 表面重建:从点云到三角网格
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Advancing_front_surface_reconstruction.h>
#include <CGAL/Poisson_reconstruction_function.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
#include <CGAL/jet_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
#include <iostream>
#include <vector>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef K::Vector_3 Vector_3;
typedef CGAL::Point_with_normal_3<K> Point_with_normal;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
// 点云处理参数
struct PointCloudProcessingParams {
int normal_estimation_neighbors; // 法向量估计邻居数
double poisson_depth; // Poisson重建深度
double smoothing_iterations; // 平滑迭代次数
double hole_filling_factor; // 孔洞填充因子
};
// 表面重建器
class SurfaceReconstructor {
public:
// 估计法向量
std::vector<Point_with_normal> estimate_normals(
const std::vector<Point_3>& points,
int num_neighbors = 18) {
std::vector<Point_with_normal> points_with_normals;
// 使用Jet拟合估计法向量
// 注意:完整实现需要CGAL的jet_estimate_normals函数
// 简化的法向量估计(使用局部PCA)
for (size_t i = 0; i < points.size(); ++i) {
// 找到k近邻
std::vector<Point_3> neighbors = find_k_nearest(points, i, num_neighbors);
// 计算局部法向量(简化)
Vector_3 normal = estimate_normal_from_neighbors(neighbors);
points_with_normals.push_back(Point_with_normal(points[i], normal));
}
// 法向量定向
orient_normals(points_with_normals);
return points_with_normals;
}
// Poisson表面重建
Surface_mesh poisson_reconstruction(
const std::vector<Point_with_normal>& points,
double depth = 8) {
Surface_mesh mesh;
// Poisson重建步骤:
// 1. 构建八叉树
// 2. 定义指示函数
// 3. 求解Poisson方程
// 4. 提取等值面
// 简化的实现说明
std::cout << "Poisson重建参数:" << std::endl;
std::cout << " 深度: " << depth << std::endl;
std::cout << " 输入点数: " << points.size() << std::endl;
// 实际实现需要使用CGAL::Poisson_reconstruction_function
return mesh;
}
// 前进前沿表面重建
Surface_mesh advancing_front_reconstruction(
const std::vector<Point_3>& points,
double radius_ratio_bound = 5.0,
double beta = 0.52) {
Surface_mesh mesh;
// Advancing Front Surface Reconstruction
// 适用于有序或半有序点云
typedef CGAL::Advancing_front_surface_reconstruction_traits_3<K> Traits;
typedef CGAL::Advancing_front_surface_reconstruction<Traits> Reconstruction;
// 构建Delaunay三角化
// 使用优先队列推进前沿
// 生成三角网格
std::cout << "前进前沿重建参数:" << std::endl;
std::cout << " 半径比界限: " << radius_ratio_bound << std::endl;
std::cout << " Beta角度: " << beta << std::endl;
return mesh;
}
// 点云下采样
std::vector<Point_3> voxel_grid_filter(const std::vector<Point_3>& points,
double voxel_size) {
std::vector<Point_3> filtered;
std::map<std::tuple<int, int, int>, std::vector<Point_3>> voxel_map;
// 将点分配到体素网格
for (const auto& p : points) {
int vx = static_cast<int>(std::floor(p.x() / voxel_size));
int vy = static_cast<int>(std::floor(p.y() / voxel_size));
int vz = static_cast<int>(std::floor(p.z() / voxel_size));
voxel_map[{vx, vy, vz}].push_back(p);
}
// 每个体素保留一个代表点(质心)
for (const auto& voxel : voxel_map) {
const auto& voxel_points = voxel.second;
double cx = 0, cy = 0, cz = 0;
for (const auto& p : voxel_points) {
cx += p.x();
cy += p.y();
cz += p.z();
}
filtered.push_back(Point_3(
cx / voxel_points.size(),
cy / voxel_points.size(),
cz / voxel_points.size()
));
}
return filtered;
}
// 统计异常值移除
std::vector<Point_3> statistical_outlier_removal(
const std::vector<Point_3>& points,
int num_neighbors = 50,
double std_dev_threshold = 1.0) {
std::vector<Point_3> filtered;
std::vector<double> avg_distances;
// 计算每个点到邻居的平均距离
for (size_t i = 0; i < points.size(); ++i) {
auto neighbors = find_k_nearest(points, i, num_neighbors);
double avg_dist = 0.0;
for (const auto& n : neighbors) {
avg_dist += std::sqrt(CGAL::to_double(
(points[i] - n).squared_length()
));
}
avg_dist /= neighbors.size();
avg_distances.push_back(avg_dist);
}
// 计算全局均值和标准差
double mean = 0.0;
for (double d : avg_distances) mean += d;
mean /= avg_distances.size();
double variance = 0.0;
for (double d : avg_distances) variance += (d - mean) * (d - mean);
variance /= avg_distances.size();
double std_dev = std::sqrt(variance);
// 移除远离均值的点
for (size_t i = 0; i < points.size(); ++i) {
if (std::abs(avg_distances[i] - mean) < std_dev_threshold * std_dev) {
filtered.push_back(points[i]);
}
}
return filtered;
}
// 网格简化
void simplify_mesh(Surface_mesh& mesh, double target_edge_length) {
// 使用CGAL的网格简化算法
// PMP::isotropic_remeshing() 或 edge_collapse
std::cout << "网格简化目标边长: " << target_edge_length << std::endl;
std::cout << "原始网格面数: " << mesh.number_of_faces() << std::endl;
// 简化的边折叠(实际应用需要完整的简化算法)
// CGAL::Surface_mesh_simplification::edge_collapse()
std::cout << "简化后面数: " << mesh.number_of_faces() << std::endl;
}
private:
std::vector<Point_3> find_k_nearest(const std::vector<Point_3>& points,
size_t index,
int k) {
// 简化的k近邻查找
// 实际应用需要KD树加速
std::vector<std::pair<double, size_t>> distances;
for (size_t i = 0; i < points.size(); ++i) {
if (i != index) {
double dist = CGAL::to_double(
(points[index] - points[i]).squared_length()
);
distances.push_back({dist, i});
}
}
std::sort(distances.begin(), distances.end());
std::vector<Point_3> neighbors;
for (int i = 0; i < std::min(k, static_cast<int>(distances.size())); ++i) {
neighbors.push_back(points[distances[i].second]);
}
return neighbors;
}
Vector_3 estimate_normal_from_neighbors(const std::vector<Point_3>& neighbors) {
// 简化的法向量估计
// 实际应用需要PCA分析
if (neighbors.size() < 3) {
return Vector_3(0, 0, 1);
}
// 使用前三个点计算近似法向量
Vector_3 v1 = neighbors[1] - neighbors[0];
Vector_3 v2 = neighbors[2] - neighbors[0];
Vector_3 normal = CGAL::cross_product(v1, v2);
double len = std::sqrt(CGAL::to_double(normal.squared_length()));
if (len > 1e-10) {
return normal / len;
}
return Vector_3(0, 0, 1);
}
void orient_normals(std::vector<Point_with_normal>& points) {
// 简化的法向量定向
// 实际应用需要MST传播或视点定向
// 假设所有法向量朝外(对于封闭表面)
// 或使用相机位置定向
}
};
// 生成测试点云(球面)
std::vector<Point_3> generate_sphere_point_cloud(int num_points, double radius) {
std::vector<Point_3> points;
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis(0.0, 1.0);
for (int i = 0; i < num_points; ++i) {
// 均匀球面采样
double theta = 2 * M_PI * dis(gen);
double phi = std::acos(2 * dis(gen) - 1);
double x = radius * std::sin(phi) * std::cos(theta);
double y = radius * std::sin(phi) * std::sin(theta);
double z = radius * std::cos(phi);
points.push_back(Point_3(x, y, z));
}
return points;
}
// 生成重建报告
void generate_reconstruction_report(const Surface_mesh& mesh,
const std::vector<Point_3>& original_points) {
std::cout << "\n========== 表面重建报告 ==========" << std::endl;
std::cout << "原始点云点数: " << original_points.size() << std::endl;
std::cout << "重建网格顶点数: " << mesh.number_of_vertices() << std::endl;
std::cout << "重建网格面数: " << mesh.number_of_faces() << std::endl;
if (mesh.number_of_faces() > 0) {
std::cout << "网格是否封闭: " << (PMP::is_closed(mesh) ? "是" : "否") << std::endl;
if (PMP::is_closed(mesh)) {
double volume = CGAL::to_double(PMP::volume(mesh));
std::cout << "网格体积: " << std::fixed << std::setprecision(2)
<< volume << std::endl;
}
double area = CGAL::to_double(PMP::area(mesh));
std::cout << "网格表面积: " << area << std::endl;
}
}
int main() {
std::cout << "=== 表面重建管线 ===" << std::endl;
SurfaceReconstructor reconstructor;
// 生成测试点云
std::cout << "\n生成测试点云(球面)..." << std::endl;
auto points = generate_sphere_point_cloud(1000, 5.0);
std::cout << "原始点数: " << points.size() << std::endl;
// 点云预处理
std::cout << "\n点云预处理..." << std::endl;
// 下采样
auto filtered = reconstructor.voxel_grid_filter(points, 0.5);
std::cout << "体素滤波后点数: " << filtered.size() << std::endl;
// 异常值移除
auto cleaned = reconstructor.statistical_outlier_removal(filtered);
std::cout << "异常值移除后点数: " << cleaned.size() << std::endl;
// 法向量估计
auto points_with_normals = reconstructor.estimate_normals(cleaned);
std::cout << "法向量估计完成" << std::endl;
// 表面重建
std::cout << "\n执行表面重建..." << std::endl;
Surface_mesh mesh = reconstructor.poisson_reconstruction(points_with_normals);
// 生成报告
generate_reconstruction_report(mesh, points);
return 0;
}21.1.5.3 扩展资源
- PCL (Point Cloud Library): 点云处理的标准库
- Open3D: 现代3D数据处理库
- 相关CGAL包: Point Set Processing, Surface Reconstruction, Shape Detection
21.1.6 案例对比总结
| 应用领域 | 核心CGAL包 | 主要算法 | 数据规模 | 精度要求 |
|---|---|---|---|---|
| BIM建筑 | PMP, AABB Tree, Nef | 碰撞检测、布尔运算 | 10^4-10^6面 | 毫米级 |
| 3D打印 | PMP, 2D Polygon | 切片、支撑生成 | 10^4-10^5面 | 0.1mm级 |
| GIS地理 | 2D/3D Triangulation | TIN、插值、可见性 | 10^6-10^8点 | 米级 |
| 分子建模 | 3D Triangulation, Alpha Shapes | 表面计算、口袋检测 | 10^3-10^5原子 | 埃级 |
| 计算机视觉 | Point Set Processing | 配准、重建 | 10^5-10^7点 | 毫米级 |
共同技术要点
- 空间索引结构: AABB树、KD树在所有领域都有广泛应用
- 三角化技术: 2D/3D Delaunay三角化是基础工具
- 网格处理: 表面网格的修复、简化和优化
- 数值稳定性: 精确几何内核在关键计算中的重要性
21.1.7 扩展学习资源
21.1.7.1 官方资源
- CGAL官方文档: https://doc.cgal.org/
- CGAL示例代码: 每个包都包含完整的示例程序
- CGAL手册: 详细的算法说明和复杂度分析
21.1.7.2 领域特定资源
BIM与建筑
- IFC标准: buildingSMART官方文档
- OpenCASCADE: 工业级几何内核,与CGAL互补
- IfcOpenShell: IFC文件解析和几何转换
3D打印
- Cura Engine: 开源切片引擎源码
- PrusaSlicer: 高级切片算法实现
- 3MF标准: 微软主导的3D打印文件格式
GIS
- GDAL/OGR: 地理空间数据读写
- PROJ: 坐标系统转换
- PostGIS: 空间数据库
分子建模
- OpenBabel: 化学格式转换
- RDKit: 化学信息学工具包
- PDB标准: 蛋白质数据库格式
计算机视觉
- PCL: 点云处理库
- Open3D: 现代3D数据处理
- COLMAP: 多视图立体重建
21.1.7.3 学术论文
- CGAL框架: “The CGAL Kernel: A Basis for Geometric Computation”
- 表面重建: “Poisson Surface Reconstruction” (Kazhdan et al.)
- 点云配准: “Efficient Variants of the ICP Algorithm” (Rusinkiewicz)
- Alpha Shapes: “Three-Dimensional Alpha Shapes” (Edelsbrunner)
21.1.8 本章小结
本章通过五个实际应用案例,展示了CGAL在不同领域的应用方法:
-
BIM领域: 利用CGAL的几何内核和碰撞检测算法,实现建筑构件的干涉检查和空间分析。
-
3D打印: 结合网格修复、切片算法和路径规划,构建完整的3D打印前处理流程。
-
GIS系统: 使用三角化和插值算法处理地形数据,实现可视性和水文分析。
-
分子建模: 应用加权Alpha Shape和表面重建算法计算分子表面和识别结合口袋。
-
计算机视觉: 通过点云配准和表面重建,构建从扫描到模型的完整管线。
这些案例的共同特点是:CGAL提供了底层的几何算法实现,开发者需要结合领域知识构建上层应用。掌握CGAL的核心概念(内核、几何对象、算法)是应用开发的基础。
21.1.9 练习
基础练习
-
BIM碰撞检测优化: 修改21.1.1节的碰撞检测代码,使用CGAL的AABB树实现更高效的批量检测。
-
3D打印切片: 实现一个完整的切片算法,将立方体模型转换为G-code。
-
地形插值: 使用CGAL的自然邻域插值,实现从离散高程点到连续DEM的转换。
进阶练习
-
分子表面可视化: 结合CGAL和OpenGL,实现分子表面的实时渲染。
-
点云配准改进: 在ICP算法中加入点到平面的误差度量,提高配准精度。
-
多算法对比: 对比Poisson重建和Advancing Front重建在不同点云上的效果。
综合项目
-
建筑扫描处理: 设计一个从激光扫描数据到BIM模型的完整处理流程。
-
蛋白质分析工具: 开发一个命令行工具,计算PDB文件的分子表面属性和口袋信息。
-
无人机地形测绘: 结合可见性分析和地形插值,实现无人机航迹规划系统。
附录:代码编译说明
本章所有代码示例均使用CGAL 5.x版本。编译命令示例:
# 基础编译
g++ -std=c++17 -I/path/to/cgal/include -I/path/to/boost/include \
-lmpfr -lgmp -o program program.cpp
# 使用CMake(推荐)
# CMakeLists.txt示例:
cmake_minimum_required(VERSION 3.12)
find_package(CGAL REQUIRED)
add_executable(program program.cpp)
target_link_libraries(program CGAL::CGAL)依赖库安装
# Ubuntu/Debian
sudo apt-get install libcgal-dev libgmp-dev libmpfr-dev
# macOS
brew install cgal gmp mpfr
# Windows (vcpkg)
vcpkg install cgal gmp mpfr运行环境要求
- C++17兼容编译器(GCC 7+, Clang 5+, MSVC 2017+)
- CGAL 5.0或更高版本
- GMP和MPFR库
- Boost库(部分示例需要)