Developer Documentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Curvature.cc
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 * $Revision$ *
45 * $LastChangedBy$ *
46 * $Date$ *
47 * *
48 \*===========================================================================*/
49 
50 
51 
52 
53 //=============================================================================
54 //
55 // IMPLEMENTATION
56 //
57 //=============================================================================
58 
59 #define CURVATURE_C
60 
61 //== INCLUDES =================================================================
62 
63 #include <ACG/Geometry/Algorithms.hh>
64 #include "Math_Tools/Math_Tools.hh"
65 
66 #include <iostream>
67 #include <OpenMesh/Core/Geometry/MathDefs.hh>
68 
69 #include <cmath>
70 
71 //== NAMESPACES ===============================================================
72 
73 namespace curvature {
74 
75 //== IMPLEMENTATION ==========================================================
76 
79 template< typename MeshT >
80 double
81 gauss_curvature(MeshT& _mesh, const typename MeshT::VertexHandle& _vh) {
82  if (_mesh.status(_vh).deleted())
83  return 0.0;
84 
85  double gauss_curv = 2.0 * M_PI;
86 
87  /*
88 
89  TODO: Check the boundary case.
90 
91  If the vertex is a boundary vertex
92  if ( _mesh.is_boundary(_vh) )
93  gauss_curv = M_PI;
94 
95  */
96 
97  const typename MeshT::Point p0 = _mesh.point(_vh);
98 
99  typename MeshT::CVOHIter voh_it( _mesh.cvoh_iter(_vh));
100  typename MeshT::CVOHIter n_voh_it = voh_it;
101 
102  if ( ! voh_it->is_valid() )
103  return 0.0;
104 
105  // move to next
106  ++n_voh_it;
107 
108  for(; voh_it.is_valid(); ++voh_it, ++n_voh_it)
109  {
110  typename MeshT::Point p1 = _mesh.point(_mesh.to_vertex_handle( *voh_it));
111  typename MeshT::Point p2 = _mesh.point(_mesh.to_vertex_handle( *n_voh_it));
112 
113  gauss_curv -= acos(OpenMesh::sane_aarg( ((p1-p0).normalize() | (p2-p0).normalize()) ));
114  }
115 
116  return gauss_curv;
117 }
118 
119 
120 template<class MeshT, class VectorT, class REALT>
121 void discrete_mean_curv_op( const MeshT& _m,
122  const typename MeshT::VertexHandle& _vh,
123  VectorT& _n,
124  REALT& _area )
125 {
126  _n = VectorT(0,0,0);
127  _area = 0.0;
128 
129  typename MeshT::ConstVertexOHalfedgeIter voh_it = _m.cvoh_iter(_vh);
130 
131  if ( ! voh_it->is_valid() )
132  return;
133 
134  for(; voh_it.is_valid(); ++voh_it)
135  {
136  if ( _m.is_boundary( _m.edge_handle( *voh_it ) ) )
137  continue;
138 
139  const typename MeshT::Point p0 = _m.point( _vh );
140  const typename MeshT::Point p1 = _m.point( _m.to_vertex_handle( *voh_it));
141 // const typename MeshT::Point p2 = _m.point( _m.to_vertex_handle( _m.next_halfedge_handle( *voh_it)));
142  const typename MeshT::Point p2 = _m.point( _m.from_vertex_handle( _m.prev_halfedge_handle( *voh_it)));
143  const typename MeshT::Point p3 = _m.point( _m.to_vertex_handle( _m.next_halfedge_handle( _m.opposite_halfedge_handle(*voh_it))));
144 
145  const REALT alpha = acos( OpenMesh::sane_aarg((p0-p2).normalize() | (p1-p2).normalize()) );
146  const REALT beta = acos( OpenMesh::sane_aarg((p0-p3).normalize() | (p1-p3).normalize()) );
147 
148  REALT cotw = 0.0;
149 
150  if ( !OpenMesh::is_eq(alpha,M_PI/2.0) )
151  cotw += (REALT(1.0))/tan(alpha);
152 
153  if ( !OpenMesh::is_eq(beta,M_PI/2.0) )
154  cotw += (REALT(1.0))/tan(beta);
155 
156 #ifdef WIN32
157  if ( OpenMesh::is_zero(cotw) )
158  continue;
159 #else
160  if ( OpenMesh::is_zero(cotw) || std::isinf(cotw) )
161  continue;
162 #endif
163 
164  // calculate area
165  const int obt = ACG::Geometry::isObtuse(p0,p1,p2);
166  if(obt == 0)
167  {
168  REALT gamma = acos( OpenMesh::sane_aarg((p0-p1).normalize() | (p2-p1).normalize()) );
169 
170  REALT tmp = 0.0;
171  if ( !OpenMesh::is_eq(alpha,M_PI/2.0) )
172  tmp += (p0-p1).sqrnorm()*1.0/tan(alpha);
173 
174  if ( !OpenMesh::is_eq(gamma,M_PI/2.0) )
175  tmp += (p0-p2).sqrnorm()*1.0/tan(gamma);
176 
177 #ifdef WIN32
178  if ( OpenMesh::is_zero(tmp) )
179  continue;
180 #else
181  if ( OpenMesh::is_zero(tmp) || std::isinf(tmp) )
182  continue;
183 #endif
184 
185  _area += 0.125*( tmp );
186  }
187  else
188  {
189  if(obt == 1)
190  {
191  _area += ((p1-p0) % (p2-p0)).norm() * 0.5 * 0.5;
192  }
193  else
194  {
195  _area += ((p1-p0) % (p2-p0)).norm() * 0.5 * 0.25;
196  }
197  }
198 
199  _n += ((p0-p1)*cotw);
200 
201  // error handling
202  //if(_area < 0) std::cerr << "error: triangle area < 0\n";
203 // if(isnan(_area))
204 // {
205 // REALT gamma = acos( ((p0-p1).normalize() | (p2-p1).normalize()) );
206 
207 /*
208  std::cerr << "***************************\n";
209  std::cerr << "error : trianlge area = nan\n";
210  std::cerr << "alpha : " << alpha << std::endl;
211  std::cerr << "beta : " << beta << std::endl;
212  std::cerr << "gamma : " << gamma << std::endl;
213  std::cerr << "cotw : " << cotw << std::endl;
214  std::cerr << "normal: " << _n << std::endl;
215  std::cerr << "p0 : " << p0 << std::endl;
216  std::cerr << "p1 : " << p1 << std::endl;
217  std::cerr << "p2 : " << p2 << std::endl;
218  std::cerr << "p3 : " << p3 << std::endl;
219  std::cerr << "***************************\n";
220 */
221 // }
222  }
223 
224  _n /= 4.0*_area;
225 }
226 
227 
228 
229 //=============================================================================
230 } // curvature Namespace
231 //=============================================================================
bool is_zero(const T &_a, Real _eps)
Definition: MathDefs.hh:66
int isObtuse(const VectorT &_p0, const VectorT &_p1, const VectorT &_p2)
Definition: Algorithms.cc:1013
T sane_aarg(T _aarg)
Trigonometry/angles - related.
Definition: MathDefs.hh:127