songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Mean and Gaussian Curvatures + Visualizer (Still wip)

This commit is contained in:
hoskillua 2022-07-26 19:32:17 +02:00
parent 66a2624641
commit 41be3688ae
4 changed files with 259 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures.cpp
View File

@ -33,35 +33,20 @@ int main(int argc, char* argv[])
return EXIT_FAILURE;
}
vertexVectorMap_tag vnm_init;
boost::associative_property_map <vertexVectorMap_tag> 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<FaceMeasureMap_tag>
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";
}

224
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Curvatures/interpolated_corrected_curvature_measures.h
View File

@ -11,6 +11,8 @@
#include <CGAL/property_map.h>
#include <numeric>
#include <queue>
#include <unordered_set>
namespace CGAL {
@ -210,8 +212,8 @@ typename GT::FT interpolated_corrected_measure_face(const std::vector<typename G
else {
typename GT::FT mu0 = 0;
// getting barycenter of points
typename GT::Vector_3 xm =
// getting center of points
typename GT::Vector_3 xm =
std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0));
xm /= n;
@ -300,31 +302,229 @@ template<typename PolygonMesh, typename FaceMeasureMap,
vpm = choose_parameter(get_parameter(np, CGAL::vertex_point),
get_const_property_map(CGAL::vertex_point, pmesh));
// TODO - handle if vnm is not provided
VNM vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map), get(Vector_map_tag(), pmesh));
VNM vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map),
get(Vector_map_tag(), pmesh));
if (is_default_parameter<NamedParameters, internal_np::vertex_normal_map_t>::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<typename GT::Vector_3> x;
std::vector<typename GT::Vector_3> 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<GT>(x, u, mu_i));
}
}
//
//
//template<typename GT>
//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>
typename GT::FT face_in_ball_ratio_2(const std::vector<typename GT::Vector_3>& 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<typename PolygonMesh, typename FaceMeasureMap,
typename NamedParameters = parameters::Default_named_parameters>
void expand_interpolated_corrected_measure_face(const PolygonMesh& pmesh,
FaceMeasureMap fmm,
const typename boost::graph_traits<PolygonMesh>::face_descriptor f,
const NamedParameters& np = parameters::default_values())
{
using parameters::choose_parameter;
using parameters::get_parameter;
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
const typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT
r = choose_parameter(get_parameter(np, internal_np::ball_radius), 0.01);
if (r < 0.000001)
return;
typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type
vpm = choose_parameter(get_parameter(np, CGAL::vertex_point),
get_const_property_map(CGAL::vertex_point, pmesh));
std::queue<face_descriptor> bfs_q;
std::unordered_set<face_descriptor> 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<typename GT::Vector_3> 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<GT>(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<typename PolygonMesh, typename FaceMeasureMap,
// typename NamedParameters = parameters::Default_named_parameters>
// void expand_interpolated_corrected_measure_mesh(const PolygonMesh& pmesh,
// FaceMeasureMap fmm,
// const NamedParameters& np = parameters::default_values())
//{
// typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
// for (face_descriptor f : faces(pmesh))
// expand_interpolated_corrected_measure_face(pmesh, fmm, f, np);
//}
template<typename PolygonMesh, typename FaceCurvatureMap,
typename NamedParameters = parameters::Default_named_parameters>
void interpolated_corrected_mean_curvature(const PolygonMesh& pmesh,
FaceCurvatureMap fcm,
const NamedParameters& np = parameters::default_values())
{
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef std::unordered_map<face_descriptor, GT::FT> FaceMeasureMap_tag;
FaceMeasureMap_tag mu0_init, mu1_init;
boost::associative_property_map<FaceMeasureMap_tag>
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<typename PolygonMesh, typename FaceCurvatureMap,
typename NamedParameters = parameters::Default_named_parameters>
void interpolated_corrected_gaussian_curvature(const PolygonMesh& pmesh,
FaceCurvatureMap fcm,
const NamedParameters& np = parameters::default_values())
{
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef std::unordered_map<face_descriptor, GT::FT> FaceMeasureMap_tag;
FaceMeasureMap_tag mu0_init, mu2_init;
boost::associative_property_map<FaceMeasureMap_tag>
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);
}
}
}
}

116
Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp
View File

@ -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<face_descriptor, double>("f:corrected_area_density");
sm_item->face_graph()->property_map<face_descriptor, double>("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<face_descriptor, double>("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<face_descriptor, double>("f:corrected_gaussian_curvature_density");
sm_item->face_graph()->property_map<face_descriptor, double>("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<face_descriptor, double> mu0;
std::tie(mu0, found) = smesh.property_map<face_descriptor, double>("f:corrected_area_density");
SMesh::Property_map<face_descriptor, double> mean_curvature;
std::tie(mean_curvature, found) = smesh.property_map<face_descriptor, double>("f:interpolated_corrected_mean_curvature");
if (found)
{
smesh.remove_property_map(mu0);
smesh.remove_property_map(mean_curvature);
}
SMesh::Property_map<face_descriptor, double> mu1;
std::tie(mu1, found) = smesh.property_map<face_descriptor, double>("f:corrected_mean_curvature_density");
SMesh::Property_map<face_descriptor, double> gaussian_curvature;
std::tie(gaussian_curvature, found) = smesh.property_map<face_descriptor, double>("f:interpolated_corrected_gaussian_curvature");
if (found)
{
smesh.remove_property_map(mu0);
}
SMesh::Property_map<face_descriptor, double> mu2;
std::tie(mu2, found) = smesh.property_map<face_descriptor, double>("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<std::string> 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<face_descriptor, double> mu_i_map;
std::tie(mu_i_map, non_init) =
smesh.add_property_map<face_descriptor, double>(tied_map[mu_index], 0);
smesh.add_property_map<face_descriptor, double>(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<face_descriptor>(tied_map[mu_index], item->face_graph());
treat_sm_property<face_descriptor>(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<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > angles_min;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > angles_max;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > mu0_min;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > mu0_max;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > mean_curvature_min;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > mean_curvature_max;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > mu1_min;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > mu1_max;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > mu2_min;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > mu2_max;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > gaussian_curvature_min;
std::unordered_map<Scene_surface_mesh_item*, std::pair<double, SMesh::Face_index> > gaussian_curvature_max;
std::unordered_map<Scene_surface_mesh_item*, Vertex_source_map> is_source;

1
STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
View File

@ -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)