Developer Documentation

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

- File List
- File Members

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 //=============================================================================
VectorT
OpenMesh::is_zero
bool is_zero(const T &_a, Real _eps)
Definition: MathDefs.hh:66
ACG::Geometry::isObtuse
int isObtuse(const VectorT &_p0, const VectorT &_p1, const VectorT &_p2)
Definition: Algorithms.cc:1013
curvature
Definition: Curvature.cc:73
OpenMesh::sane_aarg
T sane_aarg(T _aarg)
Trigonometry/angles - related.
Definition: MathDefs.hh:127