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