OMToolkit  1.0
The polygonal mesh processing tool.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
OMCurvature.hxx
Go to the documentation of this file.
1 //==============================================================================
14 #ifndef _OM_CURVATURE_HXX_
15 #define _OM_CURVATURE_HXX_
16 
18 // constants for computations
19 const double M_PIHALF = M_PI/2;
20 const double M_TWOPI = M_PI*2;
21 
23 
25 // Constructor - create an instance of this class
27 template <class Mesh, class Scalar>
28 OMCurvature<Mesh, Scalar>::OMCurvature(Mesh *mesh, OpenMesh::VPropHandleT<Normal> &curvVectorHandle, OpenMesh::VPropHandleT<Scalar> &curvMagHandle)
29 {
30  m_mesh = mesh;
31  m_mesh->request_face_normals();
32  m_mesh->request_vertex_normals();
33  m_mesh->update_normals();
34 
35  // initialize pre-computed voronoi area
36  m_mesh->add_property(m_areas, "<VoronoiA>");
37  Mesh::FIter end = m_mesh->faces_end();
38  for (Mesh::FIter face = m_mesh->faces_begin(); face!= end; ++face)
39  m_mesh->property(m_areas, face) = AreasVector(0.0, 0.0, 0.0);
40 
41  m_curvVectorHandle = curvVectorHandle;
42  m_curvMagHandle = curvMagHandle;
43 }
44 
46 // Function returns principal curvature directions for specified vertex
47 // Computation from paper of Gabriel Taubin - Estimating the tenspr of curvature of a surface from a polyhedral approximation
48 // @param vertex Vertex handle
49 // @param max maximum curvature direction vector
50 // @param min minimum curvature direction vector
52 template <class Mesh, class Scalar>
54 {
55  // helper variables (for not allocate them in the loop)
56  Normal normal = m_mesh->normal(vertex);
57  Normal inVector;
58  Normal outVector;
59  HalfedgeHandle he;
60 
61  Scalar Kij;
62  Scalar w;
63  FaceHandle face1;
64  FaceHandle face2;
65  Matrix3x3 TijAux;
66  ColumnVector3 Tij;
67  RowVector3 TijT;
68 
69  // Helper normals in matrix structure
70  ColumnVector3 normalVector;
71  RowVector3 normalVectorT;
72  normalVector(0, 0) = normal[0];
73  normalVector(1, 0) = normal[1];
74  normalVector(2, 0) = normal[2];
75  normalVectorT = normalVector.transpose();
76  ColumnVector3 edgeVectorIn(3, 1);
77  ColumnVector3 edgeVectorOut(3, 1);
78 
79  // identity matrix
80  Matrix3x3 identity;
81  identity.setIdentity();
82 
83  Matrix3x3 M;
84  M = M.setZero();
85 
86  Matrix3x3 normalMult;
87  normalMult = normalVector * normalVectorT;
88  Scalar vecLength;
89 
90  Scalar sumArea = getVoronoiArea(vertex) * 2;
91  // for each outgoing half edges
92  for (Mesh::VOHIter outgoing = m_mesh->voh_begin(vertex); outgoing; ++outgoing)
93  {
94  // retrieve vector of a half edge and put it into matrices
95  m_mesh->calc_edge_vector( outgoing, outVector);
96 
97  edgeVectorIn(0, 0) = -outVector[0];
98  edgeVectorIn(1, 0) = -outVector[1];
99  edgeVectorIn(2, 0) = -outVector[2];
100 
101  edgeVectorOut(0, 0) = outVector[0];
102  edgeVectorOut(1, 0) = outVector[1];
103  edgeVectorOut(2, 0) = outVector[2];
104 
105  // compute Tij
106 
107  Tij = (identity - normalMult) * edgeVectorIn;
108  Tij.normalize();
109  TijT = Tij.transpose();
110 
111  // compute Kij
112  vecLength = outVector.length();
113  Kij = ((normalVectorT * edgeVectorOut)(0) * 2) / (vecLength * vecLength);
114 
115  // compute weigth
116  face1 = m_mesh->face_handle(outgoing);
117  face2 = m_mesh->opposite_face_handle(outgoing);
118 
119  w = 0;
120  if (m_mesh->is_valid_handle(face1))
121  {
122  if (isObtuse(face1, he))
123  w += getVoronoiAreaTriO(face1, vertex);
124  else w += getVoronoiAreaTri(face1, vertex);
125  }
126  if (m_mesh->is_valid_handle(face2))
127  {
128  if (isObtuse(face2, he))
129  w += getVoronoiAreaTriO(face2, vertex);
130  else w += getVoronoiAreaTri(face2, vertex);
131  }
132  w /= sumArea;
133  // add to M
134  M += w * Kij * Tij * TijT;
135  }
136 
137  // comp householder matrix (restrict matrix to normal plane only)
138  ColumnVector3 E;
139  E(0, 0) = 1.0;
140  E(1, 0) = 0.0;
141  E(2, 0) = 0.0;
142 
143  ColumnVector3 W;
144  RowVector3 WT;
145  if ((E - normalVector).norm() > (E + normalVector).norm())
146  W = E - normalVector;
147  else
148  W = E + normalVector;
149 
150  vecLength = W.norm();
151  W /= vecLength;
152  WT = W.transpose();
153 
154  Matrix3x3 Q = identity - 2 * (W * WT);
155 
156  Matrix3x3 result;
157  result = Q.transpose() * M * Q;
158 
159  Matrix2x2 subM;
160 
161  // save sub matrix of householder transformed matrix M into 2x2 matrix
162  subM(0, 0) = result(1, 1);
163  subM(0, 1) = result(1, 2);
164  subM(1, 0) = result(2, 1);
165  subM(1, 1) = result(2, 2);
166 
167  // compute Givens rotation angle (angle necessary to null element (1, 0))
168  Scalar r = sqrt(subM(0, 0) * subM(0, 0) + subM(1, 0)*subM(1, 0));
169  Scalar c = subM(0, 0)/r;
170  Scalar s = subM(1, 0)/r;
171 
172  // rotate T1 and T2 vectors and compute min and max curvature
173  Normal T1 = Normal(Q(0, 1), Q(1, 1), Q(2, 1));
174  Normal T2 = Normal(Q(0, 2), Q(1, 2), Q(2, 2));
175  //max = c * T1.operator*(c) - s*T2;
176  //min = s * T1 + c*T2;
177  c = cos(c);
178  s = sin(s);
179  max = T1.operator*(c) - T2.operator*(s);
180  min = T1.operator*(s) + T2.operator*(c);
181 }
182 
184 // Calculates minimal curvature and for all mesh vertices
185 // @return True, if calculated successfully
186 // @todo Need to implement directions
188 template <class Mesh, class Scalar>
190 {
191  if (!m_mesh->IsTriMesh) return false;
192 
193  Scalar delta;
194  Scalar aux;
195  Normal min, max;
196  Mesh::VertexIter end = m_mesh->vertices_end();
197 
198  // for each vertex of the model
199  for (Mesh::VertexIter it = m_mesh->vertices_begin(); it != m_mesh->vertices_end(); ++it)
200  {
201  // calculate only if it is not boundary
202  if (!m_mesh->is_boundary(it))
203  {
204  aux = getMeanOperator(it).length()/2;
205  delta = aux * aux - getGaussOperator(it);
206  m_mesh->property(m_curvMagHandle, it) = aux - sqrt(delta);
207  getCurvatureDirections(it, max, min);
208  min.normalize();
209  m_mesh->property(m_curvVectorHandle, it) = min;
210  }
211  // else return zeros
212  else
213  {
214  m_mesh->property(m_curvMagHandle, it) = 0.0;
215  m_mesh->property(m_curvVectorHandle, it) = Normal(0.0, 0.0, 0.0);
216  }
217  }
218  return true;
219 }
220 
222 // Calculates maximal curvature and for all mesh vertices
223 // @return True, if calculated successfully
224 // @todo Need to implement directions
226 template <class Mesh, class Scalar>
228 {
229  if (!m_mesh->IsTriMesh) return false;
230 
231  Scalar delta;
232  Scalar aux;
233  Normal min, max;
234 
235  Mesh::VertexIter end = m_mesh->vertices_end();
236  // for each vertex of the model
237  for (Mesh::VertexIter it = m_mesh->vertices_begin(); it != m_mesh->vertices_end(); ++it)
238  {
239  // calculate only if it is not boundary
240  if (!m_mesh->is_boundary(it))
241  {
242  aux = getMeanOperator(it).length()/2;
243  delta = aux * aux - getGaussOperator(it);
244  m_mesh->property(m_curvMagHandle, it) = aux + sqrt(delta);
245  getCurvatureDirections(it, max, min);
246  max.normalize();
247  m_mesh->property(m_curvVectorHandle, it) = max;
248  }
249  // else return zeros
250  else
251  {
252  m_mesh->property(m_curvMagHandle, it) = 0.0;
253  m_mesh->property(m_curvVectorHandle, it) = Normal(0.0, 0.0, 0.0);
254  }
255  }
256  return true;
257 }
258 
260 // Calculates Gauss curvature and for all mesh vertices
261 // @return True, if calculated successfully
263 template <class Mesh, class Scalar>
265 {
266  if (!m_mesh->IsTriMesh) return false;
267 
268  Mesh::VertexIter end = m_mesh->vertices_end();
269  // for each vertex of the model
270  for (Mesh::VertexIter it = m_mesh->vertices_begin(); it != m_mesh->vertices_end(); ++it)
271  {
272  // calculate only if it is not boundary
273  if (!m_mesh->is_boundary(it))
274  m_mesh->property(m_curvMagHandle, it) = getGaussOperator(it);
275  // else return zero
276  else
277  m_mesh->property(m_curvMagHandle, it) = 0.0;
278 
279  // Gauss curvature has not a normal
280  m_mesh->property(m_curvVectorHandle, it) = Normal(0.0, 0.0, 0.0);
281  }
282  return true;
283 }
284 
286 // Calculates mean curvature and its direction for all mesh vertices
287 // @return True, if calculated successfully
289 template <class Mesh, class Scalar>
291 {
292  if (!m_mesh->IsTriMesh) return false;
293 
294  Mesh::VertexIter end = m_mesh->vertices_end();
295  // for each vertex of the model
296  for (Mesh::VertexIter it = m_mesh->vertices_begin(); it != m_mesh->vertices_end(); ++it)
297  {
298  // calculate only if it is not boundary
299  if (!m_mesh->is_boundary(it))
300  {
301  Normal vector = getMeanOperator(it);
302  m_mesh->property(m_curvMagHandle, it) = vector.length() / 2.0;
303  vector.normalize();
304  m_mesh->property(m_curvVectorHandle, it) = vector;
305  }
306  // else return zeros
307  else
308  {
309  m_mesh->property(m_curvMagHandle, it) = 0.0;
310  m_mesh->property(m_curvVectorHandle, it) = Normal(0.0, 0.0, 0.0);
311  }
312  }
313  return true;
314 }
315 
317 // Predicate, is triangle obtuse? (Is one of its angles obtuse?)
318 // @param face Specified triangle
319 // @param half_edge Half edge handle pointing on found obtuse angle (if not obtuse, it is not specified)
320 // @return true, if triangle has an obtuse angle
322 template <class Mesh, class Scalar>
324 {
325  Mesh::FEIter it;
326 
327  for(it = m_mesh->fe_begin(face); it; ++it)
328  {
329  if (m_mesh->calc_sector_angle(it.current_halfedge_handle()) > M_PIHALF)
330  {
331  half_edge = it.current_halfedge_handle();
332  return true;
333  }
334  }
335  return false;
336 }
337 
339 // Computes voronoi area of obtuse triangle from given vertex
340 // @param face Given triangle
341 // @param face Vertex handle
342 // @return Area of voronoi region
344 template <class Mesh, class Scalar>
346 {
347  Mesh::FHIter it = m_mesh->fh_begin(face);
348  m_mesh->fv_begin(face);
349  unsigned int i = 0;
350  // find a halfedge, which points to a vertex and is part of this face (which is not always true with registered edge)
351  for (; it; ++it)
352  {
353  if (m_mesh->to_vertex_handle(it) == vertex)
354  break;
355  ++i;
356  }
357  // if there is a problem (it should't, but.....)
358  if (m_mesh->to_vertex_handle(it) != vertex) return 0.0;
359 
360  // if there is already computed area
361  Scalar result = m_mesh->property(m_areas, face)[i];
362  if (result != 0) return result;
363 
364  // pre-compute next and previous half edges
365  HalfedgeHandle next = ++it;
366  HalfedgeHandle prev = ++it;
367  ++it;
368 
369  Scalar angleIt = m_mesh->calc_sector_angle(it);
370  Scalar angleNext = m_mesh->calc_sector_angle(next);
371 
372  /*
373  Supose this triangle, we have three possibilities:
374  it points to A, an obtuse angle
375  it points to B, successor of A
376  it points to C, predecessor of A
377 
378  C
379  \___/__
380  \
381  \
382  \
383  _\
384  \_____________________\__
385  A B
386 
387  */
388 
389  // we most find these three cases and compute different numbers
390  if (angleIt < M_PIHALF)
391  {
392  // if computing point B
393  if (angleNext < M_PIHALF)
394  {
395  Scalar lenght = m_mesh->calc_edge_length(it);
396  Scalar result = (lenght * lenght * tan(angleIt))/8.0;
397  m_mesh->property(m_areas, face)[i] = result;
398  return result;
399  }
400  // if computing point C
401  else
402  {
403  Scalar lenght = m_mesh->calc_edge_length(next);
404  Scalar result = (lenght * lenght * tan(angleIt))/8.0;
405  m_mesh->property(m_areas, face)[i] = result;
406  return result;
407  }
408  }
409  // if computing point A
410  else
411  {
412  Scalar result = m_mesh->calc_sector_area(it);
413  Scalar lenghtB = m_mesh->calc_edge_length(next);
414  result -= (lenghtB * lenghtB * tan(angleNext))/8.0;
415 
416  Scalar lenghtC = m_mesh->calc_edge_length(it);
417  result -= ((lenghtC * lenghtC * tan(m_mesh->calc_sector_angle(prev)))/8.0);
418  m_mesh->property(m_areas, face)[i] = result;
419  return result;
420  }
421 }
422 
429 //template <class Mesh, class Scalar>
430 //Scalar OMCurvature<Mesh, Scalar>::getVoronoiAreaTriO(VertexHandle *vertices)
431 //{
432 // Mesh::Normal pointing = vertices[0] - vertices[2];
433 // Mesh::Normal next = vertices[1] - vertices[0];
434 // Mesh::Normal prev = vertices[2] - vertices[1];
435 //
436 // edges[1] = vertices[2] - vertices[1];
437 // edges[2] = vertices[0] - vertices[2];
438 //
439 // Scalar angleIt = (pointing | next) / (pointing.norm() * pointing.norm());
440 // Scalar angleNext = (next | prev) / (next.norm() * prev.norm())
441 //
442 // /*
443 // Supose this triangle, we have three possibilities:
444 // it points to A, an obtuse angle
445 // it points to B, successor of A
446 // it points to C, predecessor of A
447 //
448 // C
449 // \___/__
450 // \
451 // \
452 // \
453 // _\
454 // \_____________________\__
455 // A B
456 //
457 // */
458 //
459 // // we most find these three cases and compute different numbers
460 // if (angleIt < M_PIHALF)
461 // {
462 // // if computing point B
463 // if (angleNext < M_PIHALF)
464 // {
465 // Scalar lenght = pointing.norm();
466 // Scalar result = (lenght * lenght * tan(angleIt))/8.0;
467 // return result;
468 // }
469 // // if computing point C
470 // else
471 // {
472 // Scalar lenght = next.norm();
473 // Scalar result = (lenght * lenght * tan(angleIt))/8.0;
474 // return result;
475 // }
476 // }
477 // // if computing point A
478 // else
479 // {
480 // Scalar p = (pointing.norm() + next.norm() + prev.norm()) / 2.0;
481 // Scalar result = sqrt(p * (p - pointing.norm()) * (p * next.norm()) * (p * prev.norm()));
482 // Scalar lenghtB = next.norm();
483 // result -= (lenghtB * lenghtB * tan(angleNext))/8.0;
484 //
485 // Scalar lenghtC = pointing.norm();
486 // result -= ((lenghtC * lenghtC * tan((pointing | prev) / (pointing.norm() * prev.norm())))/8.0);
487 // return result;
488 // }
489 //}
491 // Computes voronoi area of non-obtuse triangle from given vertex
492 // @param face Given triangle
493 // @param face Vertex handle
494 // @return Area of voronoi region
496 template <class Mesh, class Scalar>
498 {
499  Mesh::FHIter it = m_mesh->fh_begin(face);
500  unsigned int i = 0;
501  // find a halfedge, which points to a vertex and is part of this face (which is not always true with registered edge)
502  for (; it; ++it)
503  {
504  if (m_mesh->to_vertex_handle(it) == vertex)
505  break;
506  ++i;
507  }
508  if (m_mesh->to_vertex_handle(it) != vertex) return 0.0;
509 
510  // if there is already computed area
511  Scalar result = m_mesh->property(m_areas, face)[i];
512  if (result != 0) return result;
513 
514  // pre-compute next and previous half edges
515  Mesh::FHIter hedge = it;
516  HalfedgeHandle next = ++it;
517  HalfedgeHandle prev = ++it;
518  /*
519  We have this situation:
520  computed vertex
521  A
522  /\_
523  / \
524  next / \ hedge
525  / \
526  |/ \
527  /________\_\
528  C B
529  prev
530 
531  Sector angle is angle between hedge and successor
532  computing: (1/8) * (lenght(AC) * cot(B) + lenght(AB) * cot(C))
533  */
534 
535  Scalar auxLenght = m_mesh->calc_edge_length(hedge);
536  Scalar angle = m_mesh->calc_sector_angle(next);
537 
538  if (angle > 0 && angle < M_PIHALF)
539  result = auxLenght * auxLenght * 1/tan(angle);
540 
541  angle= m_mesh->calc_sector_angle(prev);
542  auxLenght = m_mesh->calc_edge_length(next);
543  if (angle > 0 && angle < M_PIHALF)
544  result += auxLenght * auxLenght * 1/tan(angle);
545 
546  result /= 8.0;
547  m_mesh->property(m_areas, face)[i] = result;
548  return result;
549 }
550 
552 // Computes mean curvature operator (vector, which direction equals mean curvature direction, its norm equals 2*mean curvature) for specified vertex
553 // @param vertex Vertex handle
554 // @return mean curvature operator
556 template <class Mesh, class Scalar>
558 {
559  HalfedgeHandle opp;
560  Normal result(0.0, 0.0, 0.0);
561 
562  Normal vector;
563  Scalar aux;
564  Scalar tang;
565 
566  // for each outgoing halfedge
567  for (Mesh::VertexOHalfedgeIter it = m_mesh->voh_begin(vertex); it; ++it)
568  {
569  // compute cotangents
570  opp = m_mesh->opposite_halfedge_handle(it);
571  tang = tan(m_mesh->calc_sector_angle(m_mesh->next_halfedge_handle(it)));
572  if (tang != 0)
573  aux = 1/tang;
574  tang = tan(m_mesh->calc_sector_angle(m_mesh->next_halfedge_handle(opp)));
575  if (tang != 0)
576  aux += 1/tang;
577 
578  m_mesh->calc_edge_vector(opp, vector);
579  vector *= aux;
580  result += (Normal)vector;
581  }
582 
583  // and normalize with voronoi area region
584  aux = getVoronoiArea(vertex);
585  if (aux != 0)
586  {
587  result /= (2*aux);
588  return result;
589  }
590  else return Normal(0.0, 0.0, 0.0);
591 }
592 
594 // Computes Gauss curvature value for specified vertex
595 // @param vertex Vertex handle
596 // @return Gauss curvature magnitude
598 template<class Mesh, class Scalar>
600 {
601  HalfedgeHandle opp;
602  Normal result(0.0, 0.0, 0.0);
603 
604  Normal vector;
605  Scalar aux = 0.0;
606  Scalar area;
607 
608  // for each outgoing halfedge
609  for (Mesh::VertexIHalfedgeIter it = m_mesh->vih_begin(vertex); it; ++it)
610  aux += m_mesh->calc_sector_angle(it);
611 
612  aux = M_TWOPI - aux;
613  area = getVoronoiArea(vertex);
614  if (area != 0)
615  {
616  return aux/area;
617  }
618  else return 0;
619 
620 }
621 
623 // Computes mean curvature operator (vector, which direction equals mean curvature direction, its norm equals 2*mean curvature) for specified vertex
624 // @param vertex Vertex handle
625 // @return mean curvature operator
627 template<class Mesh, class Scalar>
629 {
630  Mesh::VFIter it, end;
631  end = m_mesh->vf_end(vertex);
632 
633  HalfedgeHandle hedge;
634  Scalar result = 0.0;
635  for (it = m_mesh->vf_begin(vertex); it; ++it)
636  {
637  if (!isObtuse(it, hedge))
638  {
639  result += getVoronoiAreaTri(it, vertex);
640  }
641  else
642  {
643  result += getVoronoiAreaTriO(it, vertex);
644  }
645  }
646  return result;
647 }
648 
649 #endif // _OM_CURVATURE_HXX_