17 #include <Eigen/Dense>
19 #include <boost/graph/graph_traits.hpp>
20 #include <boost/graph/properties.hpp>
22 #include <CGAL/boost/graph/iterator.h>
23 #include <CGAL/boost/graph/helpers.h>
38 template<
typename HalfedgeGraph,
51 template<
typename HalfedgeGraph,
59 template<
typename HalfedgeGraph,
65 typename boost::graph_traits< HalfedgeGraph >::face_descriptor fd,
66 const GeometryTraits >)
68 typedef typename boost::property_traits< PointMap >::value_type
Point3d;
71 CGAL::Halfedge_around_face_circulator< HalfedgeGraph > p_halfedge(
77 return FEVV::Operators::Geometry::triangle_area< GeometryTraits >(
84 template<
typename HalfedgeGraph,
86 typename FaceNormalMap,
88 typename GraphTraits = boost::graph_traits< HalfedgeGraph >,
93 const FaceNormalMap &f_nm,
95 double pp_matrix_sum[3][3],
96 const GeometryTraits >)
100 typedef typename GraphTraits::face_descriptor face_descriptor;
105 CGAL::Halfedge_around_target_circulator< HalfedgeGraph > p_halfedge(vd, g),
112 #if 1 // WAITING for AIF stuff - TODO : REMOVE THIS BLOC LATER...
115 Vector edge(p1[0] - p2[0], p1[1] - p2[1], p1[2] - p2[2]);
117 double len_edge = gt.length(
edge);
122 face_descriptor p_facet1 =
face(*p_halfedge, g);
123 face_descriptor p_facet2 =
face(
opposite(*p_halfedge, g), g);
124 assert(p_facet1 != p_facet2);
125 if(p_facet1 == GraphTraits::null_face() ||
126 p_facet2 == GraphTraits::null_face())
131 Vector normal1 = f_nm[p_facet1];
132 Vector normal2 = f_nm[p_facet2];
135 Vector cp(normal1[1] * normal2[2] - normal1[2] * normal2[1],
136 normal1[2] * normal2[0] - normal1[0] * normal2[2],
137 normal1[0] * normal2[1] - normal1[1] * normal2[0]);
138 double ps_cp_xedge = cp[0] *
edge[0] + cp[1] *
edge[1] + cp[2] *
edge[2];
148 double pp_matrix[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
150 p_vector_edge, pp_matrix, beta / len_edge);
152 }
while(++p_halfedge != done);
154 for(
int i = 0; i < 3; i++)
155 for(
int j = 0; j < 3; j++)
156 pp_matrix_sum[i][j] /= area;
162 template<
typename HalfedgeGraph,
164 typename FaceNormalMap,
166 typename GraphTraits = boost::graph_traits< HalfedgeGraph >,
171 const FaceNormalMap &f_nm,
173 double pp_matrix_sum[3][3],
175 const GeometryTraits >)
179 typedef typename GraphTraits::face_descriptor face_descriptor;
181 std::set< VertexDescriptor >
vertices;
183 std::stack< VertexDescriptor > s;
189 VertexDescriptor v = s.top();
192 CGAL::Halfedge_around_target_circulator< HalfedgeGraph > h(v, g), done(h);
199 #if 1 // WAITING for AIF stuff - TODO : REMOVE THIS BLOC LATER...
200 Vector vec = gt.sub_p(p2, p1);
201 Vector pm_o = gt.sub_p(p, o);
203 Vector vec(p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]);
204 Vector PmO(P[0] - O[0], P[1] - O[1], P[2] - O[2]);
207 vec[0] * pm_o[0] + vec[1] * pm_o[1] + vec[2] * pm_o[2];
217 if(!CGAL::is_border_edge(*h, g))
219 double len_edge = std::sqrt(
220 vec[0] * vec[0] + vec[1] * vec[1] +
223 face_descriptor p_facet1 =
face(*h, g);
225 assert(p_facet1 != p_facet2);
226 if(p_facet1 == GraphTraits::null_face() ||
227 p_facet2 == GraphTraits::null_face())
230 Vector normal1 = f_nm[p_facet1];
231 Vector normal2 = f_nm[p_facet2];
234 Vector cp(normal1[1] * normal2[2] - normal1[2] * normal2[1],
235 normal1[2] * normal2[0] - normal1[0] * normal2[2],
236 normal1[0] * normal2[1] - normal1[1] * normal2[0]);
237 double ps_cpx_v = cp[0] * vec[0] + cp[1] * vec[1] + cp[2] * vec[2];
245 double p_vector_edge[3] = {vec[0], vec[1], vec[2]};
246 double pp_matrix[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
248 p_vector_edge, pp_matrix, beta / len_edge);
254 VertexDescriptor w =
source(
258 assert(w != GraphTraits::null_vertex());
267 }
while(++h != done);
272 double area = M_PI * radius * radius;
274 for(
int i = 0; i < 3; i++)
275 for(
int j = 0; j < 3; j++)
276 pp_matrix_sum[i][j] /= area;
306 template<
typename HalfedgeGraph,
307 typename VertexCurvatureMap,
309 typename FaceNormalMap,
313 VertexCurvatureMap &v_cm,
315 const FaceNormalMap &f_nm,
318 double &min_nrm_min_curvature,
319 double &max_nrm_min_curvature,
320 double &min_nrm_max_curvature,
321 double &max_nrm_max_curvature,
322 const GeometryTraits >)
324 typedef boost::graph_traits< HalfedgeGraph > GraphTraits;
325 typedef typename GraphTraits::vertex_iterator vertex_iterator;
329 min_nrm_min_curvature = 100000;
330 max_nrm_min_curvature = -100000;
332 min_nrm_max_curvature = 100000;
333 max_nrm_max_curvature = -100000;
337 std::cout <<
"Asking to calculate curvature" << std::endl;
338 #ifdef DBG_calculate_curvature
339 std::cout <<
"Real radius = " << radius << std::endl;
342 unsigned int skipped = 0;
344 const auto iterator_pair =
346 vertex_iterator vi = iterator_pair.first;
347 vertex_iterator vi_end = iterator_pair.second;
348 for(; vi != vi_end; ++vi)
351 if(
halfedge(*vi, g) == GraphTraits::null_halfedge())
353 v_cm[*vi].VKmaxCurv = gt.NULL_VECTOR;
354 v_cm[*vi].VKminCurv = gt.NULL_VECTOR;
355 v_cm[*vi].KmaxCurv = 0;
356 v_cm[*vi].KminCurv = 0;
361 bool no_val_pro =
false;
363 double pp_matrix_sum[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
398 cov_mat(0, 0) = pp_matrix_sum[0][0];
399 cov_mat(0, 1) = pp_matrix_sum[0][1];
400 cov_mat(0, 2) = pp_matrix_sum[0][2];
401 cov_mat(1, 0) = pp_matrix_sum[1][0];
402 cov_mat(1, 1) = pp_matrix_sum[1][1];
403 cov_mat(1, 2) = pp_matrix_sum[1][2];
404 cov_mat(2, 0) = pp_matrix_sum[2][0];
405 cov_mat(2, 1) = pp_matrix_sum[2][1];
406 cov_mat(2, 2) = pp_matrix_sum[2][2];
408 #ifdef DBG_calculate_curvature
409 std::cout <<
"\n\nCovMat\n\n" << CovMat << std::endl;
412 if(pp_matrix_sum[0][1] == 0 && pp_matrix_sum[0][2] == 0 &&
413 pp_matrix_sum[1][0] == 0 && pp_matrix_sum[1][2] == 0 &&
414 pp_matrix_sum[2][1] == 0 && pp_matrix_sum[2][0] == 0)
416 for(
int i = 0; i < 3; i++)
417 valpro(i) = cov_mat(i, i);
418 vect_pro = EigenMatrix::Identity();
428 #ifdef DBG_calculate_curvature
429 std::cout <<
"\n\nValpro\n\n" << Valpro << std::endl;
430 std::cout <<
"\n\nVectPro\n\n" << VectPro << std::endl;
436 Vector v_kmax_curv(vect_pro(0, 1), vect_pro(1, 1), vect_pro(2, 1));
437 Vector v_kmin_curv(vect_pro(0, 0), vect_pro(1, 0), vect_pro(2, 0));
439 #ifdef DBG_calculate_curvature
440 std::cout <<
"\n\nValpro sorted\n\n" << Valpro << std::endl;
441 std::cout <<
"\n\nVectPro sorted\n\n" << VectPro << std::endl;
442 std::cout <<
"\nVKmaxCurv = " << VKmaxCurv[0] <<
" " << VKmaxCurv[1]
443 <<
" " << VKmaxCurv[2] << std::endl;
444 std::cout <<
"\nVKminCurv = " << VKminCurv[0] <<
" " << VKminCurv[1]
445 <<
" " << VKminCurv[2] << std::endl;
448 v_cm[*vi].VKmaxCurv = v_kmax_curv;
449 v_cm[*vi].VKminCurv = v_kmin_curv;
453 v_cm[*vi].KmaxCurv = valpro(0);
454 v_cm[*vi].KminCurv = valpro(1);
458 v_cm[*vi].KmaxCurv = -valpro(0);
459 v_cm[*vi].KminCurv = -valpro(1);
464 v_cm[*vi].VKmaxCurv = gt.NULL_VECTOR;
465 v_cm[*vi].VKminCurv = gt.NULL_VECTOR;
466 v_cm[*vi].KmaxCurv = 0;
467 v_cm[*vi].KminCurv = 0;
470 #ifdef DBG_calculate_curvature
471 std::cout <<
"\nv_cm[*vi].KmaxCurv = " << v_cm[*vi].KmaxCurv << std::endl;
472 std::cout <<
"\nv_cm[*vi].KminCurv = " << v_cm[*vi].KminCurv << std::endl;
475 min_nrm_min_curvature =
476 std::min< double >(min_nrm_min_curvature, v_cm[*vi].KminCurv);
477 max_nrm_min_curvature =
478 std::max< double >(max_nrm_min_curvature, v_cm[*vi].KminCurv);
480 min_nrm_max_curvature =
481 std::min< double >(min_nrm_max_curvature, v_cm[*vi].KmaxCurv);
482 max_nrm_max_curvature =
483 std::max< double >(max_nrm_max_curvature, v_cm[*vi].KmaxCurv);
488 std::cout <<
"Curvature warning: "
489 << skipped <<
" isolated vertices were skipped."
532 template<
typename HalfedgeGraph,
533 typename VertexCurvatureMap,
535 typename FaceNormalMap,
539 VertexCurvatureMap &v_cm,
541 const FaceNormalMap &f_nm,
544 double &min_nrm_min_curvature,
545 double &max_nrm_min_curvature,
546 double &min_nrm_max_curvature,
547 double &max_nrm_max_curvature)
549 GeometryTraits gt(g);
560 min_nrm_min_curvature,
561 max_nrm_min_curvature,
562 min_nrm_max_curvature,
563 max_nrm_max_curvature,