MEPP2 Project
curvature_cmdm.hpp
Go to the documentation of this file.
1 // Copyright (c) 2012-2020 University of Lyon and CNRS (France).
2 // All rights reserved.
3 //
4 // This file is part of MEPP2; you can redistribute it and/or modify
5 // it under the terms of the GNU Lesser General Public License as
6 // published by the Free Software Foundation; either version 3 of
7 // the License, or (at your option) any later version.
8 //
9 // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
10 // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
11 #pragma once
12 
13 #include <set>
14 #include <stack>
15 
16 #include <Eigen/Dense> // for Matrix::Identity()
17 #include <cmath> // for std::isnan()
18 
19 #include <boost/graph/graph_traits.hpp>
20 #include <boost/graph/properties.hpp>
21 
22 #include <CGAL/boost/graph/iterator.h> // for circulators
23 #include <CGAL/boost/graph/helpers.h> // for is_border_edge()
24 
25 #include "FEVV/Filters/Generic/Manifold/Curvature/eigen_val_vect.hpp"
26 
27 #include "FEVV/Wrappings/Geometry_traits.h"
28 
29 #include "FEVV/Tools/Math/MatrixOperations.hpp"
30 #include "FEVV/Operators/Geometry/AngleOperations.hpp"
31 #include "FEVV/Operators/Geometry/ClippingAndIntersection.hpp"
32 
33 #include "FEVV/Operators/Geometry/triangles.hpp"
34 
35 namespace FEVV {
36 namespace Filters {
37 
38 template< typename HalfedgeGraph,
39  typename GeometryTraits = FEVV::Geometry_traits< HalfedgeGraph > >
40 struct v_Curv_cmdm
41 {
42  double KminCurv;
43  double KmaxCurv;
44 
45  typename GeometryTraits::Vector
46  VKminCurv;
47  typename GeometryTraits::Vector
48  VKmaxCurv;
49 };
50 
51 template< typename HalfedgeGraph,
52  typename GeometryTraits = FEVV::Geometry_traits< HalfedgeGraph > >
53 double
54 fabs_cmdm(const v_Curv_cmdm< HalfedgeGraph > &input)
55 {
56  return std::max(::fabs(input.KminCurv), ::fabs(input.KmaxCurv));
57 }
58 
59 template< typename HalfedgeGraph,
60  typename PointMap,
61  typename GeometryTraits = FEVV::Geometry_traits< HalfedgeGraph > >
62 double
63 triangle_area_cmdm(
64  const HalfedgeGraph &g,
65  const PointMap &pm,
66  typename boost::graph_traits< HalfedgeGraph >::face_descriptor fd,
67  const GeometryTraits &gt)
68 {
69  typedef typename boost::property_traits< PointMap >::value_type Point3d;
70  //using Vector = typename GeometryTraits::Vector;
71 
72  CGAL::Halfedge_around_face_circulator< HalfedgeGraph > p_halfedge(
73  halfedge(fd, g), g);
74  Point3d p = get(pm, target(*p_halfedge, g));
75  Point3d q = get(pm, target(next(*p_halfedge, g), g));
76  Point3d r = get(pm, target(next(next(*p_halfedge, g), g), g));
77 
78  return FEVV::Operators::Geometry::triangle_area< GeometryTraits >(
79  p, q, r, gt);
80 }
81 
82 //********************************************
83 // 1-ring neighborhood curvature for a vertex
84 //********************************************
85 template< typename HalfedgeGraph,
86  typename PointMap,
87  typename FaceNormalMap,
88  typename GeometryTraits = FEVV::Geometry_traits< HalfedgeGraph >,
89  typename GraphTraits = boost::graph_traits< HalfedgeGraph >,
90  typename VertexDescriptor = typename GraphTraits::vertex_descriptor >
91 void
92 principal_curvature_per_vert_cmdm(const HalfedgeGraph &g,
93  const PointMap &pm,
94  const FaceNormalMap &f_nm,
95  VertexDescriptor vd,
96  double pp_matrix_sum[3][3],
97  const GeometryTraits &gt)
98 {
99  typedef typename GeometryTraits::Point Point3d;
100  using Vector = typename GeometryTraits::Vector;
101  typedef typename GraphTraits::face_descriptor face_descriptor;
102 
103  double area = 0;
104 
105  // iterate over all edges
106  CGAL::Halfedge_around_target_circulator< HalfedgeGraph > p_halfedge(vd, g),
107  done(p_halfedge);
108  do
109  {
110  // build edge vector and comput its norm
111  Point3d p1 = get(pm, target(*p_halfedge, g));
112  Point3d p2 = get(pm, source(*p_halfedge, g));
113 #if 1 // WAITING for AIF stuff - TODO : REMOVE THIS BLOC LATER...
114  Vector edge = gt.sub_p(p1, p2);
115 #else
116  Vector edge(p1[0] - p2[0], p1[1] - p2[1], p1[2] - p2[2]);
117 #endif
118  double len_edge = gt.length(edge);
119  if(len_edge == 0) // avoid divide by zero
120  continue;
121 
122  // compute (signed) angle between two incident faces, if exists
123  face_descriptor p_facet1 = face(*p_halfedge, g);
124  face_descriptor p_facet2 = face(opposite(*p_halfedge, g), g);
125  assert(p_facet1 != p_facet2);
126  if(p_facet1 == GraphTraits::null_face() ||
127  p_facet2 == GraphTraits::null_face())
128  continue; // border edge
129 
130  Vector normal1 = f_nm[p_facet1];
131  Vector normal2 = f_nm[p_facet2];
132 
133  if(std::isnan(normal1[0]) || std::isnan(normal1[1]) ||
134  std::isnan(normal1[2]) ||
135  std::isnan(normal2[0]) || std::isnan(normal2[1]) ||
136  std::isnan(normal2[2]))
137  continue; // Normal vector is set to be (NaN, NaN, NaN) if 2 edges of
138  // the face are collinear -- [YN]
139 
140  area += triangle_area_cmdm(g, pm, p_facet1, gt);
141 
142  // ---
143  Vector cp(normal1[1] * normal2[2] - normal1[2] * normal2[1],
144  normal1[2] * normal2[0] - normal1[0] * normal2[2],
145  normal1[0] * normal2[1] - normal1[1] * normal2[0]);
146  double ps_cp_xedge = cp[0] * edge[0] + cp[1] * edge[1] + cp[2] * edge[2];
147 
148  // ---
149  double sine =
150  ps_cp_xedge /
151  len_edge; // = (CGAL::cross_product(normal1,normal2)*edge)/len_edge;
152  // TODO BETTER
153  double beta = FEVV::Operators::Geometry::asin(sine);
154 
155  // compute edge * edge^t * coeff, and add it to current matrix
156  double p_vector_edge[3] = {edge[0], edge[1], edge[2]};
157  double pp_matrix[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
158  FEVV::Math::Matrix::Square::vector_times_transpose_mult(
159  p_vector_edge, pp_matrix, beta / len_edge);
160  FEVV::Math::Matrix::Square::add(pp_matrix, pp_matrix_sum);
161  } while(++p_halfedge != done);
162 
163  if(area != 0)
164  for(int i = 0; i < 3; i++)
165  for(int j = 0; j < 3; j++)
166  pp_matrix_sum[i][j] /= area;
167  else
168  std::cout << "Unable to compute curvature." << std::endl;
169 }
170 
171 //**********************************************
172 // geodesic neighborhood curvature for a vertex
173 //**********************************************
174 template< typename HalfedgeGraph,
175  typename PointMap,
176  typename FaceNormalMap,
177  typename GeometryTraits = FEVV::Geometry_traits< HalfedgeGraph >,
178  typename GraphTraits = boost::graph_traits< HalfedgeGraph >,
179  typename VertexDescriptor = typename GraphTraits::vertex_descriptor >
180 void
181 geodes_principal_curvature_per_vert_cmdm(const HalfedgeGraph &g,
182  const PointMap &pm,
183  const FaceNormalMap &f_nm,
184  VertexDescriptor vd,
185  double pp_matrix_sum[3][3],
186  double radius,
187  const GeometryTraits &gt)
188 {
189  typedef typename GeometryTraits::Point Point3d;
190  using Vector = typename GeometryTraits::Vector;
191  typedef typename GraphTraits::face_descriptor face_descriptor;
192 
193  std::set< VertexDescriptor > vertices;
194  Point3d o = get(pm, vd);
195  std::stack< VertexDescriptor > s;
196  s.push(vd);
197  vertices.insert(vd);
198  int iter = 0;
199  while(!s.empty())
200  {
201  VertexDescriptor v = s.top();
202  s.pop();
203  Point3d p = get(pm, v);
204  CGAL::Halfedge_around_target_circulator< HalfedgeGraph > h(v, g), done(h);
205  do
206  {
207  Point3d p1 = get(pm, target(*h, g));
208  Point3d p2 = get(pm, source(*h, g));
209 
210 #if 1 // WAITING for AIF stuff - TODO : REMOVE THIS BLOC LATER...
211  Vector vec = gt.sub_p(p2, p1);
212  Vector pm_o = gt.sub_p(p, o);
213 #else
214 
215  Vector vec(p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]);
216  Vector PmO(P[0] - O[0], P[1] - O[1], P[2] - O[2]);
217 #endif
218  double ps_vx_pm_o =
219  vec[0] * pm_o[0] + vec[1] * pm_o[1] + vec[2] * pm_o[2];
220 
221  // ---
222  if(v == vd || ps_vx_pm_o >= 0.0) // if (v==pVertex || vec * (P - O) >=
223  // 0.0) // TODO BETTER ; [>"=" (YN)]
224  {
225  bool normal1_nan = false, normal2_nan = false; // [YN]
226  bool pt_face1_superpose = false, pt_face2_superpose = false; // [YN]
227 
228  bool isect = FEVV::Operators::Geometry::sphere_clip_vector(
229  o, radius, p, vec, gt);
230 
231  if(!CGAL::is_border_edge(*h, g) && isect) // ["&& isect" (YN)]
232  {
233  double len_edge =
234  std::sqrt(vec[0] * vec[0] + vec[1] * vec[1] +
235  vec[2] * vec[2]); // std::sqrt(vec*vec); TODO BETTER
236 
237  // compute (signed) angle between two incident faces, if exists
238  face_descriptor p_facet1 = face(*h, g);
239  face_descriptor p_facet2 = face(opposite(*h, g), g);
240  assert(p_facet1 != p_facet2);
241 
242  std::vector< Point3d > face1_vertices; // [YN]
243  BOOST_FOREACH(VertexDescriptor vi, CGAL::vertices_around_face(*h, g))
244  face1_vertices.push_back(get(pm, vi));
245 
246  if(face1_vertices[0][0] == face1_vertices[1][0] &&
247  face1_vertices[1][0] == face1_vertices[2][0] &&
248  face1_vertices[0][0] == face1_vertices[2][0] &&
249  face1_vertices[0][1] == face1_vertices[1][1] &&
250  face1_vertices[1][1] == face1_vertices[2][1] &&
251  face1_vertices[0][1] == face1_vertices[2][1] &&
252  face1_vertices[0][2] == face1_vertices[1][2] &&
253  face1_vertices[1][2] == face1_vertices[2][2] &&
254  face1_vertices[0][2] == face1_vertices[2][2])
255  pt_face1_superpose = true;
256 
257  std::vector< Point3d > face2_vertices; // [YN]
258  BOOST_FOREACH(VertexDescriptor vi,
259  CGAL::vertices_around_face(opposite(*h, g), g))
260  face2_vertices.push_back(get(pm, vi));
261 
262  if(face2_vertices[0][0] == face2_vertices[1][0] &&
263  face2_vertices[1][0] == face2_vertices[2][0] &&
264  face2_vertices[0][0] == face2_vertices[2][0] &&
265  face2_vertices[0][1] == face2_vertices[1][1] &&
266  face2_vertices[1][1] == face2_vertices[2][1] &&
267  face2_vertices[0][1] == face2_vertices[2][1] &&
268  face2_vertices[0][2] == face2_vertices[1][2] &&
269  face2_vertices[1][2] == face2_vertices[2][2] &&
270  face2_vertices[0][2] == face2_vertices[2][2])
271  pt_face2_superpose = true;
272 
273  if(p_facet1 != GraphTraits::null_face() &&
274  p_facet2 != GraphTraits::null_face())
275  {
276  Vector normal1 = f_nm[p_facet1];
277  Vector normal2 = f_nm[p_facet2];
278 
279  // Check if face1 and face2 contain collinear edges [YN]
280  if(std::isnan(normal1[0]) || std::isnan(normal1[1]) ||
281  std::isnan(normal1[2]))
282  normal1_nan = true;
283  if(std::isnan(normal2[0]) || std::isnan(normal2[1]) ||
284  std::isnan(normal2[2]))
285  normal2_nan = true;
286 
287  // ---
288  double beta = 0;
289 
290  if(!pt_face2_superpose && !pt_face1_superpose)
291  {
292  // Check for common cases encountered in a mesh distorted by a
293  // uniform geometric quantization [YN]
294  if(!normal1_nan && !normal2_nan)
295  {
296  // Case 1: the 2 normals are colinear in opposite direction
297  double ps = normal1[0] * normal2[0] + normal1[1] * normal2[1] +
298  normal1[2] * normal2[2]; // scalar product
299  double acos = ps / std::sqrt((normal1[0] * normal1[0] +
300  normal1[1] * normal1[1] +
301  normal1[2] * normal1[2]) *
302  (normal2[0] * normal2[0] +
303  normal2[1] * normal2[1] +
304  normal2[2] * normal2[2]));
305 
306  if(acos == -1) // cos=180 degrees => the 2 normals oppose
307  {
308  //double beta =
309  // M_PI /
310  // 2; // set to maximum curvature achievable (in radian)
311  }
312  else
313  {
314  // Normal and most common case for all types of distortion
315  Vector cp(normal1[1] * normal2[2] - normal1[2] * normal2[1],
316  normal1[2] * normal2[0] - normal1[0] * normal2[2],
317  normal1[0] * normal2[1] - normal1[1] * normal2[0]);
318  double ps_cpx_v =
319  cp[0] * vec[0] + cp[1] * vec[1] + cp[2] * vec[2];
320 
321  double sine =
322  ps_cpx_v /
323  len_edge; // (CGAL::cross_product(normal1, normal2)*vec)
324  // / len_edge; // TODO BETTER
325  beta = FEVV::Operators::Geometry::asin(sine);
326  }
327  }
328  else if((normal1_nan && !normal2_nan) ||
329  (normal2_nan && !normal1_nan))
330  {
331  // Case 2: face2 contains collinear edges
332  if(!normal1_nan)
333  {
334  // consider all haledge of the source vertex
335  bool found_face = false;
336  Point3d destination = get(pm, target(*h, g));
337  VertexDescriptor t1 = source(*h, g);
338 
339  std::stack< VertexDescriptor > s2;
340  std::set< VertexDescriptor > s22;
341  s22.insert(t1);
342  s2.push(t1);
343  while(!s2.empty() && !found_face)
344  {
345  VertexDescriptor t = s2.top();
346  s2.pop();
347  CGAL::Halfedge_around_target_circulator< HalfedgeGraph > h2(
348  t, g),
349  done2(h2);
350  do
351  {
352  if(get(pm, source(*h2, g)) == destination)
353  {
354  face_descriptor newface2 = face(*h2, g);
355  if(newface2 != GraphTraits::null_face())
356  {
357  normal2 = f_nm[newface2];
358  if(!std::isnan(normal2[0]) &&
359  !std::isnan(normal2[1]) &&
360  !std::isnan(normal2[2]) && newface2 != p_facet1)
361  {
362  Vector cp(normal1[1] * normal2[2] -
363  normal1[2] * normal2[1],
364  normal1[2] * normal2[0] -
365  normal1[0] * normal2[2],
366  normal1[0] * normal2[1] -
367  normal1[1] * normal2[0]);
368  double ps_cpx_v = cp[0] * vec[0] + cp[1] * vec[1] +
369  cp[2] * vec[2];
370 
371  double sine = ps_cpx_v / len_edge;
372  beta = FEVV::Operators::Geometry::asin(sine);
373  found_face = true;
374  break;
375  }
376  else if(s22.find(source(*h2, g)) == s22.end())
377  {
378  s2.push(source(*h2, g));
379  s22.insert(source(*h2, g));
380  }
381  }
382  }
383  h2++;
384  } while(h2 != done2);
385  destination = get(pm, source(*h, g));
386  }
387  }
388 
389  // Case 3: face1 contains collinear edges
390  if(!normal2_nan)
391  {
392  // consider all haledge of the source vertex
393  Point3d destination = get(pm, target(*h, g));
394  bool found_face = false;
395  VertexDescriptor v1 = source(*h, g);
396 
397  std::stack< VertexDescriptor > s1;
398  std::set< VertexDescriptor > s11;
399  s11.insert(v1);
400  s1.push(v1);
401  while(!s1.empty() && !found_face)
402  {
403  VertexDescriptor s = s1.top();
404  s1.pop();
405  CGAL::Halfedge_around_target_circulator< HalfedgeGraph > h1(
406  s, g),
407  done1(h1);
408  do
409  {
410  if(get(pm, source(*h1, g)) == destination)
411  {
412  face_descriptor newface1 = face(*h1, g);
413  if(newface1 != GraphTraits::null_face())
414  {
415  normal1 = f_nm[newface1]; // face or opposite face??
416  if(!std::isnan(normal1[0]) &&
417  !std::isnan(normal1[1]) &&
418  !std::isnan(normal1[2]) && newface1 != p_facet2)
419  {
420  Vector cp(normal1[1] * normal2[2] -
421  normal1[2] * normal2[1],
422  normal1[2] * normal2[0] -
423  normal1[0] * normal2[2],
424  normal1[0] * normal2[1] -
425  normal1[1] * normal2[0]);
426  double ps_cpx_v = cp[0] * vec[0] + cp[1] * vec[1] +
427  cp[2] * vec[2];
428 
429  double sine = ps_cpx_v / len_edge;
430  beta = FEVV::Operators::Geometry::asin(sine);
431  found_face = true;
432  break;
433  }
434  else if(s11.find(source(*h1, g)) == s11.end())
435  {
436  s1.push(source(*h1, g));
437  s11.insert(source(*h1, g));
438  }
439  }
440  }
441  h1++;
442  } while(h1 != done1);
443  destination = get(pm, source(*h, g));
444  }
445  }
446  }
447  else
448  {
449  // Case 4: 2 colinear faces
450  if(normal2_nan &&
451  normal1_nan) // TODO BETTER: if(...) can be removed
452  beta = 0;
453  }
454 
455  // ---
456  // compute edge * edge^t * coeff, and add it to current matrix
457  double p_vector_edge[3] = {vec[0], vec[1], vec[2]};
458  double pp_matrix[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
459  FEVV::Math::Matrix::Square::vector_times_transpose_mult(
460  p_vector_edge, pp_matrix, beta / len_edge);
461  FEVV::Math::Matrix::Square::add(pp_matrix, pp_matrix_sum);
462  }
463  }
464  }
465 
466  if(!isect)
467  {
468  VertexDescriptor w = source(
469  *h, g); // previously : iterator here, NO !!! -> error from GL
470  // ??? Vertex_iterator w=h->opposite()->vertex();
471  assert(w != GraphTraits::null_vertex());
472  if(vertices.find(w) == vertices.end())
473  {
474  vertices.insert(w);
475  s.push(w);
476  }
477  }
478  }
479  } while(++h != done);
480 
481  iter++;
482  }
483 
484  double area = M_PI * radius * radius;
485 
486  for(int i = 0; i < 3; i++)
487  for(int j = 0; j < 3; j++)
488  pp_matrix_sum[i][j] /= area;
489 }
490 
491 
492 //#define DBG_calculate_curvature_cmdm
493 
494 
517 template< typename HalfedgeGraph,
518  typename VertexCurvatureMap,
519  typename PointMap,
520  typename FaceNormalMap,
521  typename GeometryTraits = FEVV::Geometry_traits< HalfedgeGraph > >
522 void
523 calculate_curvature_cmdm(const HalfedgeGraph &g,
524  VertexCurvatureMap &v_cm,
525  const PointMap &pm,
526  const FaceNormalMap &f_nm,
527  bool is_geod,
528  double radius,
529  double &min_nrm_min_curvature,
530  double &max_nrm_min_curvature,
531  double &min_nrm_max_curvature,
532  double &max_nrm_max_curvature,
533  const GeometryTraits &gt)
534 {
535  typedef boost::graph_traits< HalfedgeGraph > GraphTraits;
536  typedef typename GraphTraits::vertex_iterator vertex_iterator;
537  //typedef typename GraphTraits::vertex_descriptor Vertex_Descriptor;
538  using Vector = typename GeometryTraits::Vector;
539 
540  min_nrm_min_curvature = 100000;
541  max_nrm_min_curvature = -100000;
542 
543  min_nrm_max_curvature = 100000;
544  max_nrm_max_curvature = -100000;
545 
546  // ---
547 
548  std::cout << "Asking to calculate curvature" << std::endl;
549 #ifdef DBG_calculate_curvature_cmdm
550  std::cout << "Real radius = " << radius << std::endl;
551 #endif
552  auto iterator_pair = vertices(g); // vertices() returns a vertex_iterator pair
553  vertex_iterator vi = iterator_pair.first;
554  vertex_iterator vi_end = iterator_pair.second;
555 
556  // #pragma omp parallel for
557  for(; vi != vi_end; ++vi)
558  {
559  bool no_val_pro = false;
560 
561  double pp_matrix_sum[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
562 
563  if(is_geod == true) // geodesic neighborhood
564  geodes_principal_curvature_per_vert_cmdm<
565  HalfedgeGraph,
566  PointMap,
567  FaceNormalMap,
568  GeometryTraits,
569  GraphTraits,
570  typename GraphTraits::vertex_descriptor >(g, pm, f_nm, *vi,
571  pp_matrix_sum, radius, gt);
572 
573  else // 1-ring neighborhood
574  principal_curvature_per_vert_cmdm<
575  HalfedgeGraph,
576  PointMap,
577  FaceNormalMap,
578  GeometryTraits,
579  GraphTraits,
580  typename GraphTraits::vertex_descriptor >(g, pm, f_nm, *vi,
581  pp_matrix_sum, gt);
582 
583  // Eigen values/vectors
584  EigenMatrix cov_mat;
585  EigenvectorsType vect_pro;
586  EigenvalueType valpro;
587 
588  cov_mat(0, 0) = pp_matrix_sum[0][0];
589  cov_mat(0, 1) = pp_matrix_sum[0][1];
590  cov_mat(0, 2) = pp_matrix_sum[0][2];
591  cov_mat(1, 0) = pp_matrix_sum[1][0];
592  cov_mat(1, 1) = pp_matrix_sum[1][1];
593  cov_mat(1, 2) = pp_matrix_sum[1][2];
594  cov_mat(2, 0) = pp_matrix_sum[2][0];
595  cov_mat(2, 1) = pp_matrix_sum[2][1];
596  cov_mat(2, 2) = pp_matrix_sum[2][2];
597 
598 #ifdef DBG_calculate_curvature_cmdm
599  std::cout << "\n\nCovMat\n\n" << CovMat << std::endl;
600 #endif
601 
602  if(pp_matrix_sum[0][1] == 0 && pp_matrix_sum[0][2] == 0 &&
603  pp_matrix_sum[1][0] == 0 && pp_matrix_sum[1][2] == 0 &&
604  pp_matrix_sum[2][1] == 0 && pp_matrix_sum[2][0] == 0)
605  {
606  for(int i = 0; i < 3; i++)
607  valpro(i) = cov_mat(i, i);
608  vect_pro = EigenMatrix::Identity();
609  }
610  else
611  {
612  if(eigen_val_vect_compute(cov_mat, vect_pro, valpro) == -1)
613  {
614  no_val_pro = true;
615  }
616  }
617 
618 #ifdef DBG_calculate_curvature_cmdm
619  std::cout << "\n\nValpro\n\n" << Valpro << std::endl;
620  std::cout << "\n\nVectPro\n\n" << VectPro << std::endl;
621 #endif
622 
623  if(!no_val_pro)
624  {
625  eigen_val_vect_sort(valpro, vect_pro);
626  Vector v_kmax_curv(vect_pro(0, 1), vect_pro(1, 1), vect_pro(2, 1));
627  Vector v_kmin_curv(vect_pro(0, 0), vect_pro(1, 0), vect_pro(2, 0));
628 
629 #ifdef DBG_calculate_curvature_cmdm
630  std::cout << "\n\nValpro sorted\n\n" << Valpro << std::endl;
631  std::cout << "\n\nVectPro sorted\n\n" << VectPro << std::endl;
632  std::cout << "\nVKmaxCurv = " << VKmaxCurv[0] << " " << VKmaxCurv[1]
633  << " " << VKmaxCurv[2] << std::endl;
634  std::cout << "\nVKminCurv = " << VKminCurv[0] << " " << VKminCurv[1]
635  << " " << VKminCurv[2] << std::endl;
636 #endif
637 
638  v_cm[*vi].VKmaxCurv = v_kmax_curv;
639  v_cm[*vi].VKminCurv = v_kmin_curv;
640 
641  if(valpro(0) > 0)
642  {
643  v_cm[*vi].KmaxCurv = valpro(0);
644  v_cm[*vi].KminCurv = valpro(1);
645  }
646  else
647  {
648  v_cm[*vi].KmaxCurv = -valpro(0);
649  v_cm[*vi].KminCurv = -valpro(1);
650  }
651  }
652  else
653  {
654  v_cm[*vi].VKmaxCurv = gt.NULL_VECTOR;
655  v_cm[*vi].VKminCurv = gt.NULL_VECTOR;
656  v_cm[*vi].KmaxCurv = 0;
657  v_cm[*vi].KminCurv = 0;
658  }
659 
660 #ifdef DBG_calculate_curvature_cmdm
661  std::cout << "\nv_cm[*vi].KmaxCurv = " << v_cm[*vi].KmaxCurv << std::endl;
662  std::cout << "\nv_cm[*vi].KminCurv = " << v_cm[*vi].KminCurv << std::endl;
663 #endif
664 
665  min_nrm_min_curvature =
666  std::min< double >(min_nrm_min_curvature, v_cm[*vi].KminCurv);
667  max_nrm_min_curvature =
668  std::max< double >(max_nrm_min_curvature, v_cm[*vi].KminCurv);
669 
670  min_nrm_max_curvature =
671  std::min< double >(min_nrm_max_curvature, v_cm[*vi].KmaxCurv);
672  max_nrm_max_curvature =
673  std::max< double >(max_nrm_max_curvature, v_cm[*vi].KmaxCurv);
674  }
675 
676  // ---
677 
678  /*
679  int cpt = 0;
680  auto iterator_pair_log = vertices(g); // vertices() returns a
681  vertex_iterator pair vertex_iterator vi_log = iterator_pair_log.first;
682  vertex_iterator vend_log = iterator_pair_log.second; for (; vi_log !=
683  vend_log; ++vi_log)
684  {
685  printf("cpt: %4d - Kmin: %.3lf - Kmax: %.3lf | VKmin: (%.3lf, %.3lf,
686  %.3lf) | VKmax: (%.3lf, %.3lf, %.3lf)\n", cpt, v_cm[*vi_log].KminCurv,
687  v_cm[*vi_log].KmaxCurv, v_cm[*vi_log].VKminCurv[0],
688  v_cm[*vi_log].VKminCurv[1], v_cm[*vi_log].VKminCurv[2],
689  v_cm[*vi_log].VKmaxCurv[0], v_cm[*vi_log].VKmaxCurv[1],
690  v_cm[*vi_log].VKmaxCurv[2]); cpt++;
691  }
692  */
693 }
694 
715 template< typename HalfedgeGraph,
716  typename VertexCurvatureMap,
717  typename PointMap,
718  typename FaceNormalMap,
719  typename GeometryTraits = FEVV::Geometry_traits< HalfedgeGraph > >
720 void
721 calculate_curvature_cmdm(const HalfedgeGraph &g,
722  VertexCurvatureMap &v_cm,
723  const PointMap &pm,
724  const FaceNormalMap &f_nm,
725  bool is_geod,
726  double radius,
727  double &min_nrm_min_curvature,
728  double &max_nrm_min_curvature,
729  double &min_nrm_max_curvature,
730  double &max_nrm_max_curvature)
731 {
732  GeometryTraits gt(g);
733  calculate_curvature_cmdm< HalfedgeGraph,
734  VertexCurvatureMap,
735  PointMap,
736  FaceNormalMap,
737  GeometryTraits >(g,
738  v_cm,
739  pm,
740  f_nm,
741  is_geod,
742  radius,
743  min_nrm_min_curvature,
744  max_nrm_min_curvature,
745  min_nrm_max_curvature,
746  max_nrm_max_curvature,
747  gt);
748 }
749 
750 } // namespace Filters
751 } // namespace FEVV
FEVV::DataStructures::AIF::vertices
std::pair< typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::vertex_iterator, typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::vertex_iterator > vertices(const FEVV::DataStructures::AIF::AIFMesh &sm)
Returns the iterator range of the vertices of the mesh.
Definition: Graph_traits_aif.h:172
EigenMatrix
Eigen::Matrix3d EigenMatrix
Definition: eigen_val_vect.hpp:19
EigenvectorsType
Eigen::Matrix3d EigenvectorsType
Definition: eigen_val_vect.hpp:20
EigenvalueType
Eigen::Vector3d EigenvalueType
Definition: eigen_val_vect.hpp:21
eigen_val_vect_sort
void eigen_val_vect_sort(EigenvalueType &ad, EigenvectorsType &u)
Definition: eigen_val_vect.hpp:62
FEVV::Filters::v_Curv_cmdm::VKminCurv
GeometryTraits::Vector VKminCurv
The minimum curvature direction.
Definition: curvature_cmdm.hpp:46
FEVV::DataStructures::AIF::next
boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::halfedge_descriptor next(typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::halfedge_descriptor h, const FEVV::DataStructures::AIF::AIFMesh &sm)
Returns the next halfedge around its face.
Definition: Graph_traits_aif.h:599
Vector
AIFMesh::Vector Vector
Definition: Graph_properties_aif.h:22
FEVV::Math::Matrix::Square::vector_times_transpose_mult
static void vector_times_transpose_mult(const ElementType *p_vector, ElementType pp_matrix[][DIM], ElementType coeff)
Compute coef * V x V^t (V is a vertical DIM-D vector)
Definition: MatrixOperations.hpp:976
FEVV::Filters::triangle_area_cmdm
double triangle_area_cmdm(const HalfedgeGraph &g, const PointMap &pm, typename boost::graph_traits< HalfedgeGraph >::face_descriptor fd, const GeometryTraits &gt)
Definition: curvature_cmdm.hpp:63
FEVV::Filters::calculate_curvature_cmdm
void calculate_curvature_cmdm(const HalfedgeGraph &g, VertexCurvatureMap &v_cm, const PointMap &pm, const FaceNormalMap &f_nm, bool is_geod, double radius, double &min_nrm_min_curvature, double &max_nrm_min_curvature, double &min_nrm_max_curvature, double &max_nrm_max_curvature, const GeometryTraits &gt)
Calculate the curvature for a mesh.
Definition: curvature_cmdm.hpp:523
eigen_val_vect.hpp
FEVV::DataStructures::AIF::opposite
boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::halfedge_descriptor opposite(typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::halfedge_descriptor h, const FEVV::DataStructures::AIF::AIFMesh &sm)
Returns the halfedge with source and target swapped.
Definition: Graph_traits_aif.h:625
FEVV::Geometry_traits< HalfedgeGraph >
Point
AIFMesh::Point Point
Definition: Graph_properties_aif.h:21
FEVV::DataStructures::AIF::source
boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::vertex_descriptor source(typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::edge_descriptor e, const FEVV::DataStructures::AIF::AIFMesh &)
Returns the source vertex of e.
Definition: Graph_traits_aif.h:387
FEVV::DataStructures::AIF::edge
boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::edge_descriptor edge(typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::halfedge_descriptor h, const FEVV::DataStructures::AIF::AIFMesh &sm)
Returns the edge corresponding to h and opposite(h).
Definition: Graph_traits_aif.h:345
FEVV::Filters::fabs_cmdm
double fabs_cmdm(const v_Curv_cmdm< HalfedgeGraph > &input)
Definition: curvature_cmdm.hpp:54
FEVV::Filters::v_Curv_cmdm::VKmaxCurv
GeometryTraits::Vector VKmaxCurv
The minimum curvature direction.
Definition: curvature_cmdm.hpp:48
FEVV::Operators::Geometry::sphere_clip_vector
static bool sphere_clip_vector(const typename GeometryTraits::Point &center, double r, const typename GeometryTraits::Point &p, typename GeometryTraits::Vector &v, const GeometryTraits &gt)
Compute the intersection of a sphere (center + radius) with a ray/line (starting point + direction ve...
Definition: ClippingAndIntersection.hpp:39
FEVV::get
FEVV::PCLPointCloudPointMap::value_type get(const FEVV::PCLPointCloudPointMap &pm, FEVV::PCLPointCloudPointMap::key_type key)
Specialization of get(point_map, key) for PCLPointCloud.
Definition: Graph_properties_pcl_point_cloud.h:117
FEVV
Interfaces for plugins These interfaces will be used for different plugins.
Definition: Assert.h:16
FEVV::Filters::principal_curvature_per_vert_cmdm
void principal_curvature_per_vert_cmdm(const HalfedgeGraph &g, const PointMap &pm, const FaceNormalMap &f_nm, VertexDescriptor vd, double pp_matrix_sum[3][3], const GeometryTraits &gt)
Definition: curvature_cmdm.hpp:92
FEVV::Filters::geodes_principal_curvature_per_vert_cmdm
void geodes_principal_curvature_per_vert_cmdm(const HalfedgeGraph &g, const PointMap &pm, const FaceNormalMap &f_nm, VertexDescriptor vd, double pp_matrix_sum[3][3], double radius, const GeometryTraits &gt)
Definition: curvature_cmdm.hpp:181
FEVV::Filters::v_Curv_cmdm::KmaxCurv
double KmaxCurv
The minimum curvature.
Definition: curvature_cmdm.hpp:43
AngleOperations.hpp
FEVV::DataStructures::AIF::halfedge
boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::halfedge_descriptor halfedge(typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::vertex_descriptor v, const FEVV::DataStructures::AIF::AIFMesh &sm)
Returns a halfedge with target v.
Definition: Graph_traits_aif.h:296
FEVV::DataStructures::AIF::target
boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::vertex_descriptor target(typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::edge_descriptor e, const FEVV::DataStructures::AIF::AIFMesh &)
Returns the target vertex of e.
Definition: Graph_traits_aif.h:400
Geometry_traits.h
FEVV::Filters::v_Curv_cmdm
Definition: curvature_cmdm.hpp:41
triangles.hpp
FEVV::DataStructures::AIF::AIFVector
Definition: AIFProperties.h:173
ClippingAndIntersection.hpp
FEVV::Operators::Geometry::acos
T acos(T cosine)
Safe call to the std::acos function.
Definition: AngleOperations.hpp:64
msdm2::vertex_descriptor
boost::graph_traits< MeshT >::vertex_descriptor vertex_descriptor
Definition: msdm2_surfacemesh.h:33
FEVV::Filters::fabs
double fabs(const v_Curv< HalfedgeGraph > &input)
Definition: curvature.hpp:54
FEVV::DataStructures::AIF::face
boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::face_descriptor face(typename boost::graph_traits< FEVV::DataStructures::AIF::AIFMesh >::halfedge_descriptor h, const FEVV::DataStructures::AIF::AIFMesh &)
Returns the face incident to halfedge h.
Definition: Graph_traits_aif.h:664
FEVV::Filters::v_Curv_cmdm::KminCurv
double KminCurv
The minimum curvature.
Definition: curvature_cmdm.hpp:42
Point3d
EnrichedPolyhedron::Point_3 Point3d
Definition: boolops_cpolyhedron_builder.hpp:27
eigen_val_vect_compute
int eigen_val_vect_compute(const EigenMatrix &a, EigenvectorsType &u, EigenvalueType &ad)
Definition: eigen_val_vect.hpp:33
FEVV::Operators::Geometry::asin
static T asin(T sine)
Safe call to the std::asin function.
Definition: AngleOperations.hpp:31
MatrixOperations.hpp
FEVV::Math::Matrix::Square::add
static void add(const ElementType p_matrix[][DIM], ElementType p_matrix_sum[][DIM])
add two matrices
Definition: MatrixOperations.hpp:1038
Generated by doxygen 1.8.20