Developer Documentation
CurvatureT_impl.hh
1 /*===========================================================================*\
2 * *
3 * OpenFlipper *
4  * Copyright (c) 2001-2015, RWTH-Aachen University *
5  * Department of Computer Graphics and Multimedia *
6  * All rights reserved. *
7  * www.openflipper.org *
8  * *
9  *---------------------------------------------------------------------------*
10  * This file is part of OpenFlipper. *
11  *---------------------------------------------------------------------------*
12  * *
13  * Redistribution and use in source and binary forms, with or without *
14  * modification, are permitted provided that the following conditions *
15  * are met: *
16  * *
17  * 1. Redistributions of source code must retain the above copyright notice, *
18  * this list of conditions and the following disclaimer. *
19  * *
20  * 2. Redistributions in binary form must reproduce the above copyright *
21  * notice, this list of conditions and the following disclaimer in the *
22  * documentation and/or other materials provided with the distribution. *
23  * *
24  * 3. Neither the name of the copyright holder nor the names of its *
25  * contributors may be used to endorse or promote products derived from *
26  * this software without specific prior written permission. *
27  * *
28  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
29  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
30  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A *
31  * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
32  * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
33  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
34  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
35  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
36  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
37  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
38  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
39 * *
40 \*===========================================================================*/
41 
42 
43 
44 
45 
46 
47 //=============================================================================
48 //
49 // IMPLEMENTATION
50 //
51 //=============================================================================
52 
53 #define CURVATURE_C
54 
55 //== INCLUDES =================================================================
56 
57 #include <ACG/Geometry/Algorithms.hh>
58 #include "Math_Tools/Math_Tools.hh"
59 
60 #include <iostream>
61 #include <OpenMesh/Core/Geometry/MathDefs.hh>
62 
63 #include <cmath>
64 
65 //== NAMESPACES ===============================================================
66 
67 namespace curvature {
68 
69 //== IMPLEMENTATION ==========================================================
70 
73 template< typename MeshT >
74 double
75 gauss_curvature(MeshT& _mesh, const typename MeshT::VertexHandle& _vh) {
76  if (_mesh.status(_vh).deleted())
77  return 0.0;
78 
79  double gauss_curv = 2.0 * M_PI;
80 
81  /*
82 
83  TODO: Check the boundary case.
84 
85  If the vertex is a boundary vertex
86  if ( _mesh.is_boundary(_vh) )
87  gauss_curv = M_PI;
88 
89  */
90 
91  const typename MeshT::Point p0 = _mesh.point(_vh);
92 
93  typename MeshT::CVOHIter voh_it( _mesh.cvoh_iter(_vh));
94  typename MeshT::CVOHIter n_voh_it = voh_it;
95 
96  if ( ! voh_it->is_valid() )
97  return 0.0;
98 
99  // move to next
100  ++n_voh_it;
101 
102  for(; voh_it.is_valid(); ++voh_it, ++n_voh_it)
103  {
104  typename MeshT::Point p1 = _mesh.point(_mesh.to_vertex_handle( *voh_it));
105  typename MeshT::Point p2 = _mesh.point(_mesh.to_vertex_handle( *n_voh_it));
106 
107  gauss_curv -= acos(OpenMesh::sane_aarg( ((p1-p0).normalize() | (p2-p0).normalize()) ));
108  }
109 
110  return gauss_curv;
111 }
112 
113 
114 template<class MeshT, class VectorT, class REALT>
115 void discrete_mean_curv_op( const MeshT& _m,
116  const typename MeshT::VertexHandle& _vh,
117  VectorT& _n,
118  REALT& _area )
119 {
120  _n = VectorT(0,0,0);
121  _area = 0.0;
122 
123  typename MeshT::ConstVertexOHalfedgeIter voh_it = _m.cvoh_iter(_vh);
124 
125  if ( ! voh_it->is_valid() )
126  return;
127 
128  for(; voh_it.is_valid(); ++voh_it)
129  {
130  if ( _m.is_boundary( _m.edge_handle( *voh_it ) ) )
131  continue;
132 
133  const typename MeshT::Point p0 = _m.point( _vh );
134  const typename MeshT::Point p1 = _m.point( _m.to_vertex_handle( *voh_it));
135 // const typename MeshT::Point p2 = _m.point( _m.to_vertex_handle( _m.next_halfedge_handle( *voh_it)));
136  const typename MeshT::Point p2 = _m.point( _m.from_vertex_handle( _m.prev_halfedge_handle( *voh_it)));
137  const typename MeshT::Point p3 = _m.point( _m.to_vertex_handle( _m.next_halfedge_handle( _m.opposite_halfedge_handle(*voh_it))));
138 
139  const REALT alpha = acos( OpenMesh::sane_aarg((p0-p2).normalize() | (p1-p2).normalize()) );
140  const REALT beta = acos( OpenMesh::sane_aarg((p0-p3).normalize() | (p1-p3).normalize()) );
141 
142  REALT cotw = 0.0;
143 
144  if ( !OpenMesh::is_eq(alpha,M_PI/2.0) )
145  cotw += (REALT(1.0))/tan(alpha);
146 
147  if ( !OpenMesh::is_eq(beta,M_PI/2.0) )
148  cotw += (REALT(1.0))/tan(beta);
149 
150 #ifdef WIN32
151  if ( OpenMesh::is_zero(cotw) )
152  continue;
153 #else
154  if ( OpenMesh::is_zero(cotw) || std::isinf(cotw) )
155  continue;
156 #endif
157 
158  // calculate area
159  const int obt = ACG::Geometry::isObtuse(p0,p1,p2);
160  if(obt == 0)
161  {
162  REALT gamma = acos( OpenMesh::sane_aarg((p0-p1).normalize() | (p2-p1).normalize()) );
163 
164  REALT tmp = 0.0;
165  if ( !OpenMesh::is_eq(alpha,M_PI/2.0) )
166  tmp += (p0-p1).sqrnorm()*1.0/tan(alpha);
167 
168  if ( !OpenMesh::is_eq(gamma,M_PI/2.0) )
169  tmp += (p0-p2).sqrnorm()*1.0/tan(gamma);
170 
171 #ifdef WIN32
172  if ( OpenMesh::is_zero(tmp) )
173  continue;
174 #else
175  if ( OpenMesh::is_zero(tmp) || std::isinf(tmp) )
176  continue;
177 #endif
178 
179  _area += 0.125*( tmp );
180  }
181  else
182  {
183  if(obt == 1)
184  {
185  _area += ((p1-p0) % (p2-p0)).norm() * 0.5 * 0.5;
186  }
187  else
188  {
189  _area += ((p1-p0) % (p2-p0)).norm() * 0.5 * 0.25;
190  }
191  }
192 
193  _n += ((p0-p1)*cotw);
194 
195  // error handling
196  //if(_area < 0) std::cerr << "error: triangle area < 0\n";
197 // if(isnan(_area))
198 // {
199 // REALT gamma = acos( ((p0-p1).normalize() | (p2-p1).normalize()) );
200 
201 /*
202  std::cerr << "***************************\n";
203  std::cerr << "error : trianlge area = nan\n";
204  std::cerr << "alpha : " << alpha << std::endl;
205  std::cerr << "beta : " << beta << std::endl;
206  std::cerr << "gamma : " << gamma << std::endl;
207  std::cerr << "cotw : " << cotw << std::endl;
208  std::cerr << "normal: " << _n << std::endl;
209  std::cerr << "p0 : " << p0 << std::endl;
210  std::cerr << "p1 : " << p1 << std::endl;
211  std::cerr << "p2 : " << p2 << std::endl;
212  std::cerr << "p3 : " << p3 << std::endl;
213  std::cerr << "***************************\n";
214 */
215 // }
216  }
217 
218  _n /= 4.0*_area;
219 }
220 
221 
222 
223 //=============================================================================
224 } // curvature Namespace
225 //=============================================================================
T sane_aarg(T _aarg)
Trigonometry/angles - related.
Definition: MathDefs.hh:122
int isObtuse(const VectorT &_p0, const VectorT &_p1, const VectorT &_p2)
Definition: Algorithms.cc:1003
double gauss_curvature(MeshT &_mesh, const typename MeshT::VertexHandle &_vh)
void discrete_mean_curv_op(const MeshT &_m, const typename MeshT::VertexHandle &_vh, VectorT &_n, REALT &_area)
bool is_zero(const T &_a, Real _eps)
Definition: MathDefs.hh:61