diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/acvd_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/acvd_example.cpp index c6a65b3b414..31c0aa848e5 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/acvd_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/acvd_example.cpp @@ -37,10 +37,18 @@ int main(int argc, char* argv[]) std::cout << "Uniform Isotropic ACVD ...." << std::endl; auto acvd_mesh = PMP::acvd_isotropic_remeshing(smesh, nb_clusters); - CGAL::IO::write_OFF("fandisk_qem-pp3000.off", acvd_mesh); + CGAL::IO::write_OFF("fandisk_acvd.off", acvd_mesh); std::cout << "Completed" << std::endl; + // With Post-Processing QEM Optimization + + std::cout << "Uniform Isotropic ACVD with QEM optimization ...." << std::endl; + + auto acvd_mesh_qem = PMP::acvd_isotropic_remeshing(smesh, nb_clusters, CGAL::parameters::post_processing_qem(true)); + CGAL::IO::write_OFF("fandisk_acvd_qem-pp.off", acvd_mesh_qem); + + std::cout << "Completed" << std::endl; /// Adaptive Isotropic ACVD diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/acvd/acvd.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/acvd/acvd.h index 6bbd2273009..d7642589ca8 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/acvd/acvd.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/acvd/acvd.h @@ -84,7 +84,7 @@ void compute_qem_vertex(std::vector> cluster_ } template -typename GT::Vector_3 compute_displacement(const Eigen::Matrix quadric, const typename GT::Vector_3& p, int& RankDeficiency) +typename GT::Vector_3 compute_displacement(const Eigen::Matrix quadric, const typename GT::Point_3& p, int& RankDeficiency) { typedef Eigen::Matrix Matrix4d; typedef Eigen::Matrix Matrix3d; @@ -169,7 +169,7 @@ typename GT::Vector_3 compute_displacement(const Eigen::Matrix -void compute_representative_point(const Eigen::Matrix& quadric, typename GT::Vector_3& p, int& RankDeficiency) +void compute_representative_point(const Eigen::Matrix& quadric, typename GT::Point_3& p, int& RankDeficiency) { // average point typename GT::Vector_3 displacement = compute_displacement(quadric, p, RankDeficiency); @@ -323,6 +323,9 @@ std::pair< choose_parameter(get_parameter(np, internal_np::vertex_principal_curvatures_and_directions_map), Default_principal_map()); + // get parameter for turning on post-processing qem point optimization + const bool post_processing_qem = choose_parameter(get_parameter(np, internal_np::post_processing_qem), false); + // if adaptive clustering if (gradation_factor > 0 && is_default_parameter::value) @@ -659,37 +662,6 @@ std::pair< std::vector> polygons; TriangleMesh simplified_mesh; - // create a point for each cluster - std::vector> cluster_quadrics(clusters.size()); - - // initialize quadrics - for (int i = 0; i < nb_clusters; i++) - cluster_quadrics[i].setZero(); - - for (Face_descriptor fd : faces(pmesh)) { - // get Vs for fd - // compute qem from Vs->"vector_3"s - // add to the 3 indices of the cluster - typename GT::Point_3 p_i = get(vpm, target(halfedge(fd, pmesh), pmesh)); - typename GT::Vector_3 vec_1 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); - p_i = get(vpm, target(next(halfedge(fd, pmesh), pmesh), pmesh)); - typename GT::Vector_3 vec_2 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); - p_i = get(vpm, target(next(next(halfedge(fd, pmesh), pmesh), pmesh), pmesh)); - typename GT::Vector_3 vec_3 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); - - int c_1 = get(vertex_cluster_pmap, target(halfedge(fd, pmesh), pmesh)); - int c_2 = get(vertex_cluster_pmap, target(next(halfedge(fd, pmesh), pmesh), pmesh)); - int c_3 = get(vertex_cluster_pmap, target(next(next(halfedge(fd, pmesh), pmesh), pmesh), pmesh)); - - if (c_1 != -1 && c_2 != -1 && c_3 != -1) - { - Eigen::Matrix q; - compute_qem_face(vec_1, vec_2, vec_3, q); - cluster_quadrics[c_1] += q; - cluster_quadrics[c_2] += q; - cluster_quadrics[c_3] += q; - } - } for (int i = 0; i < nb_clusters; i++) { @@ -698,16 +670,70 @@ std::pair< valid_cluster_map[i] = points.size(); typename GT::Vector_3 cluster_representative = clusters[i].compute_centroid(); - int RankDeficiency = 0; - //compute_representative_point(cluster_quadrics[i], cluster_representative, RankDeficiency); - - typename GT::Point_3 cluster_representative_p(cluster_representative.x(), cluster_representative.y(), cluster_representative.z()); points.push_back(cluster_representative_p); } } + if (post_processing_qem){ + // create a point for each cluster + std::vector> cluster_quadrics(clusters.size()); + + // initialize quadrics + for (int i = 0; i < nb_clusters; i++) + cluster_quadrics[i].setZero(); + + // for storing the vertex_descriptor of each face + std::vector face_vertices; + + for (Face_descriptor fd : faces(pmesh)) { + // get Vs for fd + // compute qem from Vs->"vector_3"s + // add to the 3 indices of the cluster + Halfedge_descriptor hd = halfedge(fd, pmesh); + do { + Vertex_descriptor vd = target(hd, pmesh); + face_vertices.push_back(vd); + hd = next(hd, pmesh); + } while (hd != halfedge(fd, pmesh)); + + auto p_i = get(vpm, face_vertices[0]); + typename GT::Vector_3 vec_1 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); + p_i = get(vpm, face_vertices[1]); + typename GT::Vector_3 vec_2 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); + p_i = get(vpm, face_vertices[2]); + typename GT::Vector_3 vec_3 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); + + int c_1 = get(vertex_cluster_pmap, face_vertices[0]); + int c_2 = get(vertex_cluster_pmap, face_vertices[1]); + int c_3 = get(vertex_cluster_pmap, face_vertices[2]); + + if (c_1 != -1 && c_2 != -1 && c_3 != -1) + { + Eigen::Matrix q; + compute_qem_face(vec_1, vec_2, vec_3, q); + cluster_quadrics[c_1] += q; + cluster_quadrics[c_2] += q; + cluster_quadrics[c_3] += q; + } + + face_vertices.clear(); + } + + int valid_index = 0; + + for (int i = 0; i < nb_clusters; i++) + { + if (clusters[i].weight_sum > 0) + { + int RankDeficiency = 0; + compute_representative_point(cluster_quadrics[i], points[valid_index], RankDeficiency); + valid_index++; + } + } + } + // extract boundary cycles std::vector border_hedges; extract_boundary_cycles(pmesh, std::back_inserter(border_hedges)); diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index 220d08b79b0..aad511fa1ed 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -102,6 +102,7 @@ CGAL_add_named_parameter(vertex_Gaussian_curvature_t, vertex_Gaussian_curvature, CGAL_add_named_parameter(vertex_principal_curvatures_and_directions_t, vertex_principal_curvatures_and_directions, vertex_principal_curvatures_and_directions) CGAL_add_named_parameter(ball_radius_t, ball_radius, ball_radius) CGAL_add_named_parameter(gradation_factor_t, gradation_factor, gradation_factor) +CGAL_add_named_parameter(post_processing_qem_t, post_processing_qem, post_processing_qem) CGAL_add_named_parameter(outward_orientation_t, outward_orientation, outward_orientation) CGAL_add_named_parameter(overlap_test_t, overlap_test, do_overlap_test_of_bounded_sides) CGAL_add_named_parameter(preserve_genus_t, preserve_genus, preserve_genus)