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) {
	// 初始化顶点法向量为0
	int 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 479307
	std::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 479307
	std::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 479307
	std::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;
}

日志：