diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h index 7adc37f0b4a..952a7e749ec 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h @@ -118,7 +118,7 @@ typename GT::FT interpolated_corrected_area_measure_face(const std::vector interpolated_corrected_anisotropic_measure_fa muXY[ix * 3 + iy] = 0.5 * um * (cross_product(u02[iy] * X, x01) - cross_product(u01[iy] * X, x02)); } } - // Quad: use bilinear interpolation formula + // Quad: use the bilinear interpolation formula else if (n == 4) { // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. @@ -412,6 +412,7 @@ typename GT::FT face_in_ball_ratio(const std::vector& x, std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); xm /= static_cast(n); + // computing squared distance of furthest and closest point to ball center typename GT::FT d_min = (xm - c).squared_length(); typename GT::FT d_max = d_min; @@ -422,12 +423,14 @@ typename GT::FT face_in_ball_ratio(const std::vector& x, d_min = (std::min)(d_sq, d_min); } + // if the furthest point is inside ball, return 1 if (d_max <= r * r) return 1.0; + // if the closest point is outside ball, return 0 else if (r * r <= d_min) return 0.0; + // else, approximate inclusion ratio of the triangle: d_max = sqrt(d_max); d_min = sqrt(d_min); - return (r - d_min) / (d_max - d_min); } @@ -438,17 +441,22 @@ Principal_curvatures_and_directions principal_curvatures_and_directions_from const typename GT::Vector_3 u_GT ) { + // putting anisotropic measure in matrix form Eigen::Matrix v_muXY = Eigen::Matrix::Zero(); for (std::size_t ix = 0; ix < 3; ix++) for (std::size_t iy = 0; iy < 3; iy++) v_muXY(ix, iy) = anisotropic_measure[ix * 3 + iy]; + // constant factor K to force the principal direction eigenvectors to be tangential to the surface Eigen::Matrix u(u_GT.x(), u_GT.y(), u_GT.z()); const typename GT::FT K = 1000 * v_mu0; + // symmetrizing and adding the constant term v_muXY = 0.5 * (v_muXY + v_muXY.transpose()) + K * u * u.transpose(); + + // computing eigenvalues and eigenvectors Eigen::SelfAdjointEigenSolver> eigensolver; eigensolver.computeDirect(v_muXY); @@ -462,6 +470,7 @@ Principal_curvatures_and_directions principal_curvatures_and_directions_from const typename GT::Vector_3 min_eig_vec(eig_vecs(0, 1), eig_vecs(1, 1), eig_vecs(2, 1)); const typename GT::Vector_3 max_eig_vec(eig_vecs(0, 0), eig_vecs(1, 0), eig_vecs(2, 0)); + // returning principal curvatures and directions (with the correct sign) return Principal_curvatures_and_directions( (v_mu0 != 0.0) ? -eig_vals[1] / v_mu0 : 0.0, (v_mu0 != 0.0) ? -eig_vals[0] / v_mu0 : 0.0, @@ -470,6 +479,7 @@ Principal_curvatures_and_directions principal_curvatures_and_directions_from ); } +// measures are computed for faces only if they are adjacent to the vertex template Vertex_measures interpolated_corrected_measures_one_vertex_no_radius( const PolygonMesh& pmesh, @@ -495,9 +505,11 @@ Vertex_measures interpolated_corrected_measures_one_vertex_no_radius( std::vector x; std::vector u; + // compute for each face around the vertex (except the null (boundary) face) for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { if (f != boost::graph_traits::null_face()) { + // looping over vertices in face to get point coordinates and normal vectors for (Vertex_descriptor vi : vertices_around_face(halfedge(f, pmesh), pmesh)) { Point_3 pi = get(vpm, vi); @@ -506,6 +518,7 @@ Vertex_measures interpolated_corrected_measures_one_vertex_no_radius( u.push_back(ui); } + // compute measures for selected curvatures (area is always computed) vertex_measures.area_measure += interpolated_corrected_area_measure_face(u, x); if (is_mean_curvature_selected) @@ -528,7 +541,7 @@ Vertex_measures interpolated_corrected_measures_one_vertex_no_radius( return vertex_measures; } - +// measures are computed for faces only if they are in the ball of radius 'radius' centered at the vertex template Vertex_measures interpolated_corrected_measures_one_vertex( const PolygonMesh& pmesh, @@ -547,6 +560,7 @@ Vertex_measures interpolated_corrected_measures_one_vertex( typedef typename GT::Vector_3 Vector_3; typedef typename GT::FT FT; + // the ball expansion is done using a BFS traversal from the vertex std::queue bfs_queue; std::unordered_set bfs_visited; @@ -568,8 +582,7 @@ Vertex_measures interpolated_corrected_measures_one_vertex( Face_descriptor fi = bfs_queue.front(); bfs_queue.pop(); - // looping over vertices in face to get point coordinates - + // looping over vertices in face to get point coordinates and normal vectors for (Vertex_descriptor vi : vertices_around_face(halfedge(fi, pmesh), pmesh)) { Point_3 pi = get(vpm, vi); @@ -578,8 +591,11 @@ Vertex_measures interpolated_corrected_measures_one_vertex( u.push_back(ui); } + // approximate inclusion ratio of the face in the ball const FT f_ratio = face_in_ball_ratio(x, radius, c); + // if it is not 0 (not completely outside), compute measures for selected curvatures (area is always computed) + // and add neighboring faces to the bfs queue if (f_ratio != 0.0) { vertex_measures.area_measure += f_ratio * interpolated_corrected_area_measure_face(u, x); @@ -613,6 +629,7 @@ Vertex_measures interpolated_corrected_measures_one_vertex( return vertex_measures; } +// computes selected curvatures for one specific vertex template void interpolated_corrected_curvatures_one_vertex( @@ -640,23 +657,29 @@ template::value) compute_vertex_normals(pmesh, vnm, np); + // if no radius is given, we pass -1 which will make the expansion be only on the incident faces instead of a ball typename GT::FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1); + // if the radius is 0, we use a small epsilon to expand the ball (scaled with the average edge length) if (radius == 0) radius = average_edge_length(pmesh) * EXPANDING_RADIUS_EPSILON; + // get parameters (pointers) for curvatures typename GT::FT* vertex_mean_curvature = choose_parameter(get_parameter(np, internal_np::vertex_mean_curvature), nullptr); typename GT::FT* vertex_gaussian_curvature = choose_parameter(get_parameter(np, internal_np::vertex_gaussian_curvature), nullptr); Principal_curvatures_and_directions* vertex_principal_curvatures_and_directions = choose_parameter(get_parameter(np, internal_np::vertex_principal_curvatures_and_directions), nullptr); + // determine which curvatures are selected (by checking if the pointers are not null) const bool is_mean_curvature_selected = (vertex_mean_curvature != nullptr); const bool is_gaussian_curvature_selected = (vertex_gaussian_curvature != nullptr); const bool is_principal_curvatures_and_directions_selected = (vertex_principal_curvatures_and_directions != nullptr); Vertex_measures vertex_measures; + // if the radius is negative, we do not expand the ball (only the incident faces) if (radius < 0) { vertex_measures = interpolated_corrected_measures_one_vertex_no_radius( @@ -683,6 +706,7 @@ template principal_curvatures_and_directions = principal_curvatures_and_directions_from_anisotropic_measures( vertex_measures.anisotropic_measure, @@ -783,11 +808,14 @@ private: vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map), get(Vector_map_tag(), pmesh)); + // if no normal map is given, compute normals if (is_default_parameter::value) compute_vertex_normals(pmesh, vnm, np); + // if no radius is given, we pass -1 which will make the expansion be only on the incident faces instead of a ball const FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1); + // check which curvature maps are provided by the user (determines which curvatures are computed) is_mean_curvature_selected = !is_default_parameter::value; is_gaussian_curvature_selected = !is_default_parameter::value; is_principal_curvatures_and_directions_selected = !is_default_parameter::value; @@ -800,6 +828,7 @@ private: } void set_ball_radius(const FT radius) { + // if given radius is 0, we use a small epsilon to expand the ball (scaled by the average edge length) if (radius == 0) ball_radius = average_edge_length(pmesh) * EXPANDING_RADIUS_EPSILON; else @@ -825,6 +854,8 @@ public: private: + // Computes the (selected) interpolated corrected measures for all faces + // and stores them in the property maps void interpolated_corrected_selected_measures_all_faces() { std::vector x; @@ -854,14 +885,17 @@ private: } } + // expand the measures of the faces incident to v Vertex_measures expand_interpolated_corrected_measure_vertex_no_radius(Vertex_descriptor v) { Vertex_measures vertex_measures; + // add the measures of the faces incident to v (excluding the null (boundary) face) for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { if (f == boost::graph_traits::null_face()) continue; + // only add the measures for the selected curvatures (area measure is always added) vertex_measures.area_measure += get(mu0_map, f); if (is_mean_curvature_selected) @@ -881,8 +915,10 @@ private: return vertex_measures; } + // expand the measures of the faces inside the ball of radius r around v Vertex_measures expand_interpolated_corrected_measure_vertex(Vertex_descriptor v) { + // the ball expansion is done using a BFS traversal from the vertex std::queue bfs_queue; std::unordered_set bfs_visited; @@ -891,6 +927,7 @@ private: Vertex_measures vertex_measures; + // add the measures of the faces incident to v (excluding the null (boundary) face) for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { if (f != boost::graph_traits::null_face()) { @@ -910,8 +947,11 @@ private: x.push_back(Vector_3(pi.x(), pi.y(), pi.z())); } + // compute the inclusion ratio of the face in the ball const FT f_ratio = face_in_ball_ratio(x, ball_radius, c); + // if the face is inside the ball, add the measures + // only add the measures for the selected curvatures (area measure is always added) if (f_ratio != 0.0) { vertex_measures.area_measure += f_ratio * get(mu0_map, fi); @@ -947,10 +987,13 @@ private: for (Vertex_descriptor v : vertices(pmesh)) { + // expand the computed measures (on faces) to the vertices Vertex_measures vertex_measures = (ball_radius < 0) ? expand_interpolated_corrected_measure_vertex_no_radius(v) : expand_interpolated_corrected_measure_vertex(v); + // compute the selected curvatures from the expanded measures and store them in the property maps + // if the area measure is zero, the curvature is set to zero if (is_mean_curvature_selected) { vertex_measures.area_measure != 0 ? put(mean_curvature_map, v, 0.5 * vertex_measures.mean_curvature_measure / vertex_measures.area_measure) : @@ -964,6 +1007,7 @@ private: } if (is_principal_curvatures_and_directions_selected) { + // compute the principal curvatures and directions from the anisotropic measure const Vector_3 v_normal = get(vnm, v); const Principal_curvatures_and_directions principal_curvatures_and_directions = principal_curvatures_and_directions_from_anisotropic_measures( vertex_measures.anisotropic_measure, @@ -1402,7 +1446,6 @@ template::vertex_descriptor v, const NamedParameters& np = parameters::default_values()) { - // use interpolated_corrected_curvatures_one_vertex to compute mean curvature typename GT::FT mean_curvature; interpolated_corrected_curvatures_one_vertex(pmesh, v, np.vertex_mean_curvature(&mean_curvature)); return mean_curvature; @@ -1468,7 +1511,6 @@ template::vertex_descriptor v, const NamedParameters& np = parameters::default_values()) { - // use interpolated_corrected_curvatures_one_vertex to compute gaussian curvature typename GT::FT gc; interpolated_corrected_curvatures_one_vertex(pmesh, v, np.vertex_gaussian_curvature(&gc)); return gc; @@ -1533,7 +1575,6 @@ template::vertex_descriptor v, const NamedParameters& np = parameters::default_values()) { - // use interpolated_corrected_curvatures_one_vertex to compute principal curvatures Principal_curvatures_and_directions pcd; interpolated_corrected_curvatures_one_vertex(pmesh, v, np.vertex_principal_curvatures_and_directions(&pcd)); return pcd;