diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Curvatures/interpolated_corrected_curvature_measures.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Curvatures/interpolated_corrected_curvature_measures.h index 4a50a618e66..94aaff9b8b5 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Curvatures/interpolated_corrected_curvature_measures.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Curvatures/interpolated_corrected_curvature_measures.h @@ -79,675 +79,911 @@ struct Principal_curvatures_and_directions { namespace internal { - template - typename GT::FT average_edge_length(const PolygonMesh& pmesh) { - const std::size_t n = edges(pmesh).size(); - if (n == 0) - return 0; +template +typename GT::FT average_edge_length(const PolygonMesh& pmesh) { + const std::size_t n = edges(pmesh).size(); + if (n == 0) + return 0; - typename GT::FT avg_edge_length = 0; - for (auto e : edges(pmesh)) - avg_edge_length += edge_length(e, pmesh); + typename GT::FT avg_edge_length = 0; + for (auto e : edges(pmesh)) + avg_edge_length += edge_length(e, pmesh); - avg_edge_length /= n; - return avg_edge_length; - } + avg_edge_length /= n; + return avg_edge_length; +} - enum Curvature_measure_index { - MU0_AREA_MEASURE, ///< corrected area density - MU1_MEAN_CURVATURE_MEASURE, ///< corrected mean curvature density - MU2_GAUSSIAN_CURVATURE_MEASURE ///< corrected gaussian curvature density - }; +enum Curvature_measure_index { + MU0_AREA_MEASURE, ///< corrected area density + MU1_MEAN_CURVATURE_MEASURE, ///< corrected mean curvature density + MU2_GAUSSIAN_CURVATURE_MEASURE ///< corrected gaussian curvature density +}; - template - struct Vertex_measures { - typename GT::FT area_measure = 0; - typename GT::FT mean_curvature_measure = 0; - typename GT::FT gaussian_curvature_measure = 0; - std::array anisotropic_measure = { 0, 0, 0, - 0, 0, 0, - 0, 0, 0 }; - }; +template +struct Vertex_measures { + typename GT::FT area_measure = 0; + typename GT::FT mean_curvature_measure = 0; + typename GT::FT gaussian_curvature_measure = 0; + std::array anisotropic_measure = { 0, 0, 0, + 0, 0, 0, + 0, 0, 0 }; +}; - template - typename GT::FT interpolated_corrected_area_measure_face(const std::vector& u, - const std::vector& x) +template +typename GT::FT interpolated_corrected_area_measure_face(const std::vector& u, + const std::vector& x) +{ + const std::size_t n = x.size(); + CGAL_precondition(u.size() == n); + CGAL_precondition(n >= 3); + + typename GT::Construct_cross_product_vector_3 cross_product; + + // Triangle: use triangle formula + if (n == 3) { - const std::size_t n = x.size(); - CGAL_precondition(u.size() == n); - CGAL_precondition(n >= 3); - - typename GT::Construct_cross_product_vector_3 cross_product; - - // Triangle: use triangle formula - if (n == 3) - { - const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; - return 0.5 * um * cross_product(x[1] - x[0], x[2] - x[0]); - } - // Quad: use bilinear interpolation formula - else if (n == 4) - { - // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. - // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 - return (1.0 / 36.0) * ( - (4 * u[0] + 2 * u[1] + 2 * u[3] + u[2]) * cross_product(x[1] - x[0], x[3] - x[0]) - + (2 * u[0] + 4 * u[1] + u[3] + 2 * u[2]) * cross_product(x[1] - x[0], x[2] - x[1]) - + (2 * u[0] + u[1] + 4 * u[3] + 2 * u[2]) * cross_product(x[2] - x[3], x[3] - x[0]) - + (u[0] + 2 * u[1] + 2 * u[3] + 4 * u[2]) * cross_product(x[2] - x[3], x[2] - x[1]) - ); - } - // N-gon: split into n triangles by polygon center and use triangle formula for each - else - { - typename GT::FT mu0 = 0; - - // getting center of points - typename GT::Vector_3 xc = - std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); - xc /= n; - - // getting unit average normal of points - typename GT::Vector_3 uc = - std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); - uc /= sqrt(uc * uc); - - // summing each triangle's measure after triangulation by barycenter split. - for (std::size_t i = 0; i < n; i++) - { - mu0 += interpolated_corrected_area_measure_face( - std::vector {u[i], u[(i + 1) % n], uc}, - std::vector {x[i], x[(i + 1) % n], xc} - ); - } - return mu0; - } + const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; + return 0.5 * um * cross_product(x[1] - x[0], x[2] - x[0]); } - - template - typename GT::FT interpolated_corrected_mean_curvature_measure_face(const std::vector& u, - const std::vector& x) + // Quad: use bilinear interpolation formula + else if (n == 4) { - const std::size_t n = x.size(); - CGAL_precondition(u.size() == n); - CGAL_precondition(n >= 3); - - typename GT::Construct_cross_product_vector_3 cross_product; - - // Triangle: use triangle formula - if (n == 3) - { - const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; - - return 0.5 * um * (cross_product(u[2] - u[1], x[0]) - + cross_product(u[0] - u[2], x[1]) - + cross_product(u[1] - u[0], x[2])); - } - // Quad: use bilinear interpolation formula - else if (n == 4) - { - // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. - // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 - - const typename GT::Vector_3 u02 = u[2] - u[0]; - const typename GT::Vector_3 u13 = u[3] - u[1]; - const typename GT::Vector_3 x0_cross = cross_product(u13, x[0]); - const typename GT::Vector_3 x1_cross = -cross_product(u02, x[1]); - const typename GT::Vector_3 x3_cross = cross_product(u02, x[3]); - const typename GT::Vector_3 x2_cross = -cross_product(u13, x[2]); - - return (1.0 / 12.0) * ( - u[0] * (2 * x0_cross - cross_product((u[3] + u[2]), x[1]) + cross_product((u[1] + u[2]), x[3]) + x2_cross) - + u[1] * (cross_product((u[3] + u[2]), x[0]) + 2 * x1_cross + x3_cross - cross_product((u[0] + u[3]), x[2])) - + u[3] * (-cross_product((u[1] + u[2]), x[0]) + x1_cross + 2 * x3_cross + cross_product((u[0] + u[1]), x[2])) - + u[2] * (x0_cross + cross_product((u[0] + u[3]), x[1]) - cross_product((u[0] + u[1]), x[3]) + 2 * x2_cross) - ); - } - // N-gon: split into n triangles by polygon center and use triangle formula for each - else - { - typename GT::FT mu1 = 0; - - // getting center of points - typename GT::Vector_3 xc = - std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); - xc /= n; - - // getting unit average normal of points - typename GT::Vector_3 uc = - std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); - uc /= sqrt(uc * uc); - - // summing each triangle's measure after triangulation by barycenter split. - for (std::size_t i = 0; i < n; i++) - { - mu1 += interpolated_corrected_mean_curvature_measure_face( - std::vector {u[i], u[(i + 1) % n], uc}, - std::vector {x[i], x[(i + 1) % n], xc} - ); - } - return mu1; - } - } - - template - typename GT::FT interpolated_corrected_gaussian_curvature_measure_face(const std::vector& u, - const std::vector& x = {}) - { - const std::size_t n = u.size(); - CGAL_precondition(n >= 3); - - typename GT::Construct_cross_product_vector_3 cross_product; - - // Triangle: use triangle formula - if (n == 3) - { - return 0.5 * u[0] * cross_product(u[1], u[2]); - } - // Quad: use bilinear interpolation formula - else if (n == 4) - { - // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. - // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 - return (1.0 / 36.0) * ( - (4 * u[0] + 2 * u[1] + 2 * u[3] + u[2]) * cross_product(u[1] - u[0], u[3] - u[0]) - + (2 * u[0] + 4 * u[1] + u[3] + 2 * u[2]) * cross_product(u[1] - u[0], u[2] - u[1]) - + (2 * u[0] + u[1] + 4 * u[3] + 2 * u[2]) * cross_product(u[2] - u[3], u[3] - u[0]) - + (u[0] + 2 * u[1] + 2 * u[3] + 4 * u[2]) * cross_product(u[2] - u[3], u[2] - u[1]) - ); - } - // N-gon: split into n triangles by polygon center and use triangle formula for each - else - { - typename GT::FT mu2 = 0; - - // getting unit average normal of points - typename GT::Vector_3 uc = - std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); - uc /= sqrt(uc * uc); - - // summing each triangle's measure after triangulation by barycenter split. - for (std::size_t i = 0; i < n; i++) - { - mu2 += interpolated_corrected_gaussian_curvature_measure_face( - std::vector {u[i], u[(i + 1) % n], uc} - ); - } - return mu2; - } - } - - template - std::array interpolated_corrected_anisotropic_measure_face(const std::vector& u, - const std::vector& x) - { - const std::size_t n = x.size(); - CGAL_precondition(u.size() == n); - CGAL_precondition(n >= 3); - - typename GT::Construct_cross_product_vector_3 cross_product; - std::array muXY{ 0 }; - - // Triangle: use triangle formula - if (n == 3) - { - const typename GT::Vector_3 u01 = u[1] - u[0]; - const typename GT::Vector_3 u02 = u[2] - u[0]; - const typename GT::Vector_3 x01 = x[1] - x[0]; - const typename GT::Vector_3 x02 = x[2] - x[0]; - const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; - - for (std::size_t ix = 0; ix < 3; ix++) - { - typename GT::Vector_3 X; - if (ix == 0) - X = typename GT::Vector_3(1, 0, 0); - if (ix == 1) - X = typename GT::Vector_3(0, 1, 0); - if (ix == 2) - X = typename GT::Vector_3(0, 0, 1); - - for (std::size_t iy = 0; iy < 3; iy++) - muXY[ix * 3 + iy] = 0.5 * um * (cross_product(u02[iy] * X, x01) - cross_product(u01[iy] * X, x02)); - } - } - // Quad: use bilinear interpolation formula - else if (n == 4) - { - // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. - // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 - for (std::size_t ix = 0; ix < 3; ix++) - { - typename GT::Vector_3 X; - if (ix == 0) - X = typename GT::Vector_3(1, 0, 0); - if (ix == 1) - X = typename GT::Vector_3(0, 1, 0); - if (ix == 2) - X = typename GT::Vector_3(0, 0, 1); - - const typename GT::Vector_3 u0xX = cross_product(u[0], X); - const typename GT::Vector_3 u1xX = cross_product(u[1], X); - const typename GT::Vector_3 u2xX = cross_product(u[2], X); - const typename GT::Vector_3 u3xX = cross_product(u[3], X); - - for (std::size_t iy = 0; iy < 3; iy++) - muXY[ix * 3 + iy] = (1.0 / 72.0) * ( - - u[0][iy] * (u0xX * (-x[0] - 11 * x[1] + 13 * x[3] - x[2]) - + u1xX * (-5 * x[0] - 7 * x[1] + 11 * x[3] + x[2]) - + u3xX * (x[0] - 7 * x[1] + 11 * x[3] - 5 * x[2]) - + u2xX * (-x[0] - 5 * x[1] + 7 * x[3] - x[2]) - ) - + u[1][iy] * (u0xX * (13 * x[0] - x[1] - 7 * x[3] - 5 * x[2]) - + u1xX * (17 * x[0] - 5 * x[1] - 5 * x[3] - 7 * x[2]) - + u3xX * (5 * x[0] + x[1] + x[3] - 7 * x[2]) - + u2xX * (7 * x[0] - x[1] + 5 * x[3] - 11 * x[2]) - ) - + u[2][iy] * (u0xX * (-11 * x[0] + 5 * x[1] - x[3] + 7 * x[2]) - + u1xX * (-7 * x[0] + x[1] + x[3] + 5 * x[2]) - + u3xX * (-7 * x[0] - 5 * x[1] - 5 * x[3] + 17 * x[2]) - + u2xX * (-5 * x[0] - 7 * x[1] - x[3] + 13 * x[2]) - ) - + u[3][iy] * (u0xX * (-x[0] + 7 * x[1] - 5 * x[3] - x[2]) - + u1xX * (-5 * x[0] + 11 * x[1] - 7 * x[3] + x[2]) - + u3xX * (x[0] + 11 * x[1] - 7 * x[3] - 5 * x[2]) - + u2xX * (-x[0] + 13 * x[1] - 11 * x[3] - x[2]) - ) - - ); - } - } - // N-gon: split into n triangles by polygon center and use triangle formula for each - else - { - // getting center of points - typename GT::Vector_3 xc = - std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); - xc /= n; - - // getting unit average normal of points - typename GT::Vector_3 uc = - std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); - uc /= sqrt(uc * uc); - - // summing each triangle's measure after triangulation by barycenter split. - for (std::size_t i = 0; i < n; i++) - { - std::array muXY_curr_triangle = - interpolated_corrected_anisotropic_measure_face( - std::vector {u[i], u[(i + 1) % n], uc}, - std::vector {x[i], x[(i + 1) % n], xc} - ); - - for (std::size_t ix = 0; ix < 3; ix++) - for (std::size_t iy = 0; iy < 3; iy++) - muXY[ix * 3 + iy] += muXY_curr_triangle[ix * 3 + iy]; - } - } - return muXY; - } - - //template - //typename GT::FT triangle_in_ball_ratio(const typename GT::Vector_3 x1, - // const typename GT::Vector_3 x2, - // const typename GT::Vector_3 x3, - // const typename GT::FT r, - // const typename GT::Vector_3 c, - // const std::size_t res = 3) - //{ - // const typename GT::FT R = r * r; - // const typename GT::FT acc = 1.0 / res; - // std::size_t samples_in = 0; - // for (GT::FT alpha = acc / 3; alpha < 1; alpha += acc) - // for (GT::FT beta = acc / 3; beta < 1 - alpha; beta += acc) - // { - // if ((alpha * x1 + beta * x2 + (1 - alpha - beta) * x3 - c).squared_length() < R) - // samples_in++; - // } - // return samples_in / (typename GT::FT)(res * (res + 1) / 2); - //} - - template - typename GT::FT face_in_ball_ratio(const std::vector& x, - const typename GT::FT r, - const typename GT::Vector_3 c) - { - const std::size_t n = x.size(); - - // getting center of points - typename GT::Vector_3 xm = - std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); - xm /= n; - - typename GT::FT d_min = (xm - c).squared_length(); - typename GT::FT d_max = d_min; - - for (const typename GT::Vector_3 xi : x) - { - const typename GT::FT d_sq = (xi - c).squared_length(); - d_max = (std::max)(d_sq, d_max); - d_min = (std::min)(d_sq, d_min); - } - - if (d_max <= r * r) return 1.0; - else if (r * r <= d_min) return 0.0; - - d_max = sqrt(d_max); - d_min = sqrt(d_min); - - return (r - d_min) / (d_max - d_min); - } - - template - Principal_curvatures_and_directions principal_curvatures_and_directions_from_anisotropic_measures( - const std::array anisotropic_measure, - const typename GT::FT v_mu0, - const typename GT::Vector_3 u_GT - ) - { - 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]; - - Eigen::Matrix u(u_GT.x(), u_GT.y(), u_GT.z()); - const typename GT::FT K = 1000 * v_mu0; - - v_muXY = 0.5 * (v_muXY + v_muXY.transpose()) + K * u * u.transpose(); - - Eigen::SelfAdjointEigenSolver> eigensolver; - - eigensolver.computeDirect(v_muXY); - - if (eigensolver.info() != Eigen::Success) - return Principal_curvatures_and_directions(); - - const Eigen::Matrix eig_vals = eigensolver.eigenvalues(); - const Eigen::Matrix eig_vecs = eigensolver.eigenvectors(); - - 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)); - - 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, - min_eig_vec, - max_eig_vec + // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. + // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 + return (1.0 / 36.0) * ( + (4 * u[0] + 2 * u[1] + 2 * u[3] + u[2]) * cross_product(x[1] - x[0], x[3] - x[0]) + + (2 * u[0] + 4 * u[1] + u[3] + 2 * u[2]) * cross_product(x[1] - x[0], x[2] - x[1]) + + (2 * u[0] + u[1] + 4 * u[3] + 2 * u[2]) * cross_product(x[2] - x[3], x[3] - x[0]) + + (u[0] + 2 * u[1] + 2 * u[3] + 4 * u[2]) * cross_product(x[2] - x[3], x[2] - x[1]) ); } - - template - class Interpolated_corrected_curvatures_computer + // N-gon: split into n triangles by polygon center and use triangle formula for each + else { - typedef typename GetGeomTraits::type GT; + typename GT::FT mu0 = 0; - typedef typename GT::FT FT; - typedef typename GT::Point_3 Point_3; - typedef typename GT::Vector_3 Vector_3; + // getting center of points + typename GT::Vector_3 xc = + std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); + xc /= n; - typedef typename boost::graph_traits::halfedge_descriptor Halfedge_descriptor; - typedef typename boost::graph_traits::edge_descriptor Edge_descriptor; - typedef typename boost::graph_traits::face_descriptor Face_descriptor; - typedef typename boost::graph_traits::vertex_descriptor Vertex_descriptor; + // getting unit average normal of points + typename GT::Vector_3 uc = + std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); + uc /= sqrt(uc * uc); - typedef typename GetVertexPointMap::const_type Vertex_position_map; - - typedef dynamic_vertex_property_t Vector_map_tag; - typedef typename boost::property_map::const_type Default_vector_map; - typedef typename internal_np::Lookup_named_param_def::type Vertex_normal_map; - - typedef Constant_property_map Default_scalar_map; - typedef typename internal_np::Lookup_named_param_def::type Vertex_mean_curvature_map; - - typedef typename internal_np::Lookup_named_param_def::type Vertex_gaussian_curvature_map; - - typedef Constant_property_map> Default_principal_map; - typedef typename internal_np::Lookup_named_param_def::type Vertex_principal_curvatures_and_directions_map; - - typedef typename boost::property_map>::const_type Face_scalar_measure_map; - typedef typename boost::property_map>>::const_type Face_anisotropic_measure_map; - - private: - const PolygonMesh& pmesh; - Vertex_position_map vpm; - Vertex_normal_map vnm; - FT ball_radius; - - bool is_mean_curvature_selected; - bool is_gaussian_curvature_selected; - bool is_principal_curvatures_and_directions_selected; - - Vertex_mean_curvature_map mean_curvature_map; - Vertex_gaussian_curvature_map gaussian_curvature_map; - Vertex_principal_curvatures_and_directions_map principal_curvatures_and_directions_map; - - Face_scalar_measure_map mu0_map, mu1_map, mu2_map; - Face_anisotropic_measure_map muXY_map; - - void set_property_maps() { - mu0_map = get(CGAL::dynamic_face_property_t(), pmesh); - mu1_map = get(CGAL::dynamic_face_property_t(), pmesh); - mu2_map = get(CGAL::dynamic_face_property_t(), pmesh); - muXY_map = get(CGAL::dynamic_face_property_t>(), pmesh); - - } - - void set_named_params(const NamedParameters& np) + // summing each triangle's measure after triangulation by barycenter split. + for (std::size_t i = 0; i < n; i++) { - using parameters::choose_parameter; - using parameters::get_parameter; - using parameters::is_default_parameter; - - vpm = choose_parameter(get_parameter(np, CGAL::vertex_point), - get_const_property_map(CGAL::vertex_point, pmesh)); - - vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map), - get(Vector_map_tag(), pmesh)); - - if (is_default_parameter::value) - compute_vertex_normals(pmesh, vnm, np); - - const FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1); - - 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; - - mean_curvature_map = choose_parameter(get_parameter(np, internal_np::vertex_mean_curvature_map), Default_scalar_map()); - gaussian_curvature_map = choose_parameter(get_parameter(np, internal_np::vertex_gaussian_curvature_map), Default_scalar_map()); - principal_curvatures_and_directions_map = choose_parameter(get_parameter(np, internal_np::vertex_principal_curvatures_and_directions_map), Default_principal_map()); - - set_ball_radius(radius); + mu0 += interpolated_corrected_area_measure_face( + std::vector {u[i], u[(i + 1) % n], uc}, + std::vector {x[i], x[(i + 1) % n], xc} + ); } + return mu0; + } +} - void set_ball_radius(const FT radius) { - if (radius == 0) - ball_radius = average_edge_length(pmesh) * EXPANDING_RADIUS_EPSILON; - else - ball_radius = radius; - } +template +typename GT::FT interpolated_corrected_mean_curvature_measure_face(const std::vector& u, + const std::vector& x) +{ + const std::size_t n = x.size(); + CGAL_precondition(u.size() == n); + CGAL_precondition(n >= 3); - public: + typename GT::Construct_cross_product_vector_3 cross_product; - Interpolated_corrected_curvatures_computer(const PolygonMesh& pmesh, - const NamedParameters& np = parameters::default_values() - ) : - pmesh(pmesh) + // Triangle: use triangle formula + if (n == 3) + { + const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; + + return 0.5 * um * (cross_product(u[2] - u[1], x[0]) + + cross_product(u[0] - u[2], x[1]) + + cross_product(u[1] - u[0], x[2])); + } + // Quad: use bilinear interpolation formula + else if (n == 4) + { + // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. + // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 + + const typename GT::Vector_3 u02 = u[2] - u[0]; + const typename GT::Vector_3 u13 = u[3] - u[1]; + const typename GT::Vector_3 x0_cross = cross_product(u13, x[0]); + const typename GT::Vector_3 x1_cross = -cross_product(u02, x[1]); + const typename GT::Vector_3 x3_cross = cross_product(u02, x[3]); + const typename GT::Vector_3 x2_cross = -cross_product(u13, x[2]); + + return (1.0 / 12.0) * ( + u[0] * (2 * x0_cross - cross_product((u[3] + u[2]), x[1]) + cross_product((u[1] + u[2]), x[3]) + x2_cross) + + u[1] * (cross_product((u[3] + u[2]), x[0]) + 2 * x1_cross + x3_cross - cross_product((u[0] + u[3]), x[2])) + + u[3] * (-cross_product((u[1] + u[2]), x[0]) + x1_cross + 2 * x3_cross + cross_product((u[0] + u[1]), x[2])) + + u[2] * (x0_cross + cross_product((u[0] + u[3]), x[1]) - cross_product((u[0] + u[1]), x[3]) + 2 * x2_cross) + ); + } + // N-gon: split into n triangles by polygon center and use triangle formula for each + else + { + typename GT::FT mu1 = 0; + + // getting center of points + typename GT::Vector_3 xc = + std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); + xc /= n; + + // getting unit average normal of points + typename GT::Vector_3 uc = + std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); + uc /= sqrt(uc * uc); + + // summing each triangle's measure after triangulation by barycenter split. + for (std::size_t i = 0; i < n; i++) { - set_named_params(np); + mu1 += interpolated_corrected_mean_curvature_measure_face( + std::vector {u[i], u[(i + 1) % n], uc}, + std::vector {x[i], x[(i + 1) % n], xc} + ); + } + return mu1; + } +} - if (is_mean_curvature_selected || is_gaussian_curvature_selected || is_principal_curvatures_and_directions_selected) +template +typename GT::FT interpolated_corrected_gaussian_curvature_measure_face(const std::vector& u, + const std::vector& x = {}) +{ + const std::size_t n = u.size(); + CGAL_precondition(n >= 3); + + typename GT::Construct_cross_product_vector_3 cross_product; + + // Triangle: use triangle formula + if (n == 3) + { + return 0.5 * u[0] * cross_product(u[1], u[2]); + } + // Quad: use bilinear interpolation formula + else if (n == 4) + { + // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. + // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 + return (1.0 / 36.0) * ( + (4 * u[0] + 2 * u[1] + 2 * u[3] + u[2]) * cross_product(u[1] - u[0], u[3] - u[0]) + + (2 * u[0] + 4 * u[1] + u[3] + 2 * u[2]) * cross_product(u[1] - u[0], u[2] - u[1]) + + (2 * u[0] + u[1] + 4 * u[3] + 2 * u[2]) * cross_product(u[2] - u[3], u[3] - u[0]) + + (u[0] + 2 * u[1] + 2 * u[3] + 4 * u[2]) * cross_product(u[2] - u[3], u[2] - u[1]) + ); + } + // N-gon: split into n triangles by polygon center and use triangle formula for each + else + { + typename GT::FT mu2 = 0; + + // getting unit average normal of points + typename GT::Vector_3 uc = + std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); + uc /= sqrt(uc * uc); + + // summing each triangle's measure after triangulation by barycenter split. + for (std::size_t i = 0; i < n; i++) + { + mu2 += interpolated_corrected_gaussian_curvature_measure_face( + std::vector {u[i], u[(i + 1) % n], uc} + ); + } + return mu2; + } +} + +template +std::array interpolated_corrected_anisotropic_measure_face(const std::vector& u, + const std::vector& x) +{ + const std::size_t n = x.size(); + CGAL_precondition(u.size() == n); + CGAL_precondition(n >= 3); + + typename GT::Construct_cross_product_vector_3 cross_product; + std::array muXY{ 0 }; + + // Triangle: use triangle formula + if (n == 3) + { + const typename GT::Vector_3 u01 = u[1] - u[0]; + const typename GT::Vector_3 u02 = u[2] - u[0]; + const typename GT::Vector_3 x01 = x[1] - x[0]; + const typename GT::Vector_3 x02 = x[2] - x[0]; + const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; + + for (std::size_t ix = 0; ix < 3; ix++) + { + typename GT::Vector_3 X; + if (ix == 0) + X = typename GT::Vector_3(1, 0, 0); + if (ix == 1) + X = typename GT::Vector_3(0, 1, 0); + if (ix == 2) + X = typename GT::Vector_3(0, 0, 1); + + for (std::size_t iy = 0; iy < 3; iy++) + muXY[ix * 3 + iy] = 0.5 * um * (cross_product(u02[iy] * X, x01) - cross_product(u01[iy] * X, x02)); + } + } + // Quad: use bilinear interpolation formula + else if (n == 4) + { + // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. + // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 + for (std::size_t ix = 0; ix < 3; ix++) + { + typename GT::Vector_3 X; + if (ix == 0) + X = typename GT::Vector_3(1, 0, 0); + if (ix == 1) + X = typename GT::Vector_3(0, 1, 0); + if (ix == 2) + X = typename GT::Vector_3(0, 0, 1); + + const typename GT::Vector_3 u0xX = cross_product(u[0], X); + const typename GT::Vector_3 u1xX = cross_product(u[1], X); + const typename GT::Vector_3 u2xX = cross_product(u[2], X); + const typename GT::Vector_3 u3xX = cross_product(u[3], X); + + for (std::size_t iy = 0; iy < 3; iy++) + muXY[ix * 3 + iy] = (1.0 / 72.0) * ( + + u[0][iy] * (u0xX * (-x[0] - 11 * x[1] + 13 * x[3] - x[2]) + + u1xX * (-5 * x[0] - 7 * x[1] + 11 * x[3] + x[2]) + + u3xX * (x[0] - 7 * x[1] + 11 * x[3] - 5 * x[2]) + + u2xX * (-x[0] - 5 * x[1] + 7 * x[3] - x[2]) + ) + + u[1][iy] * (u0xX * (13 * x[0] - x[1] - 7 * x[3] - 5 * x[2]) + + u1xX * (17 * x[0] - 5 * x[1] - 5 * x[3] - 7 * x[2]) + + u3xX * (5 * x[0] + x[1] + x[3] - 7 * x[2]) + + u2xX * (7 * x[0] - x[1] + 5 * x[3] - 11 * x[2]) + ) + + u[2][iy] * (u0xX * (-11 * x[0] + 5 * x[1] - x[3] + 7 * x[2]) + + u1xX * (-7 * x[0] + x[1] + x[3] + 5 * x[2]) + + u3xX * (-7 * x[0] - 5 * x[1] - 5 * x[3] + 17 * x[2]) + + u2xX * (-5 * x[0] - 7 * x[1] - x[3] + 13 * x[2]) + ) + + u[3][iy] * (u0xX * (-x[0] + 7 * x[1] - 5 * x[3] - x[2]) + + u1xX * (-5 * x[0] + 11 * x[1] - 7 * x[3] + x[2]) + + u3xX * (x[0] + 11 * x[1] - 7 * x[3] - 5 * x[2]) + + u2xX * (-x[0] + 13 * x[1] - 11 * x[3] - x[2]) + ) + + ); + } + } + // N-gon: split into n triangles by polygon center and use triangle formula for each + else + { + // getting center of points + typename GT::Vector_3 xc = + std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); + xc /= n; + + // getting unit average normal of points + typename GT::Vector_3 uc = + std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); + uc /= sqrt(uc * uc); + + // summing each triangle's measure after triangulation by barycenter split. + for (std::size_t i = 0; i < n; i++) + { + std::array muXY_curr_triangle = + interpolated_corrected_anisotropic_measure_face( + std::vector {u[i], u[(i + 1) % n], uc}, + std::vector {x[i], x[(i + 1) % n], xc} + ); + + for (std::size_t ix = 0; ix < 3; ix++) + for (std::size_t iy = 0; iy < 3; iy++) + muXY[ix * 3 + iy] += muXY_curr_triangle[ix * 3 + iy]; + } + } + return muXY; +} + +//template +//typename GT::FT triangle_in_ball_ratio(const typename GT::Vector_3 x1, +// const typename GT::Vector_3 x2, +// const typename GT::Vector_3 x3, +// const typename GT::FT r, +// const typename GT::Vector_3 c, +// const std::size_t res = 3) +//{ +// const typename GT::FT R = r * r; +// const typename GT::FT acc = 1.0 / res; +// std::size_t samples_in = 0; +// for (GT::FT alpha = acc / 3; alpha < 1; alpha += acc) +// for (GT::FT beta = acc / 3; beta < 1 - alpha; beta += acc) +// { +// if ((alpha * x1 + beta * x2 + (1 - alpha - beta) * x3 - c).squared_length() < R) +// samples_in++; +// } +// return samples_in / (typename GT::FT)(res * (res + 1) / 2); +//} + +template +typename GT::FT face_in_ball_ratio(const std::vector& x, + const typename GT::FT r, + const typename GT::Vector_3 c) +{ + const std::size_t n = x.size(); + + // getting center of points + typename GT::Vector_3 xm = + std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); + xm /= n; + + typename GT::FT d_min = (xm - c).squared_length(); + typename GT::FT d_max = d_min; + + for (const typename GT::Vector_3 xi : x) + { + const typename GT::FT d_sq = (xi - c).squared_length(); + d_max = (std::max)(d_sq, d_max); + d_min = (std::min)(d_sq, d_min); + } + + if (d_max <= r * r) return 1.0; + else if (r * r <= d_min) return 0.0; + + d_max = sqrt(d_max); + d_min = sqrt(d_min); + + return (r - d_min) / (d_max - d_min); +} + +template +Principal_curvatures_and_directions principal_curvatures_and_directions_from_anisotropic_measures( + const std::array anisotropic_measure, + const typename GT::FT v_mu0, + const typename GT::Vector_3 u_GT +) +{ + 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]; + + Eigen::Matrix u(u_GT.x(), u_GT.y(), u_GT.z()); + const typename GT::FT K = 1000 * v_mu0; + + v_muXY = 0.5 * (v_muXY + v_muXY.transpose()) + K * u * u.transpose(); + + Eigen::SelfAdjointEigenSolver> eigensolver; + + eigensolver.computeDirect(v_muXY); + + if (eigensolver.info() != Eigen::Success) + return Principal_curvatures_and_directions(); + + const Eigen::Matrix eig_vals = eigensolver.eigenvalues(); + const Eigen::Matrix eig_vecs = eigensolver.eigenvectors(); + + 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)); + + 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, + min_eig_vec, + max_eig_vec + ); +} + +template +typename Vertex_measures interpolated_corrected_measures_one_vertex_no_radius( + const PolygonMesh pmesh, + const typename boost::graph_traits::vertex_descriptor v, + const bool is_mean_curvature_selected, + const bool is_gaussian_curvature_selected, + const bool is_principal_curvatures_and_directions_selected, + const VPM vpm, + const VNM vnm +) +{ + typedef typename boost::graph_traits::face_descriptor Face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor Vertex_descriptor; + typedef typename GT::Point_3 Point_3; + typedef typename GT::Vector_3 Vector_3; + typedef typename GT::FT FT; + + std::queue bfs_queue; + std::unordered_set bfs_visited; + + typename Vertex_measures vertex_measures; + + std::vector x; + std::vector u; + + for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { + if (f != boost::graph_traits::null_face()) + { + for (Vertex_descriptor vi : vertices_around_face(halfedge(f, pmesh), pmesh)) { - set_property_maps(); + Point_3 pi = get(vpm, vi); + Vector_3 ui = get(vnm, vi); + x.push_back(Vector_3(pi.x(), pi.y(), pi.z())); + u.push_back(ui); + } - compute_selected_curvatures(); + vertex_measures.area_measure += interpolated_corrected_area_measure_face(u, x); + + if (is_mean_curvature_selected) + vertex_measures.mean_curvature_measure += interpolated_corrected_mean_curvature_measure_face(u, x); + + if (is_gaussian_curvature_selected) + vertex_measures.gaussian_curvature_measure += interpolated_corrected_gaussian_curvature_measure_face(u, x); + + if (is_principal_curvatures_and_directions_selected) + { + const std::array face_anisotropic_measure = interpolated_corrected_anisotropic_measure_face(u, x); + for (std::size_t i = 0; i < 3 * 3; i++) + vertex_measures.anisotropic_measure[i] += face_anisotropic_measure[i]; } } - private: + x.clear(); + u.clear(); + } + return vertex_measures; +} - void interpolated_corrected_selected_measures_all_faces() + +template +typename Vertex_measures interpolated_corrected_measures_one_vertex( + const PolygonMesh pmesh, + const typename boost::graph_traits::vertex_descriptor v, + const typename GT::FT radius, + const bool is_mean_curvature_selected, + const bool is_gaussian_curvature_selected, + const bool is_principal_curvatures_and_directions_selected, + const VPM vpm, + const VNM vnm +) +{ + typedef typename boost::graph_traits::face_descriptor Face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor Vertex_descriptor; + typedef typename GT::Point_3 Point_3; + typedef typename GT::Vector_3 Vector_3; + typedef typename GT::FT FT; + + std::queue bfs_queue; + std::unordered_set bfs_visited; + + typename Vertex_measures vertex_measures; + + typename Point_3 vp = get(vpm, v); + typename Vector_3 c = Vector_3(vp.x(), vp.y(), vp.z()); + + for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { + if (f != boost::graph_traits::null_face()) { + bfs_queue.push(f); + bfs_visited.insert(f); + } + } + std::vector x; + std::vector u; + while (!bfs_queue.empty()) { + Face_descriptor fi = bfs_queue.front(); + bfs_queue.pop(); + + // looping over vertices in face to get point coordinates + + for (Vertex_descriptor vi : vertices_around_face(halfedge(fi, pmesh), pmesh)) + { + Point_3 pi = get(vpm, vi); + Vector_3 ui = get(vnm, vi); + x.push_back(Vector_3(pi.x(), pi.y(), pi.z())); + u.push_back(ui); + } + + const FT f_ratio = face_in_ball_ratio(x, radius, c); + + if (f_ratio != 0.0) + { + vertex_measures.area_measure += f_ratio * interpolated_corrected_area_measure_face(u, x); + + if (is_mean_curvature_selected) + vertex_measures.mean_curvature_measure += f_ratio * interpolated_corrected_mean_curvature_measure_face(u, x); + + if (is_gaussian_curvature_selected) + vertex_measures.gaussian_curvature_measure += f_ratio * interpolated_corrected_gaussian_curvature_measure_face(u, x); + + if (is_principal_curvatures_and_directions_selected) + { + const std::array face_anisotropic_measure = interpolated_corrected_anisotropic_measure_face(u, x); + for (std::size_t i = 0; i < 3 * 3; i++) + vertex_measures.anisotropic_measure[i] += f_ratio * face_anisotropic_measure[i]; + } + + for (Face_descriptor fj : faces_around_face(halfedge(fi, pmesh), pmesh)) + { + if (bfs_visited.find(fj) == bfs_visited.end() && fj != boost::graph_traits::null_face()) + { + bfs_queue.push(fj); + bfs_visited.insert(fj); + } + } + } + + x.clear(); + u.clear(); + } + return vertex_measures; +} + +template + void interpolated_corrected_curvatures_one_vertex( + const PolygonMesh pmesh, + const typename boost::graph_traits::vertex_descriptor v, + NamedParameters& np = parameters::default_values() + ) +{ + typedef typename GetGeomTraits::type GT; + typedef typename GetVertexPointMap::const_type Vertex_position_map; + + typedef dynamic_vertex_property_t Vector_map_tag; + typedef typename boost::property_map::const_type Default_vector_map; + typedef typename internal_np::Lookup_named_param_def::type Vertex_normal_map; + + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::is_default_parameter; + + Vertex_position_map vpm = choose_parameter(get_parameter(np, CGAL::vertex_point), + get_const_property_map(CGAL::vertex_point, pmesh)); + + Vertex_normal_map vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map), + get(Vector_map_tag(), pmesh)); + + if (is_default_parameter::value) + compute_vertex_normals(pmesh, vnm, np); + + typename GT::FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1); + if (radius == 0) + radius = average_edge_length(pmesh) * EXPANDING_RADIUS_EPSILON; + else + radius = radius; + + 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); + typename Principal_curvatures_and_directions* vertex_principal_curvatures_and_directions = choose_parameter(get_parameter(np, internal_np::vertex_principal_curvatures_and_directions), nullptr); + + 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); + + std::cout << is_mean_curvature_selected << is_gaussian_curvature_selected << is_principal_curvatures_and_directions_selected << std::endl; + + Vertex_measures vertex_measures; + + if (radius < 0) + vertex_measures = interpolated_corrected_measures_one_vertex_no_radius( + pmesh, + v, + is_mean_curvature_selected, + is_gaussian_curvature_selected, + is_principal_curvatures_and_directions_selected, + vpm, + vnm + ); + else + vertex_measures = interpolated_corrected_measures_one_vertex( + pmesh, + v, + radius, + is_mean_curvature_selected, + is_gaussian_curvature_selected, + is_principal_curvatures_and_directions_selected, + vpm, + vnm + ); + + + if (is_mean_curvature_selected) { + *vertex_mean_curvature = vertex_measures.area_measure != 0 ? + 0.5 * vertex_measures.mean_curvature_measure / vertex_measures.area_measure : 0; + } + + if (is_gaussian_curvature_selected) { + *vertex_gaussian_curvature = vertex_measures.area_measure != 0 ? + vertex_measures.gaussian_curvature_measure / vertex_measures.area_measure : 0; + } + + if (is_principal_curvatures_and_directions_selected) { + const GT::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, + vertex_measures.area_measure, + v_normal + ); + *vertex_principal_curvatures_and_directions = principal_curvatures_and_directions; + } +} + + +template +class Interpolated_corrected_curvatures_computer +{ + typedef typename GetGeomTraits::type GT; + + typedef typename GT::FT FT; + typedef typename GT::Point_3 Point_3; + typedef typename GT::Vector_3 Vector_3; + + typedef typename boost::graph_traits::halfedge_descriptor Halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor Edge_descriptor; + typedef typename boost::graph_traits::face_descriptor Face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor Vertex_descriptor; + + typedef typename GetVertexPointMap::const_type Vertex_position_map; + + typedef dynamic_vertex_property_t Vector_map_tag; + typedef typename boost::property_map::const_type Default_vector_map; + typedef typename internal_np::Lookup_named_param_def::type Vertex_normal_map; + + typedef Constant_property_map Default_scalar_map; + typedef typename internal_np::Lookup_named_param_def::type Vertex_mean_curvature_map; + + typedef typename internal_np::Lookup_named_param_def::type Vertex_gaussian_curvature_map; + + typedef Constant_property_map> Default_principal_map; + typedef typename internal_np::Lookup_named_param_def::type Vertex_principal_curvatures_and_directions_map; + + typedef typename boost::property_map>::const_type Face_scalar_measure_map; + typedef typename boost::property_map>>::const_type Face_anisotropic_measure_map; + +private: + const PolygonMesh& pmesh; + Vertex_position_map vpm; + Vertex_normal_map vnm; + FT ball_radius; + + bool is_mean_curvature_selected; + bool is_gaussian_curvature_selected; + bool is_principal_curvatures_and_directions_selected; + + Vertex_mean_curvature_map mean_curvature_map; + Vertex_gaussian_curvature_map gaussian_curvature_map; + Vertex_principal_curvatures_and_directions_map principal_curvatures_and_directions_map; + + Face_scalar_measure_map mu0_map, mu1_map, mu2_map; + Face_anisotropic_measure_map muXY_map; + + void set_property_maps() { + mu0_map = get(CGAL::dynamic_face_property_t(), pmesh); + mu1_map = get(CGAL::dynamic_face_property_t(), pmesh); + mu2_map = get(CGAL::dynamic_face_property_t(), pmesh); + muXY_map = get(CGAL::dynamic_face_property_t>(), pmesh); + + } + + void set_named_params(const NamedParameters& np) + { + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::is_default_parameter; + + vpm = choose_parameter(get_parameter(np, CGAL::vertex_point), + get_const_property_map(CGAL::vertex_point, pmesh)); + + vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map), + get(Vector_map_tag(), pmesh)); + + if (is_default_parameter::value) + compute_vertex_normals(pmesh, vnm, np); + + const FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1); + + 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; + + mean_curvature_map = choose_parameter(get_parameter(np, internal_np::vertex_mean_curvature_map), Default_scalar_map()); + gaussian_curvature_map = choose_parameter(get_parameter(np, internal_np::vertex_gaussian_curvature_map), Default_scalar_map()); + principal_curvatures_and_directions_map = choose_parameter(get_parameter(np, internal_np::vertex_principal_curvatures_and_directions_map), Default_principal_map()); + + set_ball_radius(radius); + } + + void set_ball_radius(const FT radius) { + if (radius == 0) + ball_radius = average_edge_length(pmesh) * EXPANDING_RADIUS_EPSILON; + else + ball_radius = radius; + } + +public: + + Interpolated_corrected_curvatures_computer(const PolygonMesh& pmesh, + const NamedParameters& np = parameters::default_values() + ) : + pmesh(pmesh) + { + set_named_params(np); + + if (is_mean_curvature_selected || is_gaussian_curvature_selected || is_principal_curvatures_and_directions_selected) + { + set_property_maps(); + + compute_selected_curvatures(); + } + } + +private: + + void interpolated_corrected_selected_measures_all_faces() + { + std::vector x; + std::vector u; + + for (Face_descriptor f : faces(pmesh)) + { + for (Vertex_descriptor v : vertices_around_face(halfedge(f, pmesh), pmesh)) + { + Point_3 p = get(vpm, v); + x.push_back(Vector_3(p.x(), p.y(), p.z())); + u.push_back(get(vnm, v)); + } + put(mu0_map, f, interpolated_corrected_area_measure_face(u, x)); + + if (is_mean_curvature_selected) + put(mu1_map, f, interpolated_corrected_mean_curvature_measure_face(u, x)); + + if (is_gaussian_curvature_selected) + put(mu2_map, f, interpolated_corrected_gaussian_curvature_measure_face(u)); + + if (is_principal_curvatures_and_directions_selected) + put(muXY_map, f, interpolated_corrected_anisotropic_measure_face(u, x)); + + x.clear(); + u.clear(); + } + } + + Vertex_measures expand_interpolated_corrected_measure_vertex_no_radius(Vertex_descriptor v) + { + Vertex_measures vertex_measures; + + for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { + if (f == boost::graph_traits::null_face()) + continue; + + vertex_measures.area_measure += get(mu0_map, f); + + if (is_mean_curvature_selected) + vertex_measures.mean_curvature_measure += get(mu1_map, f); + + if (is_gaussian_curvature_selected) + vertex_measures.gaussian_curvature_measure += get(mu2_map, f); + + if (is_principal_curvatures_and_directions_selected) + { + const std::array face_anisotropic_measure = get(muXY_map, f); + for (std::size_t i = 0; i < 3 * 3; i++) + vertex_measures.anisotropic_measure[i] += face_anisotropic_measure[i]; + } + } + + return vertex_measures; + } + + Vertex_measures expand_interpolated_corrected_measure_vertex(Vertex_descriptor v) + { + std::queue bfs_queue; + std::unordered_set bfs_visited; + + Point_3 vp = get(vpm, v); + Vector_3 c = Vector_3(vp.x(), vp.y(), vp.z()); + + Vertex_measures vertex_measures; + + for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { + if (f != boost::graph_traits::null_face()) + { + bfs_queue.push(f); + bfs_visited.insert(f); + } + } + while (!bfs_queue.empty()) { + Face_descriptor fi = bfs_queue.front(); + bfs_queue.pop(); + + // looping over vertices in face to get point coordinates std::vector x; - std::vector u; - - for (Face_descriptor f : faces(pmesh)) + for (Vertex_descriptor vi : vertices_around_face(halfedge(fi, pmesh), pmesh)) { - for (Vertex_descriptor v : vertices_around_face(halfedge(f, pmesh), pmesh)) - { - Point_3 p = get(vpm, v); - x.push_back(Vector_3(p.x(), p.y(), p.z())); - u.push_back(get(vnm, v)); - } - put(mu0_map, f, interpolated_corrected_area_measure_face(u, x)); - - if (is_mean_curvature_selected) - put(mu1_map, f, interpolated_corrected_mean_curvature_measure_face(u, x)); - - if (is_gaussian_curvature_selected) - put(mu2_map, f, interpolated_corrected_gaussian_curvature_measure_face(u)); - - if (is_principal_curvatures_and_directions_selected) - put(muXY_map, f, interpolated_corrected_anisotropic_measure_face(u, x)); - - x.clear(); - u.clear(); + Point_3 pi = get(vpm, vi); + x.push_back(Vector_3(pi.x(), pi.y(), pi.z())); } - } - Vertex_measures expand_interpolated_corrected_measure_vertex_no_radius(Vertex_descriptor v) - { - Vertex_measures vertex_measures; + const FT f_ratio = face_in_ball_ratio(x, ball_radius, c); - for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { - if (f == boost::graph_traits::null_face()) - continue; - - vertex_measures.area_measure += get(mu0_map, f); + if (f_ratio != 0.0) + { + vertex_measures.area_measure += f_ratio * get(mu0_map, fi); if (is_mean_curvature_selected) - vertex_measures.mean_curvature_measure += get(mu1_map, f); + vertex_measures.mean_curvature_measure += f_ratio * get(mu1_map, fi); if (is_gaussian_curvature_selected) - vertex_measures.gaussian_curvature_measure += get(mu2_map, f); + vertex_measures.gaussian_curvature_measure += f_ratio * get(mu2_map, fi); if (is_principal_curvatures_and_directions_selected) { - const std::array face_anisotropic_measure = get(muXY_map, f); + const std::array face_anisotropic_measure = get(muXY_map, fi); for (std::size_t i = 0; i < 3 * 3; i++) - vertex_measures.anisotropic_measure[i] += face_anisotropic_measure[i]; + vertex_measures.anisotropic_measure[i] += f_ratio * face_anisotropic_measure[i]; + } + + for (Face_descriptor fj : faces_around_face(halfedge(fi, pmesh), pmesh)) + { + if (bfs_visited.find(fj) == bfs_visited.end() && fj != boost::graph_traits::null_face()) + { + bfs_queue.push(fj); + bfs_visited.insert(fj); + } } } - - return vertex_measures; } + return vertex_measures; + } - Vertex_measures expand_interpolated_corrected_measure_vertex(Vertex_descriptor v) + void compute_selected_curvatures() { + interpolated_corrected_selected_measures_all_faces(); + + for (Vertex_descriptor v : vertices(pmesh)) { - std::queue bfs_queue; - std::unordered_set bfs_visited; + Vertex_measures vertex_measures = (ball_radius < 0) ? + expand_interpolated_corrected_measure_vertex_no_radius(v) : + expand_interpolated_corrected_measure_vertex(v); - Point_3 vp = get(vpm, v); - Vector_3 c = Vector_3(vp.x(), vp.y(), vp.z()); - - Vertex_measures vertex_measures; - - for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { - if (f != boost::graph_traits::null_face()) - { - bfs_queue.push(f); - bfs_visited.insert(f); - } + 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) : + put(mean_curvature_map, v, 0); } - while (!bfs_queue.empty()) { - Face_descriptor fi = bfs_queue.front(); - bfs_queue.pop(); - // looping over vertices in face to get point coordinates - std::vector x; - for (Vertex_descriptor vi : vertices_around_face(halfedge(fi, pmesh), pmesh)) - { - Point_3 pi = get(vpm, vi); - x.push_back(Vector_3(pi.x(), pi.y(), pi.z())); - } - - const FT f_ratio = face_in_ball_ratio(x, ball_radius, c); - - if (f_ratio != 0.0) - { - vertex_measures.area_measure += f_ratio * get(mu0_map, fi); - - if (is_mean_curvature_selected) - vertex_measures.mean_curvature_measure += f_ratio * get(mu1_map, fi); - - if (is_gaussian_curvature_selected) - vertex_measures.gaussian_curvature_measure += f_ratio * get(mu2_map, fi); - - if (is_principal_curvatures_and_directions_selected) - { - const std::array face_anisotropic_measure = get(muXY_map, fi); - for (std::size_t i = 0; i < 3 * 3; i++) - vertex_measures.anisotropic_measure[i] += f_ratio * face_anisotropic_measure[i]; - } - - for (Face_descriptor fj : faces_around_face(halfedge(fi, pmesh), pmesh)) - { - if (bfs_visited.find(fj) == bfs_visited.end() && fj != boost::graph_traits::null_face()) - { - bfs_queue.push(fj); - bfs_visited.insert(fj); - } - } - } + if (is_gaussian_curvature_selected) { + vertex_measures.area_measure != 0 ? + put(gaussian_curvature_map, v, vertex_measures.gaussian_curvature_measure / vertex_measures.area_measure) : + put(gaussian_curvature_map, v, 0); } - return vertex_measures; - } - void compute_selected_curvatures() { - interpolated_corrected_selected_measures_all_faces(); - - for (Vertex_descriptor v : vertices(pmesh)) - { - Vertex_measures vertex_measures = (ball_radius < 0) ? - expand_interpolated_corrected_measure_vertex_no_radius(v) : - expand_interpolated_corrected_measure_vertex(v); - - 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) : - put(mean_curvature_map, v, 0); - } - - if (is_gaussian_curvature_selected) { - vertex_measures.area_measure != 0 ? - put(gaussian_curvature_map, v, vertex_measures.gaussian_curvature_measure / vertex_measures.area_measure) : - put(gaussian_curvature_map, v, 0); - } - - if (is_principal_curvatures_and_directions_selected) { - 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, - vertex_measures.area_measure, - v_normal - ); - put(principal_curvatures_and_directions_map, v, principal_curvatures_and_directions); - } + if (is_principal_curvatures_and_directions_selected) { + 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, + vertex_measures.area_measure, + v_normal + ); + put(principal_curvatures_and_directions_map, v, principal_curvatures_and_directions); } } - }; + } +}; } // namespace internal @@ -1030,6 +1266,50 @@ template(pmesh, np); } +template + typename GT::FT interpolated_corrected_mean_curvature_at_vertex(const PolygonMesh& pmesh, + typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values()) +{ + // use interpolated_corrected_curvatures_at_vertex to compute mean curvature + typename GT::FT* mean_curvature = new typename GT::FT(); + interpolated_corrected_curvatures_at_vertex(pmesh, v, np.vertex_mean_curvature(mean_curvature)); + return *mean_curvature; +} + +template + typename GT::FT interpolated_corrected_gaussian_curvature_at_vertex(const PolygonMesh& pmesh, + typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values()) +{ + // use interpolated_corrected_curvatures_at_vertex to compute gaussian curvature + typename GT::FT* gc = new typename GT::FT(); + interpolated_corrected_curvatures_at_vertex(pmesh, v, np.vertex_gaussian_curvature(gc)); + return *gc; +} + +template + Principal_curvatures_and_directions interpolated_corrected_principal_curvatures_and_directions_at_vertex(const PolygonMesh& pmesh, + typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values()) +{ + // use interpolated_corrected_curvatures_at_vertex to compute principal curvatures + Principal_curvatures_and_directions* pcd = new Principal_curvatures_and_directions(); + interpolated_corrected_curvatures_at_vertex(pmesh, v, np.vertex_principal_curvatures_and_directions(pcd)); + return *pcd; +} + +template + void interpolated_corrected_curvatures_at_vertex(const PolygonMesh& pmesh, + typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values()) +{ + internal::interpolated_corrected_curvatures_one_vertex(pmesh, v, np); +} } // namespace Polygon_mesh_processing } // namespace CGAL 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 21b8227e86c..0d546977561 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -90,6 +90,9 @@ CGAL_add_named_parameter(nb_points_per_distance_unit_t, nb_points_per_distance_u CGAL_add_named_parameter(vertex_mean_curvature_map_t, vertex_mean_curvature_map, vertex_mean_curvature_map) CGAL_add_named_parameter(vertex_gaussian_curvature_map_t, vertex_gaussian_curvature_map, vertex_gaussian_curvature_map) CGAL_add_named_parameter(vertex_principal_curvatures_and_directions_map_t, vertex_principal_curvatures_and_directions_map, vertex_principal_curvatures_and_directions_map) +CGAL_add_named_parameter(vertex_mean_curvature_t, vertex_mean_curvature, vertex_mean_curvature) +CGAL_add_named_parameter(vertex_gaussian_curvature_t, vertex_gaussian_curvature, 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(outward_orientation_t, outward_orientation, outward_orientation) CGAL_add_named_parameter(overlap_test_t, overlap_test, do_overlap_test_of_bounded_sides)