diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures.cpp index 505d4328645..4480366bc2b 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures.cpp @@ -33,35 +33,20 @@ int main(int argc, char* argv[]) return EXIT_FAILURE; } - vertexVectorMap_tag vnm_init; - boost::associative_property_map vnm(vnm_init); - - PMP::compute_vertex_normals(g1, vnm); - - FaceMeasureMap_tag mu0_init, mu1_init, mu2_init; + FaceMeasureMap_tag mean_curvature_init, gaussian_curvature_init; boost::associative_property_map - mu0_map(mu0_init), mu1_map(mu1_init), mu2_map(mu2_init); + mean_curvature_map(mean_curvature_init), gaussian_curvature_map(gaussian_curvature_init); - PMP::interpolated_corrected_measure_mesh( + PMP::interpolated_corrected_mean_curvature( g1, - mu0_map, - PMP::MU0_AREA_MEASURE); - - PMP::interpolated_corrected_measure_mesh( + mean_curvature_map + ); + PMP::interpolated_corrected_gaussian_curvature( g1, - mu1_map, - PMP::MU1_MEAN_CURVATURE_MEASURE, - CGAL::parameters::vertex_normal_map(vnm)); - - PMP::interpolated_corrected_measure_mesh( - g1, - mu2_map, - PMP::MU2_GAUSSIAN_CURVATURE_MEASURE, - CGAL::parameters::vertex_normal_map(vnm)); + gaussian_curvature_map + ); for (face_descriptor f: g1.faces()) - std::cout << f.idx() << ": " - << get(mu0_map, f) << ", " - << get(mu1_map, f) << ", " - << get(mu2_map, f) << ", " << "\n"; + std::cout << f.idx() << ": HC = " << get(mean_curvature_map, f) + << ", GC = " << get(gaussian_curvature_map, f) << "\n"; } 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 2409eab42b7..737a2525ebd 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 @@ -11,6 +11,8 @@ #include #include +#include +#include namespace CGAL { @@ -210,8 +212,8 @@ typename GT::FT interpolated_corrected_measure_face(const std::vector::value) compute_vertex_normals(pmesh, vnm, np); for (face_descriptor f : faces(pmesh)) { - halfedge_descriptor h_start = pmesh.halfedge(f); - halfedge_descriptor h_iter = h_start; - std::vector x; std::vector u; - // looping over vertices in face - do { - vertex_descriptor v = source(h_iter, pmesh); + for (vertex_descriptor v : vertices_around_face(halfedge(f, pmesh), pmesh)) + { typename GT::Point_3 p = get(vpm, v); x.push_back(typename GT::Vector_3(p.x(), p.y(), p.z())); u.push_back(get(vnm, v)); - h_iter = next(h_iter, pmesh); - } while (h_iter != h_start); + } put(fmm, f, interpolated_corrected_measure_face(x, u, mu_i)); } } + +// +// +//template +//typename GT::FT triangle_in_ball_ratio_1(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_2(const std::vector& x, + const typename GT::FT r, + const typename GT::Vector_3 c) +{ + 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 + void expand_interpolated_corrected_measure_face(const PolygonMesh& pmesh, + FaceMeasureMap fmm, + const typename boost::graph_traits::face_descriptor f, + const NamedParameters& np = parameters::default_values()) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + typedef typename GetGeomTraits::type GT; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + const typename GetGeomTraits::type::FT + r = choose_parameter(get_parameter(np, internal_np::ball_radius), 0.01); + + if (r < 0.000001) + return; + + typename GetVertexPointMap::const_type + vpm = choose_parameter(get_parameter(np, CGAL::vertex_point), + get_const_property_map(CGAL::vertex_point, pmesh)); + + + std::queue bfs_q; + std::unordered_set bfs_v; + + //get face center c + typename GT::Vector_3 c; + for (vertex_descriptor v : vertices_around_face(halfedge(f, pmesh), pmesh)) + { + typename GT::Point_3 p = get(vpm, v); + c += typename GT::Vector_3(p.x(), p.y(), p.z()); + } + c /= degree(f, pmesh); + + GT::FT corrected_mui = 0; + + bfs_q.push(f); + bfs_v.insert(f); + + while (!bfs_q.empty()) { + face_descriptor fi = bfs_q.front(); + bfs_q.pop(); + + // looping over vertices in face to get point coordinates + std::vector x; + for (vertex_descriptor v : vertices_around_face(halfedge(fi, pmesh), pmesh)) + { + typename GT::Point_3 p = get(vpm, v); + x.push_back(typename GT::Vector_3(p.x(), p.y(), p.z())); + } + + const typename GT::FT f_ratio = face_in_ball_ratio_2(x, r, c); + + if (f_ratio > 0.000001) + { + corrected_mui += f_ratio * get(fmm, fi); + for (face_descriptor fj : faces_around_face(halfedge(fi, pmesh), pmesh)) + { + if (bfs_v.find(fj) == bfs_v.end()) + { + bfs_q.push(fj); + bfs_v.insert(fj); + } + } + } + } + + put(fmm, f, corrected_mui); +} + +//template +// void expand_interpolated_corrected_measure_mesh(const PolygonMesh& pmesh, +// FaceMeasureMap fmm, +// const NamedParameters& np = parameters::default_values()) +//{ +// typedef typename boost::graph_traits::face_descriptor face_descriptor; +// for (face_descriptor f : faces(pmesh)) +// expand_interpolated_corrected_measure_face(pmesh, fmm, f, np); +//} + +template + void interpolated_corrected_mean_curvature(const PolygonMesh& pmesh, + FaceCurvatureMap fcm, + const NamedParameters& np = parameters::default_values()) +{ + typedef typename GetGeomTraits::type GT; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef std::unordered_map FaceMeasureMap_tag; + + FaceMeasureMap_tag mu0_init, mu1_init; + boost::associative_property_map + mu0_map(mu0_init), mu1_map(mu1_init); + + interpolated_corrected_measure_mesh(pmesh, mu0_map, MU0_AREA_MEASURE); + interpolated_corrected_measure_mesh(pmesh, mu1_map, MU1_MEAN_CURVATURE_MEASURE); + + for (face_descriptor f : faces(pmesh)) + { + expand_interpolated_corrected_measure_face(pmesh, mu0_map, f, np); + expand_interpolated_corrected_measure_face(pmesh, mu1_map, f, np); + + GT::FT f_mu0 = get(mu0_map, f); + if (f_mu0 > 0.000001) + put(fcm, f, 0.5 * get(mu1_map, f) / get(mu0_map, f)); + else + put(fcm, f, 0); + } + +} + +template + void interpolated_corrected_gaussian_curvature(const PolygonMesh& pmesh, + FaceCurvatureMap fcm, + const NamedParameters& np = parameters::default_values()) +{ + typedef typename GetGeomTraits::type GT; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef std::unordered_map FaceMeasureMap_tag; + + FaceMeasureMap_tag mu0_init, mu2_init; + boost::associative_property_map + mu0_map(mu0_init), mu2_map(mu2_init); + + interpolated_corrected_measure_mesh(pmesh, mu0_map, MU0_AREA_MEASURE); + interpolated_corrected_measure_mesh(pmesh, mu2_map, MU2_GAUSSIAN_CURVATURE_MEASURE); + + for (face_descriptor f : faces(pmesh)) + { + expand_interpolated_corrected_measure_face(pmesh, mu0_map, f, np); + expand_interpolated_corrected_measure_face(pmesh, mu2_map, f, np); + + GT::FT f_mu0 = get(mu0_map, f); + if(f_mu0 > 0.000001) + put(fcm, f, get(mu2_map, f) / f_mu0); + else + put(fcm, f, 0); + } + +} + + + } } diff --git a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp index 1ea1394649c..88d86966f1d 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp @@ -533,9 +533,8 @@ private Q_SLOTS: dock_widget->propertyBox->addItem("Scaled Jacobian"); dock_widget->propertyBox->addItem("Heat Intensity"); dock_widget->propertyBox->addItem("Heat Intensity (Intrinsic Delaunay)"); - dock_widget->propertyBox->addItem("corrected area density"); - dock_widget->propertyBox->addItem("corrected mean curvature density"); - dock_widget->propertyBox->addItem("corrected gaussian curvature density"); + dock_widget->propertyBox->addItem("Interpolated Corrected Mean Curvature"); + dock_widget->propertyBox->addItem("Interpolated Corrected Gaussian Curvature"); detectSMScalarProperties(sm_item->face_graph()); } @@ -610,15 +609,11 @@ private Q_SLOTS: return; sm_item->setRenderingMode(Gouraud); break; - case 4: // corrected area density - displayInterpolatedCurvatureMeasure(sm_item, PMP::MU0_AREA_MEASURE); - sm_item->setRenderingMode(Flat); - break; - case 5: // corrected mean curvature density + case 4: // Interpolated Corrected Mean Curvature displayInterpolatedCurvatureMeasure(sm_item, PMP::MU1_MEAN_CURVATURE_MEASURE); sm_item->setRenderingMode(Flat); break; - case 6: // corrected gaussian curvature density + case 5: // Interpolated Corrected Gaussian Curvature displayInterpolatedCurvatureMeasure(sm_item, PMP::MU2_GAUSSIAN_CURVATURE_MEASURE); sm_item->setRenderingMode(Flat); break; @@ -657,15 +652,11 @@ private Q_SLOTS: if(does_exist) sm_item->face_graph()->remove_property_map(pmap); std::tie(pmap, does_exist) = - sm_item->face_graph()->property_map("f:corrected_area_density"); + sm_item->face_graph()->property_map("f:interpolated_corrected_mean_curvature"); if (does_exist) sm_item->face_graph()->remove_property_map(pmap); std::tie(pmap, does_exist) = - sm_item->face_graph()->property_map("f:corrected_mean_curvature_density"); - if (does_exist) - sm_item->face_graph()->remove_property_map(pmap); - std::tie(pmap, does_exist) = - sm_item->face_graph()->property_map("f:corrected_gaussian_curvature_density"); + sm_item->face_graph()->property_map("f:interpolated_corrected_gaussian_curvature"); if (does_exist) sm_item->face_graph()->remove_property_map(pmap); }); @@ -707,16 +698,12 @@ private Q_SLOTS: dock_widget->zoomToMaxButton->setEnabled(jacobian_max.count(sm_item)>0); break; case 4: - dock_widget->zoomToMinButton->setEnabled(mu0_min.count(sm_item) > 0); - dock_widget->zoomToMaxButton->setEnabled(mu0_max.count(sm_item) > 0); + dock_widget->zoomToMinButton->setEnabled(mean_curvature_min.count(sm_item) > 0); + dock_widget->zoomToMaxButton->setEnabled(mean_curvature_max.count(sm_item) > 0); break; case 5: - dock_widget->zoomToMinButton->setEnabled(mu1_min.count(sm_item) > 0); - dock_widget->zoomToMaxButton->setEnabled(mu1_max.count(sm_item) > 0); - break; - case 6: - dock_widget->zoomToMinButton->setEnabled(mu2_min.count(sm_item) > 0); - dock_widget->zoomToMaxButton->setEnabled(mu2_max.count(sm_item) > 0); + dock_widget->zoomToMinButton->setEnabled(gaussian_curvature_min.count(sm_item) > 0); + dock_widget->zoomToMaxButton->setEnabled(gaussian_curvature_max.count(sm_item) > 0); break; default: break; @@ -744,23 +731,17 @@ private Q_SLOTS: { smesh.remove_property_map(angles); } - SMesh::Property_map mu0; - std::tie(mu0, found) = smesh.property_map("f:corrected_area_density"); + SMesh::Property_map mean_curvature; + std::tie(mean_curvature, found) = smesh.property_map("f:interpolated_corrected_mean_curvature"); if (found) { - smesh.remove_property_map(mu0); + smesh.remove_property_map(mean_curvature); } - SMesh::Property_map mu1; - std::tie(mu1, found) = smesh.property_map("f:corrected_mean_curvature_density"); + SMesh::Property_map gaussian_curvature; + std::tie(gaussian_curvature, found) = smesh.property_map("f:interpolated_corrected_gaussian_curvature"); if (found) { - smesh.remove_property_map(mu0); - } - SMesh::Property_map mu2; - std::tie(mu2, found) = smesh.property_map("f:corrected_gaussian_curvature_density"); - if (found) - { - smesh.remove_property_map(mu0); + smesh.remove_property_map(gaussian_curvature); } } @@ -804,20 +785,20 @@ private Q_SLOTS: void displayInterpolatedCurvatureMeasure(Scene_surface_mesh_item* item, PMP::Curvature_measure_index mu_index) { - std::vector tied_map = { - "f:corrected_area_density", - "f:corrected_mean_curvature_density", - "f:corrected_gaussian_curvature_density" - }; + std::string tied_string = (mu_index == PMP::MU1_MEAN_CURVATURE_MEASURE)? + "f:interpolated_corrected_mean_curvature": "f:interpolated_corrected_gaussian_curvature"; SMesh& smesh = *item->face_graph(); //compute once and store the value per face bool non_init; SMesh::Property_map mu_i_map; std::tie(mu_i_map, non_init) = - smesh.add_property_map(tied_map[mu_index], 0); + smesh.add_property_map(tied_string, 0); if (non_init) { - PMP::interpolated_corrected_measure_mesh(smesh, mu_i_map, mu_index); + if (mu_index == PMP::MU1_MEAN_CURVATURE_MEASURE) + PMP::interpolated_corrected_mean_curvature(smesh, mu_i_map); + else + PMP::interpolated_corrected_gaussian_curvature(smesh, mu_i_map); double res_min = ARBITRARY_DBL_MAX, res_max = -ARBITRARY_DBL_MAX; @@ -837,24 +818,20 @@ private Q_SLOTS: } switch (mu_index) { - case PMP::MU0_AREA_MEASURE: - mu0_max[item] = std::make_pair(res_max, max_index); - mu0_min[item] = std::make_pair(res_min, min_index); - break; case PMP::MU1_MEAN_CURVATURE_MEASURE: - mu1_max[item] = std::make_pair(res_max, max_index); - mu1_min[item] = std::make_pair(res_min, min_index); + mean_curvature_max[item] = std::make_pair(res_max, max_index); + mean_curvature_min[item] = std::make_pair(res_min, min_index); break; case PMP::MU2_GAUSSIAN_CURVATURE_MEASURE: - mu2_max[item] = std::make_pair(res_max, max_index); - mu2_min[item] = std::make_pair(res_min, min_index); + gaussian_curvature_max[item] = std::make_pair(res_max, max_index); + gaussian_curvature_min[item] = std::make_pair(res_min, min_index); break; } connect(item, &Scene_surface_mesh_item::itemChanged, this, &DisplayPropertyPlugin::resetProperty); } - treat_sm_property(tied_map[mu_index], item->face_graph()); + treat_sm_property(tied_string, item->face_graph()); } @@ -1159,7 +1136,6 @@ private Q_SLOTS: case 1: case 4: case 5: - case 6: dock_widget->groupBox-> setEnabled(true); dock_widget->groupBox_3->setEnabled(true); @@ -1223,7 +1199,7 @@ private Q_SLOTS: case 4: { ::zoomToId(*item->face_graph(), - QString("f%1").arg(mu0_min[item].second), + QString("f%1").arg(mean_curvature_min[item].second), getActiveViewer(), dummy_fd, dummy_p); @@ -1232,20 +1208,12 @@ private Q_SLOTS: case 5: { ::zoomToId(*item->face_graph(), - QString("f%1").arg(mu1_min[item].second), + QString("f%1").arg(gaussian_curvature_min[item].second), getActiveViewer(), dummy_fd, dummy_p); } break; - case 6: - { - ::zoomToId(*item->face_graph(), - QString("f%1").arg(mu2_min[item].second), - getActiveViewer(), - dummy_fd, - dummy_p); - } break; break; default: @@ -1284,7 +1252,7 @@ private Q_SLOTS: case 4: { ::zoomToId(*item->face_graph(), - QString("f%1").arg(mu0_max[item].second), + QString("f%1").arg(mean_curvature_max[item].second), getActiveViewer(), dummy_fd, dummy_p); @@ -1293,21 +1261,12 @@ private Q_SLOTS: case 5: { ::zoomToId(*item->face_graph(), - QString("f%1").arg(mu1_max[item].second), + QString("f%1").arg(gaussian_curvature_max[item].second), getActiveViewer(), dummy_fd, dummy_p); } break; - case 6: - { - ::zoomToId(*item->face_graph(), - QString("f%1").arg(mu2_max[item].second), - getActiveViewer(), - dummy_fd, - dummy_p); - } - break; default: break; } @@ -1598,14 +1557,11 @@ private: std::unordered_map > angles_min; std::unordered_map > angles_max; - std::unordered_map > mu0_min; - std::unordered_map > mu0_max; + std::unordered_map > mean_curvature_min; + std::unordered_map > mean_curvature_max; - std::unordered_map > mu1_min; - std::unordered_map > mu1_max; - - std::unordered_map > mu2_min; - std::unordered_map > mu2_max; + std::unordered_map > gaussian_curvature_min; + std::unordered_map > gaussian_curvature_max; std::unordered_map is_source; 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 c224b498251..dd4785ab7e8 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -87,6 +87,7 @@ CGAL_add_named_parameter(number_of_points_per_edge_t, number_of_points_per_edge, CGAL_add_named_parameter(number_of_points_on_edges_t, number_of_points_on_edges, number_of_points_on_edges) CGAL_add_named_parameter(nb_points_per_area_unit_t, nb_points_per_area_unit, number_of_points_per_area_unit) CGAL_add_named_parameter(nb_points_per_distance_unit_t, nb_points_per_distance_unit, number_of_points_per_distance_unit) +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) CGAL_add_named_parameter(preserve_genus_t, preserve_genus, preserve_genus)