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::VertexOHalfedgeIter
Kernel::VertexOHalfedgeIter VertexOHalfedgeIter
Circulator.
Definition: PolyMeshT.hh:163
OpenMesh::PolyMeshT::Point
Kernel::Point Point
Coordinate type.
Definition: PolyMeshT.hh:112
OpenMesh::PolyMeshT::Scalar
Kernel::Scalar Scalar
Scalar type.
Definition: PolyMeshT.hh:110
OpenMesh::PolyMeshT::VertexVertexIter
Kernel::VertexVertexIter VertexVertexIter
Circulator.
Definition: PolyMeshT.hh:162
OpenMesh::BaseHandle::is_valid
bool is_valid() const
The handle is valid iff the index is not negative.
Definition: Handles.hh:72