Developer Documentation

DiffGeoT_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 // CLASS DiffGeoT - IMPLEMENTATION
46 //
47 //=============================================================================
48
49 #define DIFFGEOT_C
50
51 //== INCLUDES =================================================================
52
53 //#include <mb_base/mb_base.hh>
54 #include "DiffGeoT.hh"
55
56 #ifdef WIN32
57 #undef min
58 #undef max
59 #endif
60
61
62 //== NAMESPACES ===============================================================
63
64 namespace Remeshing {
65
66 //== IMPLEMENTATION ==========================================================
67
68
69 template <class Mesh>
70 DiffGeoT<Mesh>::
71 DiffGeoT(Mesh& _mesh)
72 : mesh_(_mesh),
73 weights_computed_(false),
74 area_computed_(false)
75 {
76 }
77
78
79 //-----------------------------------------------------------------------------
80
81
82 template <class Mesh>
83 DiffGeoT<Mesh>::
84 ~DiffGeoT()
85 {
86 mesh_.remove_property(area_);
87 mesh_.remove_property(gauss_curvature_);
88 mesh_.remove_property(mean_curvature_);
89 mesh_.remove_property(edge_weight_);
90 }
91
92
93 //-----------------------------------------------------------------------------
94
95
96 template <class Mesh>
97 void
98 DiffGeoT<Mesh>::
99 compute(unsigned int _post_smoothing_iters)
100 {
101 compute_edge_weights();
102 compute_area();
103 compute_gauss_curvature();
104 compute_mean_curvature();
105 post_smoothing(_post_smoothing_iters);
106 }
107
108
109 //-----------------------------------------------------------------------------
110
111
112 template <class Mesh>
113 void
114 DiffGeoT<Mesh>::
115 compute_edge_weights()
116 {
117 if (!edge_weight_.is_valid())
118 mesh_.add_property(edge_weight_);
119
120
121 typename Mesh::EdgeIter e_it, e_end(mesh_.edges_end());
122 typename Mesh::HalfedgeHandle heh2;
123 typename Mesh::Scalar weight;
124
125
126 for (e_it=mesh_.edges_begin(); e_it!=e_end; ++e_it)
127 {
128 const typename Mesh::HalfedgeHandle heh0 = mesh_.halfedge_handle(*e_it, 0);
129 const typename Mesh::Point& p0 = mesh_.point(mesh_.to_vertex_handle(heh0));
130
131 const typename Mesh::HalfedgeHandle heh1 = mesh_.halfedge_handle(*e_it, 1);
132 const typename Mesh::Point& p1 = mesh_.point(mesh_.to_vertex_handle(heh1));
133
134 heh2 = mesh_.next_halfedge_handle(heh0);
135 const typename Mesh::Point& p2 = mesh_.point(mesh_.to_vertex_handle(heh2));
136 const typename Mesh::Point d0 = (p0 - p2).normalize();
137 const typename Mesh::Point d1 = (p1 - p2).normalize();
138 weight = 1.0 / tan(acos(d0|d1));
139
140 heh2 = mesh_.next_halfedge_handle(heh1);
141 const typename Mesh::Point& p3 = mesh_.point(mesh_.to_vertex_handle(heh2));
142 const typename Mesh::Point d2 = (p0 - p3).normalize();
143 const typename Mesh::Point d3 = (p1 - p3).normalize();
144 weight += 1.0 / tan(acos(d2|d3));
145
146 mesh_.property(edge_weight_, *e_it) = weight;
147 }
148
149
150 weights_computed_ = true;
151 }
152
153
154 //-----------------------------------------------------------------------------
155
156
157 template <class Mesh>
158 void
159 DiffGeoT<Mesh>::
160 compute_area()
161 {
162 if (!area_.is_valid())
163 mesh_.add_property(area_);
164
165
166 typename Mesh::VertexIter v_it, v_end(mesh_.vertices_end());
167
168
169 for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
170 mesh_.property(area_, *v_it) = compute_area(*v_it);
171
172
173 area_computed_ = true;
174 }
175
176
177 //-----------------------------------------------------------------------------
178
179
180 template <class Mesh>
181 typename Mesh::Scalar
182 DiffGeoT<Mesh>::
183 compute_area(VertexHandle _vh) const
184 {
185 typename Mesh::HalfedgeHandle heh0, heh1, heh2;
186 typename Mesh::VertexOHalfedgeIter voh_it;
187
188 typename Mesh::Point P, Q, R, PQ, QR, PR;
189 double normPQ, normQR, normPR;
190 double angleP, angleQ, angleR;
191 double area;
192 double ub(0.999), lb(-0.999);
193 const double epsilon( 1e-12 );
194
195 area = 0.0;
196
197 for (voh_it=mesh_.voh_iter(_vh); voh_it.is_valid(); ++voh_it)
198 {
199 heh0 = *voh_it;
200 heh1 = mesh_.next_halfedge_handle(heh0);
201 heh2 = mesh_.next_halfedge_handle(heh1);
202
203 if (mesh_.is_boundary(heh0)) continue;
204
205 P = mesh_.point(mesh_.to_vertex_handle(heh2));
206 Q = mesh_.point(mesh_.to_vertex_handle(heh0));
207 R = mesh_.point(mesh_.to_vertex_handle(heh1));
208
209 (PQ = Q) -= P;
210 (QR = R) -= Q;
211 (PR = R) -= P;
212
213 normPQ = PQ.norm();
214 normQR = QR.norm();
215 normPR = PR.norm();
216
217 if( normPQ <= epsilon || normQR <= epsilon || normPR <= epsilon )
218 continue;
219
220 angleP = (PQ | PR) / normPQ / normPR;
221 angleQ = -(PQ | QR) / normPQ / normQR;
222 angleR = (QR | PR) / normQR / normPR;
223
224 if (angleP > ub) angleP = ub; else if (angleP < lb) angleP = lb;
225 if (angleQ > ub) angleQ = ub; else if (angleQ < lb) angleQ = lb;
226 if (angleR > ub) angleR = ub; else if (angleR < lb) angleR = lb;
227
228 angleP = acos(angleP);
229 angleQ = acos(angleQ);
230 angleR = acos(angleR);
231
232 if (angleP >= M_PI_2 || angleQ >= M_PI_2 || angleR >= M_PI_2)
233 {
234 if (angleP >= M_PI_2)
235 area += 0.25 * (PQ % PR).norm();
236 else
237 area += 0.125 * (PQ % PR).norm();
238 }
239 else
240 {
241 area += 0.125 * (normPR * normPR / tan(angleQ) +
242 normPQ * normPQ / tan(angleR));
243 }
244 }
245
246 return area;
247 }
248
249
250 //-----------------------------------------------------------------------------
251
252
253 template <class Mesh>
254 void
255 DiffGeoT<Mesh>::
256 compute_gauss_curvature()
257 {
258 if (!gauss_curvature_.is_valid())
259 mesh_.add_property(gauss_curvature_);
260
261 if (!area_computed_)
262 compute_area();
263
264
265 typename Mesh::VertexIter v_it, v_end(mesh_.vertices_end());
266 typename Mesh::VertexVertexIter vv_it, vv_it2;
267 typename Mesh::Point d0, d1;
268 typename Mesh::Scalar curv, count;
269 typename Mesh::Scalar angles, cos_angle;
270 typename Mesh::Scalar lb(-1.0), ub(1.0);
271
272
273 // compute for all non-boundary vertices
274 for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
275 {
276 if (!mesh_.is_boundary(*v_it))
277 {
278 angles = 0.0;
279
280 for (vv_it=mesh_.vv_iter(*v_it); vv_it.is_valid(); ++vv_it)
281 {
282 vv_it2 = vv_it; ++vv_it2;
283
284 d0 = (mesh_.point(*vv_it) - mesh_.point(*v_it)).normalize();
285 d1 = (mesh_.point(*vv_it2) - mesh_.point(*v_it)).normalize();
286
287 cos_angle = std::max(lb, std::min(ub, (d0 | d1)));
288 angles += acos(cos_angle);
289 }
290
291 mesh_.property(gauss_curvature_, *v_it) = (2*M_PI-angles) / mesh_.property(area_, *v_it);
292 }
293 }
294
295
296 // boundary vertices: get from neighbors
297 for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
298 {
299 if (mesh_.is_boundary(*v_it))
300 {
301 curv = count = 0.0;
302
303 for (vv_it=mesh_.vv_iter(*v_it); vv_it.is_valid(); ++vv_it)
304 if (!mesh_.is_boundary(*vv_it))
305 {
306 curv += mesh_.property(gauss_curvature_, *vv_it);
307 ++count;
308 }
309
310 if (count)
311 mesh_.property(gauss_curvature_, *v_it) = curv / count;
312 }
313 }
314 }
315
316
317 //-----------------------------------------------------------------------------
318
319
320 template <class Mesh>
321 void
322 DiffGeoT<Mesh>::
323 compute_mean_curvature()
324 {
325 if (!mean_curvature_.is_valid())
326 mesh_.add_property(mean_curvature_);
327
328 if (!weights_computed_)
329 compute_edge_weights();
330
331 if (!area_computed_)
332 compute_area();
333
334
335 typename Mesh::VertexIter v_it, v_end(mesh_.vertices_end());
336 typename Mesh::VertexVertexIter vv_it;
337 typename Mesh::VertexOHalfedgeIter voh_it;
338 typename Mesh::Scalar weight;
339 typename Mesh::Point umbrella;
340 typename Mesh::EdgeHandle eh;
341 typename Mesh::Scalar curv, count;
342
343
344 // compute for all non-boundary vertices
345 for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
346 {
347 if (!mesh_.is_boundary(*v_it))
348 {
349 umbrella[0] = umbrella[1] = umbrella[2] = 0.0;
350
351 for (voh_it=mesh_.voh_iter(*v_it); voh_it.is_valid(); ++voh_it)
352 {
353 eh = mesh_.edge_handle(*voh_it);
354 weight = mesh_.property(edge_weight_, eh);
355 umbrella += (mesh_.point(*v_it) - mesh_.point(mesh_.to_vertex_handle(*voh_it))) * weight;
356 }
357
358 mesh_.property(mean_curvature_, *v_it) =
359 umbrella.norm() / (4.0 * mesh_.property(area_, *v_it));
360 }
361 }
362
363
364 // boundary vertices: get from neighbors
365 for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
366 {
367 if (mesh_.is_boundary(*v_it))
368 {
369 curv = count = 0.0;
370
371 for (vv_it=mesh_.vv_iter(*v_it); vv_it.is_valid(); ++vv_it)
372 if (!mesh_.is_boundary(*vv_it))
373 {
374 curv += mesh_.property(mean_curvature_, *vv_it);
375 ++count;
376 }
377
378 if (count)
379 mesh_.property(mean_curvature_, *v_it) = curv / count;
380 }
381 }
382 }
383
384
385 //-----------------------------------------------------------------------------
386
387
388 template <class Mesh>
389 void
390 DiffGeoT<Mesh>::
391 post_smoothing(unsigned int _iters)
392 {
393 // early out
394 if (!_iters) return;
395
396
397 // something should already be there...
398 if (! (weights_computed_ &&
399 area_computed_ &&
400 mean_curvature_.is_valid() &&
401 gauss_curvature_.is_valid()) )
402 {
403 std::cerr << "DiffGeoT::post_smoothing: something is wrong!\n";
404 return;
405 }
406
407
408
409 typename Mesh::VertexIter v_it, v_end(mesh_.vertices_end());
410 typename Mesh::VertexOHalfedgeIter voh_it;
411 typename Mesh::Scalar w, ww;
412 typename Mesh::Scalar gc, mc;
413 typename Mesh::EdgeHandle eh;
414
415
416
417 // add helper properties to store new values
418 OpenMesh::VPropHandleT<Scalar> new_gauss, new_mean;
419 mesh_.add_property(new_gauss);
420 mesh_.add_property(new_mean);
421
422
423
424 for (unsigned int i=0; i<_iters; ++i)
425 {
426
427 // compute new value
428 for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
429 {
430 if (!mesh_.is_boundary(*v_it))
431 {
432 gc = mc = ww = 0.0;
433
434 for (voh_it=mesh_.voh_iter(*v_it); voh_it.is_valid(); ++voh_it)
435 {
436 eh = mesh_.edge_handle(*voh_it);
437 ww += (w = mesh_.property(edge_weight_, eh));
438 mc += w * mesh_.property(mean_curvature_, mesh_.to_vertex_handle(*voh_it));
439 gc += w * mesh_.property(gauss_curvature_, mesh_.to_vertex_handle(*voh_it));
440 }
441
442 if (ww)
443 {
444 mesh_.property(new_mean, *v_it) = mc / ww;
445 mesh_.property(new_gauss, *v_it) = gc / ww;
446 }
447 }
448 }
449
450
451 // replace old by new value
452 for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
453 {
454 if (!mesh_.is_boundary(*v_it))
455 {
456 mesh_.property(mean_curvature_, *v_it) = mesh_.property(new_mean, *v_it);
457
458 mesh_.property(gauss_curvature_, *v_it) = mesh_.property(new_gauss, *v_it);
459 }
460 }
461 }
462
463
464 // remove helper properties
465 mesh_.remove_property(new_gauss);
466 mesh_.remove_property(new_mean);
467 }
468
469
470 //=============================================================================
471 } // namespace Remeshing
472 //=============================================================================
OpenMesh::PolyMeshT::VertexVertexIter
Kernel::VertexVertexIter VertexVertexIter
Circulator.
Definition: PolyMeshT.hh:162
OpenMesh::VPropHandleT< Scalar >
OpenMesh::PolyMeshT::Point
Kernel::Point Point
Coordinate type.
Definition: PolyMeshT.hh:112
OpenMesh::TriMesh_ArrayKernelT
Definition: PolyMesh_ArrayKernelT.hh:73
OpenMesh::BaseHandle::is_valid
bool is_valid() const
The handle is valid iff the index is not negative.
Definition: Handles.hh:72
OpenMesh::PolyMeshT::VertexOHalfedgeIter
Kernel::VertexOHalfedgeIter VertexOHalfedgeIter
Circulator.
Definition: PolyMeshT.hh:163
OpenMesh::PolyMeshT::Scalar
Kernel::Scalar Scalar
Scalar type.
Definition: PolyMeshT.hh:110
Remeshing
Definition: AdaptiveRemesherT.hh:64