给定地面实况分割,怎么使用标准软件CGAL(计算几何算法库)2提取线粒体中心线
时间: 2024-06-08 13:07:50 浏览: 15
首先,需要将地面实况分割转化为一个三维点云。然后,使用CGAL库中的Alpha Shapes算法提取点云的凸壳。接着,使用CGAL库中的Straight Skeleton算法计算凸壳的直线骨架,并将其作为线粒体中心线输出。
下面是一个简单的示例代码:
```c++
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Straight_skeleton_builder_3.h>
#include <CGAL/IO/Ply_reader.h>
#include <CGAL/IO/Ply_writer.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Iso_cuboid_3 Iso_cuboid_3;
typedef CGAL::Alpha_shape_vertex_base_3<Kernel> Vb;
typedef CGAL::Alpha_shape_cell_base_3<Kernel> Fb;
typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
typedef CGAL::Alpha_shape_3<Tds> Alpha_shape_3;
typedef Alpha_shape_3::Alpha_iterator Alpha_iterator;
typedef Alpha_shape_3::Cell_handle Cell_handle;
typedef Alpha_shape_3::Vertex_handle Vertex_handle;
typedef Alpha_shape_3::Facet Facet;
typedef Alpha_shape_3::Cell Cell;
typedef Alpha_shape_3::NT NT;
typedef CGAL::Surface_mesh<Point_3> Mesh;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef Polyhedron::HalfedgeDS HalfedgeDS;
typedef CGAL::Straight_skeleton_builder_traits_3<Kernel> SsTraits;
typedef CGAL::Straight_skeleton_builder_3<SsTraits> SsBuilder;
int main(int argc, char** argv)
{
// 读入点云
std::vector<Point_3> points;
std::ifstream in(argv[1]);
CGAL::read_xyz_points(in, std::back_inserter(points));
in.close();
// 计算点云的凸壳
Alpha_shape_3 alpha(points.begin(), points.end(), 0.0, Alpha_shape_3::GENERAL);
std::vector<std::pair<NT, Cell_handle>> alpha_cells;
for(Alpha_iterator it = alpha.alpha_begin(); it != alpha.alpha_end(); ++it)
{
alpha_cells.clear();
alpha.get_alpha_shape_cells(std::back_inserter(alpha_cells), *it);
if(alpha_cells.size() > 0)
{
Mesh mesh;
std::map<Cell_handle, int> cell_to_index;
std::vector<std::vector<int>> cell_faces(alpha_cells.size());
for(int i = 0; i < alpha_cells.size(); i++)
{
Cell_handle ch = alpha_cells[i].second;
int index = mesh.add_vertex(ch->voronoi());
cell_to_index[ch] = index;
for(int j = 0; j < 4; j++)
{
Cell_handle neighbor = ch->neighbor(j);
if(!alpha.is_infinite(neighbor) && alpha.classify(ch, j) != Alpha_shape_3::EXTERIOR)
{
if(cell_to_index.find(neighbor) == cell_to_index.end())
{
int index2 = mesh.add_vertex(neighbor->voronoi());
cell_to_index[neighbor] = index2;
}
int v1 = cell_to_index[ch];
int v2 = cell_to_index[neighbor];
cell_faces[i].push_back(mesh.add_face(v1, v2, -1));
}
}
}
// 三角剖分面
std::vector<std::vector<int>> face_indices;
CGAL::Polygon_mesh_processing::triangulate_faces(mesh, std::back_inserter(face_indices));
// 计算法向量
std::vector<Point_3> face_normals(face_indices.size());
for(int i = 0; i < face_indices.size(); i++)
{
std::vector<Point_3> face_points;
for(int j = 0; j < face_indices[i].size(); j++)
{
int index = face_indices[i][j];
face_points.push_back(mesh.point(mesh.vertex_handle(index)));
}
CGAL::Polygon_mesh_processing::compute_face_normal(face_points[0], face_points[1], face_points[2], face_normals[i]);
}
// 翻转法向量
for(int i = 0; i < face_indices.size(); i++)
{
if(CGAL::Polygon_mesh_processing::orientation(face_indices[i], mesh) == CGAL::Polygon_mesh_processing::CLOCKWISE)
{
std::reverse(face_indices[i].begin(), face_indices[i].end());
face_normals[i] = -face_normals[i];
}
}
// 输出到PLY文件
std::ofstream out(argv[2]);
CGAL::write_PLY(out, mesh, CGAL::make_normal_of_point_with_normal_pmap(face_normals));
out.close();
break;
}
}
// 计算直线骨架
Polyhedron polyhedron;
std::ifstream in2(argv[2]);
CGAL::read_PLY(in2, polyhedron);
in2.close();
SsBuilder builder(polyhedron.facets_begin(), polyhedron.facets_end(), polyhedron.points_begin(), polyhedron.points_end());
Polyhedron skeleton_polyhedron;
builder(skeleton_polyhedron);
// 输出到PLY文件
std::ofstream out2(argv[3]);
CGAL::write_PLY(out2, skeleton_polyhedron);
out2.close();
return 0;
}
```
这个示例代码包含了以下步骤:
1. 读入点云数据。
2. 使用Alpha Shapes算法计算点云的凸壳。
3. 将凸壳转化为一个三角网格,并进行三角剖分。
4. 计算每个三角形面的法向量,并翻转面的法向量,使其指向凸壳内部。
5. 将三角网格输出到PLY文件中。
6. 使用Straight Skeleton算法计算凸壳的直线骨架。
7. 将直线骨架输出到PLY文件中。
注意,这个示例代码只适用于简单的点云数据,对于更复杂的点云数据,可能需要进行一些调整。同时,这个示例代码仅供参考,实际使用时需要根据具体情况进行修改。
相关推荐
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)