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

debugging curvature flow

This commit is contained in:
konstantinos katrioplas 2017-08-12 22:58:00 +03:00
parent e6d3c9224c
commit 42e7d12df6
3 changed files with 139 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

156
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h
View File

@ -20,6 +20,7 @@
#include <utility>
#include <math.h>
#include <CGAL/Monge_via_jet_fitting.h>
namespace CGAL {
@ -41,6 +42,7 @@ class Curvature_flow
typedef typename GeomTraits::Point_3 Point;
typedef typename GeomTraits::Vector_3 Vector;
typedef typename GeomTraits::Triangle_3 Triangle;
typedef typename GeomTraits::FT FT;
typedef std::vector<Triangle> Triangle_list;
typedef CGAL::AABB_triangle_primitive<GeomTraits, typename Triangle_list::iterator> AABB_Primitive;
@ -48,15 +50,32 @@ class Curvature_flow
typedef CGAL::AABB_tree<AABB_Traits> Tree;
// <one halfedge around v, pair of incident halfedges to this halfedge around v>
typedef std::pair<halfedge_descriptor, halfedge_descriptor> He_pair;
typedef std::map<halfedge_descriptor, He_pair> Edges_around_map;
typedef std::pair<halfedge_descriptor, halfedge_descriptor> he_pair;
typedef std::map<halfedge_descriptor, he_pair> Edges_around_map;
typedef typename CGAL::Monge_via_jet_fitting<GeomTraits> Monge_via_jet_fitting;
typedef typename Monge_via_jet_fitting::Monge_form Monge_form;
public:
Curvature_flow(PolygonMesh& pmesh, VertexPointMap& vpmap, VertexConstraintMap& vcmap, EdgeConstraintMap& ecmap) :
mesh_(pmesh), vpmap_(vpmap), vcmap_(vcmap), ecmap_(ecmap), cot_calculator_(pmesh, vpmap)
{}
{
/*
std::size_t num_points = vertices(mesh_).size();
std::size_t min_num_of_points = 6; // (d+1)(d+2)/2, for d=2
if(num_points < min_num_of_points)
{
CGAL_error_msg("Find curvature: Not enough points in the mesh.");
}
*/
}
template<typename FaceRange>
void init_remeshing(const FaceRange& face_range)
@ -72,6 +91,8 @@ public:
tree_ptr_ = new Tree(input_triangles_.begin(), input_triangles_.end());
tree_ptr_->accelerate_distance_queries();
//mean_k_ = compute_mean_curvature();
}
std::size_t remove_degenerate_faces()
@ -88,6 +109,8 @@ public:
return nb_removed_faces;
}
void curvature_smoothing()
{
std::map<vertex_descriptor, Point> barycenters;
@ -104,45 +127,54 @@ public:
.geom_traits(traits_));
n_map[v] = vn;
// find incident halfedges
Edges_around_map he_map;
typename Edges_around_map::iterator it;
BOOST_FOREACH(halfedge_descriptor hi, halfedges_around_source(v, mesh_))
he_map[hi] = He_pair( next(hi, mesh_), prev(opposite(hi, mesh_), mesh_) );
he_map[hi] = he_pair( next(hi, mesh_), prev(opposite(hi, mesh_), mesh_) );
Vector weighted_move = CGAL::NULL_VECTOR;
// calculate movement
Vector curvature_normal = CGAL::NULL_VECTOR;
double sum_cot_weights = 0;
for(it = he_map.begin(); it!= he_map.end(); ++it)
{
halfedge_descriptor hi = it->first; //main_he
halfedge_descriptor hi = it->first;
he_pair incd_edges = it->second;
// weight
double weight = cot_angles(hi, it->second); // incd_edges
double weight = cot_angles(hi, incd_edges);
sum_cot_weights += weight;
// displacement vector
Point Xi = get(vpmap_, source(hi, mesh_));
Point Xj = get(vpmap_, target(hi, mesh_));
Vector vec(Xj, Xi);
//Vector vec(Xj, Xi);
Vector vec(Xi, Xj);
// add weight
vec *= weight;
// sum vecs
weighted_move += vec;
curvature_normal += vec;
}
// divide with total weight - if there is actually weight
if(sum_cot_weights != 0)
weighted_move /= sum_cot_weights;
curvature_normal /= sum_cot_weights;
Point weighted_barycenter = get(vpmap_, v) - weighted_move;
//Point weighted_barycenter = get(vpmap_, v) - curvature_normal;
Point weighted_barycenter = get(vpmap_, v) + curvature_normal;
barycenters[v] = weighted_barycenter;
} // not on border
} // all vertices
typedef typename std::map<vertex_descriptor, Point>::value_type VP;
/*
// compute locations on tangent plane
typedef typename std::map<vertex_descriptor, Point>::value_type VP;
std::map<vertex_descriptor, Point> new_locations;
@ -155,10 +187,13 @@ public:
new_locations[vp.first] = q + ( n * Vector(q, p) ) * n ;
}
// update location
BOOST_FOREACH(const VP& vp, new_locations)
{
*/
// update location
//BOOST_FOREACH(const VP& vp, new_locations) // uncomment for tangent plane
BOOST_FOREACH(const VP& vp, barycenters)
{
#define CGAL_PMP_SMOOTHING_DEBUG
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "from: "<< get(vpmap_, vp.first);
#endif
@ -194,46 +229,55 @@ public:
//take one halfedge whose target is v
BOOST_FOREACH(halfedge_descriptor ht, halfedges_around_target(v, mesh_))
{
// is it ok if it is degenerate?
// is it ok if a face is degenerate?
A = area(faces_around_target(ht, mesh_), mesh_);
continue;
}
// find incident halfedges
Edges_around_map he_map;
typename Edges_around_map::iterator it;
BOOST_FOREACH(halfedge_descriptor hi, halfedges_around_source(v, mesh_))
he_map[hi] = He_pair( next(hi, mesh_), prev(opposite(hi, mesh_), mesh_) );
he_map[hi] = he_pair( next(hi, mesh_), prev(opposite(hi, mesh_), mesh_) );
// calculate movement
Vector weighted_move = CGAL::NULL_VECTOR;
for(it = he_map.begin(); it!= he_map.end(); ++it)
{
halfedge_descriptor hi = it->first; //main_he
halfedge_descriptor hi = it->first;
he_pair incd_edges = it->second;
// weight
double weight = cot_angles(hi, it->second); // incd_edges
double weight = cot_angles(hi, incd_edges);
// displacement vector
Point xi = get(vpmap_, source(hi, mesh_));
Point xj = get(vpmap_, target(hi, mesh_));
Vector vec(xj, xi);
Vector vec(xi, xj);
// add weight
vec *= weight; // comment out for just laplacian
// add weight.
vec *= weight;
// sum vecs
weighted_move += vec;
std::cout<<"hi";
}
Vector curvature_normal = CGAL::NULL_VECTOR;
if (A != 0)
{
curvature_normal = weighted_move / (4 * A);
//curvature_normal = weighted_move / he_map.size(); // just laplacian
}
Point old_point = get(vpmap_, v);
Point new_point = get(vpmap_, v) + curvature_normal;
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "new location before projection: "<< new_point << std::endl;
#endif
barycenters[v] = new_point;
} // not on border
@ -241,7 +285,7 @@ public:
} // all vertices
/*
// calculate location on tangential plane
std::map<vertex_descriptor, Point> new_locations;
BOOST_FOREACH(const auto& vp, barycenters)
{
@ -251,12 +295,11 @@ public:
new_locations[vp.first] = q + ( n * Vector(q, p) ) * n ;
}
*/
// update location
BOOST_FOREACH(const auto& vp, barycenters)
//BOOST_FOREACH(const auto& vp, new_locations)
//BOOST_FOREACH(const auto& vp, barycenters)
BOOST_FOREACH(const auto& vp, new_locations)
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
@ -312,7 +355,7 @@ private:
return sqlength(halfedge(e, mesh_));
}
double cot_angles(const halfedge_descriptor& main_he, const He_pair& incd_edges)
double cot_angles(const halfedge_descriptor& main_he, const he_pair& incd_edges)
{
vertex_descriptor vs = source(main_he, mesh_);
vertex_descriptor vt = target(main_he, mesh_);
@ -325,10 +368,11 @@ private:
Point p2 = get(vpmap_, v2);
Point pt = get(vpmap_, vt);
Point ps = get(vpmap_, vs);
// avoid degenerate cases
Vector edge1(pt, p1);
Vector edge2(pt, p2);
Vector vec_main_he(pt, ps);
double tolerance = 1e-3;
if ( edge1.squared_length() < tolerance ||
edge2.squared_length() < tolerance ||
@ -348,6 +392,61 @@ private:
return a1 + a2;
}
/*
double compute_mean_curvature()
{
// gather points around v
std::vector<Point> incident_points = gather_all_points();
Monge_form monge_form;
Monge_via_jet_fitting monge_fit;
std::size_t d_fit = 2; // d_fit >= d_monge
std::size_t d_monge = 2; // need 2 principal coeeficients
std::size_t Nd = (d_fit + 1)*(d_fit + 1) / 2.0;
CGAL_assertion(incident_points.size() >= Nd);
monge_form = monge_fit(incident_points.begin(), incident_points.end(), d_fit, d_monge);
const double k1 = monge_form.principal_curvatures(0);
const double k2 = monge_form.principal_curvatures(1);
return (k1 + k2) / 2.0;
}
std::vector<Point> points_around_vertex(vertex_descriptor v)
{
std::vector<Point> incident_vertices;
for(halfedge_descriptor h : halfedges_around_target(v, mesh_))
{
vertex_descriptor vs = source(h, mesh_);
incident_vertices.push_back(get(vpmap_, vs));
}
// temp assertion
std::vector<Point> incident_vertices2;
for(vertex_descriptor vi : vertices_around_target(v, mesh_))
{
incident_vertices2.push_back(get(vpmap_, vi));
}
CGAL_assertion(incident_vertices.size() == incident_vertices2.size());
return incident_vertices;
}
std::vector<Point> gather_all_points()
{
//TODO SEE IF THERE IS SOMTEHTING alREADY FOR POINTS
std::vector<Point> points;
for(vertex_descriptor v : vertices(mesh_))
{
points.push_back(get(vpmap_, v)); // todo: preallocate it and fill it by pointing
}
return points;
}
*/
bool is_constrained(const edge_descriptor& e)
{
return get(ecmap_, e);
@ -381,6 +480,7 @@ private:
Tree* tree_ptr_;
GeomTraits traits_;
double min_sq_edge_len_;
//double mean_k_;
// from Weights.h
CGAL::internal::Cotangent_value_Meyer_secure<PolygonMesh, VertexPointMap> cot_calculator_;
//CGAL::internal::Cotangent_value_clamped_2<PolygonMesh, VertexPointMap> cot_calculator_;

12
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_impl.h
View File

@ -49,8 +49,8 @@ class Compatible_remesher
typedef CGAL::AABB_tree<AABB_Traits> Tree;
// <one halfedge around v, pair of incident halfedges to this halfedge around v>
typedef std::pair<halfedge_descriptor, halfedge_descriptor> He_pair;
typedef std::map<halfedge_descriptor, He_pair> Edges_around_map;
typedef std::pair<halfedge_descriptor, halfedge_descriptor> he_pair;
typedef std::map<halfedge_descriptor, he_pair> Edges_around_map;
public:
@ -112,7 +112,7 @@ public:
Edges_around_map he_map;
typename Edges_around_map::iterator it;
BOOST_FOREACH(halfedge_descriptor hi, halfedges_around_source(v, mesh_))
he_map[hi] = He_pair( next(hi, mesh_), prev(opposite(hi, mesh_), mesh_) );
he_map[hi] = he_pair( next(hi, mesh_), prev(opposite(hi, mesh_), mesh_) );
// calculate movement
Vector move = CGAL::NULL_VECTOR;
@ -120,7 +120,7 @@ public:
for(it = he_map.begin(); it != he_map.end(); ++it)
{
halfedge_descriptor main_he = it->first;
He_pair incident_pair = it->second;
he_pair incident_pair = it->second;
Vector rotated_edge = rotate_edge(main_he, incident_pair);
@ -434,7 +434,7 @@ private:
return move_flag;
}
Vector rotate_edge(const halfedge_descriptor& main_he, const He_pair& incd_edges)
Vector rotate_edge(const halfedge_descriptor& main_he, const he_pair& incd_edges)
{
// get common vertex around which the edge is rotated
Point pt = get(vpmap_, target(main_he, mesh_)); // CORRECT SYNATX vt
@ -502,7 +502,7 @@ private:
return bisector;
}
double get_angle(const He_pair& incd_edges, const Vector& vn)
double get_angle(const he_pair& incd_edges, const Vector& vn)
{
Point s = get(vpmap_, source(incd_edges.first, mesh_));
Point p1 = get(vpmap_, target(incd_edges.first, mesh_));

10
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smoothing.h
View File

@ -2,8 +2,8 @@
#define CGAL_POLYGON_MESH_PROCESSING_SMOOTHING_H
#include <CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/smoothing_impl.h>
#include <CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/curvature_flow_impl.h>
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_impl.h>
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h>
#include <CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
@ -326,11 +326,11 @@ void curvature_flow(PolygonMesh& pmesh, const FaceRange& faces, const NamedParam
std::cout << " * Iteration " << (i + 1) << " *" << std::endl;
#endif
curvature_remesher.curvature_smoothing();
//curvature_remesher.good_curvature_smoothing();
curvature_remesher.curvature_smoothing(); // normalized version sec 5.5
//curvature_remesher.good_curvature_smoothing(); // formula 14
//for now
//curvature_remesher.project_to_surface();
curvature_remesher.project_to_surface();
}