CGAL边折叠edge_collapse的问题

使用edge_collapse对一个模型简化，之后回收垃圾，collect_garbage

处理之前的顶点和三角形数量：

number_of_vertices: 955730 number_of_faces: 1903410
num_vertices: 955730 num_faces: 1903410

处理之后的顶点和三角形数量：
number_of_vertices: 479309 number_of_faces: 950633
num_vertices: 479309 num_faces: 950633
可以看到顶点数和三角形数都减少了

然后遍历每个三角形，计算每个点的法线：

for (auto face : mesh.faces()) {

       auto h = mesh.halfedge(face);
       auto v0 = mesh.target(h);
       auto v1 = mesh.target(mesh.next(h));
       auto v2 = mesh.target(mesh.prev(h));

}

奇怪的事情发生了，当读取到ID为5623的面时，其第一个顶点v0 值为v953550，超出了number_of_vertices的大小，引起程序崩溃

测试程序如下：


#include <CGAL/Simple_cartesian.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/IO/polygon_mesh_io.h>
#include <CGAL/IO/STL.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_2.h>#include <CGAL/Polygon_mesh_processing/compute_normal.h>#include <CGAL/Aff_transformation_3.h>
#include <CGAL/Aff_transformation_2.h>#include <CGAL/IO/polygon_soup_io.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>#include <CGAL/Filtered_kernel/internal/Static_filters/Is_degenerate_3.h>#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/stitch_borders.h>
#include <CGAL/Polygon_mesh_processing/repair_degeneracies.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_stop_predicate.h>
#include <CGAL/Surface_mesh_simplification/Edge_collapse_visitor_base.h>#include <CGAL/Simple_cartesian.h>  // 或者你使用的其他内核
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>typedef std::function<void(std::size_t, std::size_t)> _callback;namespace SMS = CGAL::Surface_mesh_simplification;typedef std::function<void(std::size_t, std::size_t)> _callback;typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel;
typedef K::Point_2 Point2;
typedef K::Point_3 Point3;
typedef K::Triangle_3 Triangle3;
typedef K::Vector_3 Vector3;
typedef K::Vector_2 Vector2;
typedef K::Plane_3 Plane3;
typedef CGAL::Bbox_3 Box3;
typedef CGAL::Bbox_2 Box2;
typedef K::Segment_2 Segment2;
typedef K::Segment_3 Segment3;
typedef CGAL::Polygon_2<K> Polygon2;
typedef CGAL::Surface_mesh<Point3> Mesh;
typedef std::shared_ptr<Mesh>	Mesh_ptr;typedef boost::graph_traits<Mesh> GT;
typedef typename GT::face_descriptor face_descriptor;
typedef typename GT::edge_descriptor edge_descriptor;
typedef typename GT::halfedge_descriptor halfedge_descriptor;
typedef typename GT::vertex_descriptor vertex_descriptor;
typedef typename GT::face_iterator face_iterator;
typedef typename GT::edge_iterator edge_iterator;
typedef typename GT::vertex_iterator vertex_iterator;void compute_vertex_normals_manual(const Mesh& mesh,std::vector<Vector3>& vertex_normals) {// 初始化顶点法向量为0int size = mesh.num_vertices();vertex_normals.resize(mesh.num_vertices(), Vector3(0, 0, 0));std::cout << "vertex_normals size " << size << std::endl;// 用于记录每个顶点相邻的面数std::vector<int> vertex_face_count(mesh.num_vertices(), 0);// 遍历所有面片std::ofstream ofs("compute_vertex_normals_manual.txt");ofs << "num_faces " << (int)mesh.num_faces() << std::endl;ofs << "num_vertices " << (int)mesh.num_vertices() << std::endl;for (auto face : mesh.faces()) {ofs << face.id() << std::endl;if (5623 == face.id()) ofs << "a" << std::endl;if (5623 == face.id()) {bool del = mesh.is_removed(face);ofs << "is_removed " << del << std::endl;}// 获取面片的三个顶点auto h = mesh.halfedge(face);if (5623 == face.id()) ofs << "halfedge " << h << std::endl;if (5623 == face.id()) ofs << "halfedge is_removed " << mesh.is_removed(h)<< std::endl;auto v0 = mesh.target(h);auto v1 = mesh.target(mesh.next(h));auto v2 = mesh.target(mesh.prev(h));if (5623 == face.id()) ofs << "mesh.point" << std::endl;if (5623 == face.id()) ofs << v0 << " " << v1 << " " << v2 << std::endl;if (5623 == face.id()) {/*bool del1 = mesh.is_removed(v0);bool del2 = mesh.is_removed(v1);bool del3 = mesh.is_removed(v2);ofs << "is_removed del1 " << del1 << std::endl;ofs << "is_removed del2 " << del2 << std::endl;ofs << "is_removed del3 " << del3 << std::endl;
*/}// 计算面片法向量Point3 p0 = mesh.point(v0);Point3 p1 = mesh.point(v1);Point3 p2 = mesh.point(v2);if (5623 == face.id()) ofs << "p0" << p0.x() << " " << p0.y() << " " << p0.z() << std::endl;if (5623 == face.id()) ofs << "p1" << p1.x() << " " << p1.y() << " " << p1.z() << std::endl;if (5623 == face.id()) ofs << "p2" << p2.x() << " " << p2.y() << " " << p2.z() << std::endl;if (5623 == face.id()) ofs << "c" << std::endl;Vector3 edge1 = p1 - p0;Vector3 edge2 = p2 - p0;Vector3 face_normal = CGAL::cross_product(edge1, edge2);if (5623 == face.id()) ofs << "d" << std::endl;// 归一化面片法向量（可选）double len = CGAL::sqrt(face_normal.squared_length());if (len > 0) face_normal = face_normal / len;if (5623 == face.id()) ofs << "e" << std::endl;// 累加到顶点法向量vertex_normals[v0.idx()] = vertex_normals[v0.idx()] + face_normal;vertex_normals[v1.idx()] = vertex_normals[v1.idx()] + face_normal;vertex_normals[v2.idx()] = vertex_normals[v2.idx()] + face_normal;if (5623 == face.id()) ofs << "f" << std::endl;vertex_face_count[v0.idx()]++;vertex_face_count[v1.idx()]++;vertex_face_count[v2.idx()]++;if (5623 == face.id()) ofs << "g" << std::endl;}// 平均法向量并归一化for (size_t i = 0; i < vertex_normals.size(); ++i) {if (vertex_face_count[i] > 0) {vertex_normals[i] = vertex_normals[i] / vertex_face_count[i];double len = CGAL::sqrt(vertex_normals[i].squared_length());if (len > 0) vertex_normals[i] = vertex_normals[i] / len;}}
}int main()
{typedef CGAL::Point_3<CGAL::Simple_cartesian<float>> STL_point;typedef CGAL::array<int, 3> STL_tringle;std::cout << "begin" << std::endl;std::string file_path = "F:\\WORK\\STL\\Print-test\\test-bug\\1\\Groot_Planter_-_v2Groot_Plante.stl";std::vector<STL_point> points;std::vector<STL_tringle> faces;Mesh_ptr mesh = std::make_shared<Mesh>();//box_top = Box3();try {CGAL::IO::read_polygon_soup(file_path, points, faces);/*STL_FILE::read_STL(_file_path, points, faces, [&](std::size_t a, std::size_t b) {aggregator.updateSectionProgress(a, b);});
*/CGAL::Polygon_mesh_processing::repair_polygon_soup(points, faces);std::cout << "repair_polygon_soup" << std::endl;std::vector <Mesh::vertex_index>vHandles;vHandles.reserve(points.size());mesh->reserve(points.size(), points.size() + faces.size(), faces.size());for (auto& pt : points) {auto vh = mesh->add_vertex(Point3((double)pt[0], (double)pt[1], (double)pt[2]));vHandles.push_back(vh);//box_top += Box3(pt[0], pt[1], pt[2], pt[0], pt[1], pt[2]);}int i = 0;for (auto& tri : faces) {int v0 = tri[0];int v1 = tri[1];int v2 = tri[2];if (v0 == v1 || v1 == v2 || v0 == v2)continue;auto vh0 = vHandles[v0];auto vh1 = vHandles[v1];auto vh2 = vHandles[v2];auto fh = mesh->add_face(vh0, vh1, vh2);if (!fh.is_valid()) {Point3 p0 = mesh->point(vh0);Point3 p1 = mesh->point(vh1);Point3 p2 = mesh->point(vh2);auto vh3 = mesh->add_vertex(p0);auto vh4 = mesh->add_vertex(p1);auto vh5 = mesh->add_vertex(p2);mesh->add_face(vh3, vh4, vh5);}++i;}}catch (const std::exception& e) {std::cerr << "load_STL" << e.what() << std::endl;}Mesh& _mesh = *mesh;size_t total_edges = _mesh.number_of_edges();CGAL::Surface_mesh_simplification::Edge_count_stop_predicate<Mesh> stop(total_edges / 2 - 1); // 目标边数std::cout << "************start***************" << std::endl;//num_vertices 955730 number_of_vertices 479307std::cout << "number_of_vertices: " << _mesh.number_of_vertices() << "  number_of_faces: " << _mesh.number_of_faces() << std::endl;std::cout << "num_vertices: " << _mesh.num_vertices() << "  num_faces: " << _mesh.num_faces() << std::endl;SMS::edge_collapse(_mesh, stop);if (!CGAL::is_valid(_mesh)) {//CGAL::Polygon_mesh_processing::remove_degenerate_faces(_mesh);//CGAL::Polygon_mesh_processing::stitch_borders(_mesh);}std::cout << "************edge_collapse***************"   << std::endl;//num_vertices 955730 number_of_vertices 479307std::cout << "number_of_vertices: " << _mesh.number_of_vertices() << "  number_of_faces: " << _mesh.number_of_faces() << std::endl;std::cout << "num_vertices: " << _mesh.num_vertices() << "  num_faces: " << _mesh.num_faces() << std::endl;_mesh.collect_garbage();//_mesh.collect_garbage();std::cout << "************collect_garbage***************" << std::endl;//num_vertices 955730 number_of_vertices 479307std::cout << "number_of_vertices: " << _mesh.number_of_vertices() << "  number_of_faces: " << _mesh.number_of_faces() << std::endl;std::cout << "num_vertices: " << _mesh.num_vertices() << "  num_faces: " << _mesh.num_faces() << std::endl;////auto vertex_normals = _mesh.add_property_map<vertex_descriptor, Vector3>("v:normal").first;//CGAL::Polygon_mesh_processing::compute_vertex_normals(_mesh, vertex_normals);//std::vector<Vector3> vertex_normals;compute_vertex_normals_manual(_mesh, vertex_normals);std::cout << "compute_vertex_normals" << std::endl;
}

日志：