OMToolkit 1.0 The polygonal mesh processing tool.

- Main Page
- Related Pages
- Namespaces
- Classes
- Files

- File List
- File Members

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>
53 void OMCurvature<Mesh, Scalar>::getCurvatureDirections(VertexHandle vertex, Normal &max, Normal &min)
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>
189 bool OMCurvature<Mesh, Scalar>::calcMinCurvature()
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>
227 bool OMCurvature<Mesh, Scalar>::calcMaxCurvature()
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>
264 bool OMCurvature<Mesh, Scalar>::calcGaussCurvature()
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>
290 bool OMCurvature<Mesh, Scalar>::calcMeanCurvature()
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>
323 bool OMCurvature<Mesh, Scalar>::isObtuse(FaceHandle face, HalfedgeHandle &half_edge)
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>
345 Scalar OMCurvature<Mesh, Scalar>::getVoronoiAreaTriO(FaceHandle face, VertexHandle vertex)
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>
497 Scalar OMCurvature<Mesh, Scalar>::getVoronoiAreaTri(FaceHandle face, VertexHandle vertex)
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>
557 typename OMCurvature<Mesh, Scalar>::Normal OMCurvature<Mesh, Scalar>::getMeanOperator(VertexHandle vertex)
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>
599 Scalar OMCurvature<Mesh, Scalar>::getGaussOperator(VertexHandle vertex)
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>
628 Scalar OMCurvature<Mesh, Scalar>::getVoronoiArea(VertexHandle vertex)
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_