Developer Documentation
ModHausdorffT.cc
Go to the documentation of this file.
1 /* ========================================================================= *
2  * *
3  * OpenMesh *
4  * Copyright (c) 2001-2015, RWTH-Aachen University *
5  * Department of Computer Graphics and Multimedia *
6  * All rights reserved. *
7  * www.openmesh.org *
8  * *
9  *---------------------------------------------------------------------------*
10  * This file is part of OpenMesh. *
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  * $Date$ *
46  * *
47 \*===========================================================================*/
48 
53 //=============================================================================
54 //
55 // CLASS ModHausdorffT - IMPLEMENTATION
56 //
57 //=============================================================================
58 
59 #define OPENMESH_DECIMATER_MODHAUSDORFFT_C
60 
61 
62 //== INCLUDES =================================================================
63 
64 #include "ModHausdorffT.hh"
65 
66 
67 //== NAMESPACES ===============================================================
68 
69 namespace OpenMesh {
70 namespace Decimater {
71 
72 //== IMPLEMENTATION ==========================================================
73 
74 template <class MeshT>
75 typename ModHausdorffT<MeshT>::Scalar
77 distPointTriangleSquared( const Point& _p,
78  const Point& _v0,
79  const Point& _v1,
80  const Point& _v2 )
81 {
82  const Point v0v1 = _v1 - _v0;
83  const Point v0v2 = _v2 - _v0;
84  const Point n = v0v1 % v0v2; // not normalized !
85  const double d = n.sqrnorm();
86 
87 
88  // Check if the triangle is degenerated
89  if (d < FLT_MIN && d > -FLT_MIN) {
90  return -1.0;
91  }
92  const double invD = 1.0 / d;
93 
94  // these are not needed for every point, should still perform
95  // better with many points against one triangle
96  const Point v1v2 = _v2 - _v1;
97  const double inv_v0v2_2 = 1.0 / v0v2.sqrnorm();
98  const double inv_v0v1_2 = 1.0 / v0v1.sqrnorm();
99  const double inv_v1v2_2 = 1.0 / v1v2.sqrnorm();
100 
101 
102  Point v0p = _p - _v0;
103  Point t = v0p % n;
104  typename Point::value_type s01, s02, s12;
105  const double a = (t | v0v2) * -invD;
106  const double b = (t | v0v1) * invD;
107 
108  if (a < 0)
109  {
110  // Calculate the distance to an edge or a corner vertex
111  s02 = ( v0v2 | v0p ) * inv_v0v2_2;
112  if (s02 < 0.0)
113  {
114  s01 = ( v0v1 | v0p ) * inv_v0v1_2;
115  if (s01 <= 0.0) {
116  v0p = _v0;
117  } else if (s01 >= 1.0) {
118  v0p = _v1;
119  } else {
120  v0p = _v0 + v0v1 * s01;
121  }
122  } else if (s02 > 1.0) {
123  s12 = ( v1v2 | ( _p - _v1 )) * inv_v1v2_2;
124  if (s12 >= 1.0) {
125  v0p = _v2;
126  } else if (s12 <= 0.0) {
127  v0p = _v1;
128  } else {
129  v0p = _v1 + v1v2 * s12;
130  }
131  } else {
132  v0p = _v0 + v0v2 * s02;
133  }
134  } else if (b < 0.0) {
135  // Calculate the distance to an edge or a corner vertex
136  s01 = ( v0v1 | v0p ) * inv_v0v1_2;
137  if (s01 < 0.0)
138  {
139  // const Point n = v0v1 % v0v2; // not normalized !
140  s02 = ( v0v2 | v0p ) * inv_v0v2_2;
141  if (s02 <= 0.0) {
142  v0p = _v0;
143  } else if (s02 >= 1.0) {
144  v0p = _v2;
145  } else {
146  v0p = _v0 + v0v2 * s02;
147  }
148  } else if (s01 > 1.0) {
149  s12 = ( v1v2 | ( _p - _v1 )) * inv_v1v2_2;
150  if (s12 >= 1.0) {
151  v0p = _v2;
152  } else if (s12 <= 0.0) {
153  v0p = _v1;
154  } else {
155  v0p = _v1 + v1v2 * s12;
156  }
157  } else {
158  v0p = _v0 + v0v1 * s01;
159  }
160  } else if (a+b > 1.0) {
161  // Calculate the distance to an edge or a corner vertex
162  s12 = ( v1v2 | ( _p - _v1 )) * inv_v1v2_2;
163  if (s12 >= 1.0) {
164  s02 = ( v0v2 | v0p ) * inv_v0v2_2;
165  if (s02 <= 0.0) {
166  v0p = _v0;
167  } else if (s02 >= 1.0) {
168  v0p = _v2;
169  } else {
170  v0p = _v0 + v0v2*s02;
171  }
172  } else if (s12 <= 0.0) {
173  s01 = ( v0v1 | v0p ) * inv_v0v1_2;
174  if (s01 <= 0.0) {
175  v0p = _v0;
176  } else if (s01 >= 1.0) {
177  v0p = _v1;
178  } else {
179  v0p = _v0 + v0v1 * s01;
180  }
181  } else {
182  v0p = _v1 + v1v2 * s12;
183  }
184  } else {
185  // Calculate the distance to an interior point of the triangle
186  return ( (_p - n*((n|v0p) * invD)) - _p).sqrnorm();
187  }
188 
189  return (v0p - _p).sqrnorm();
190 }
191 
192 
193 template <class MeshT>
194 void
197 {
198  typename Mesh::FIter f_it(mesh_.faces_begin()), f_end(mesh_.faces_end());
199 
200  for (; f_it!=f_end; ++f_it)
201  mesh_.property(points_, *f_it).clear();
202 }
203 
204 
205 //-----------------------------------------------------------------------------
206 
207 
208 template <class MeshT>
209 float
212 {
213  std::vector<FaceHandle> faces; faces.reserve(20);
214  typename Mesh::VertexFaceIter vf_it;
215  typename Mesh::FaceHandle fh;
216  const typename Mesh::Scalar sqr_tolerace = tolerance_*tolerance_;
217  typename Mesh::CFVIter fv_it;
218  bool ok;
219 
220  // Clear the temporary point storage
221  tmp_points_.clear();
222 
223  // collect all points to be tested
224  // collect all faces to be tested against
225  for (vf_it=mesh_.vf_iter(_ci.v0); vf_it.is_valid(); ++vf_it) {
226  fh = *vf_it;
227 
228  if (fh != _ci.fl && fh != _ci.fr)
229  faces.push_back(fh);
230 
231  Points& pts = mesh_.property(points_, fh);
232  std::copy(pts.begin(), pts.end(), std::back_inserter(tmp_points_));
233  }
234 
235  // add point to be removed
236  tmp_points_.push_back(_ci.p0);
237 
238  // setup iterators
239  typename std::vector<FaceHandle>::iterator fh_it, fh_end(faces.end());
240  typename Points::const_iterator p_it, p_end(tmp_points_.end());
241 
242  // simulate collapse
243  mesh_.set_point(_ci.v0, _ci.p1);
244 
245  // for each point: try to find a face such that error is < tolerance
246  ok = true;
247 
248  for (p_it=tmp_points_.begin(); ok && p_it!=p_end; ++p_it) {
249  ok = false;
250 
251  for (fh_it=faces.begin(); !ok && fh_it!=fh_end; ++fh_it) {
252  fv_it=mesh_.cfv_iter(*fh_it);
253  const Point& p0 = mesh_.point(*fv_it);
254  const Point& p1 = mesh_.point(*(++fv_it));
255  const Point& p2 = mesh_.point(*(++fv_it));
256 
257  if ( distPointTriangleSquared(*p_it, p0, p1, p2) <= sqr_tolerace)
258  ok = true;
259  }
260  }
261 
262  // undo simulation changes
263  mesh_.set_point(_ci.v0, _ci.p0);
264 
265  return ( ok ? Base::LEGAL_COLLAPSE : Base::ILLEGAL_COLLAPSE );
266 }
267 
268 //-----------------------------------------------------------------------------
269 
270 template<class MeshT>
272  if (_factor >= 0.0 && _factor <= 1.0) {
273  // the smaller the factor, the smaller tolerance gets
274  // thus creating a stricter constraint
275  // division by error_tolerance_factor_ is for normalization
276  Scalar tolerance = tolerance_ * Scalar(_factor / this->error_tolerance_factor_);
277  set_tolerance(tolerance);
278  this->error_tolerance_factor_ = _factor;
279  }
280 }
281 
282 //-----------------------------------------------------------------------------
283 
284 template <class MeshT>
285 void
288 {
289  typename Mesh::VertexFaceIter vf_it;
290  FaceHandle fh;
291  std::vector<FaceHandle> faces;
292 
293 
294  // collect points & neighboring triangles
295 
296  tmp_points_.clear();
297  faces.reserve(20);
298 
299  // collect active faces and their points
300  for (vf_it=mesh_.vf_iter(_ci.v1); vf_it.is_valid(); ++vf_it) {
301  fh = *vf_it;
302  faces.push_back(fh);
303 
304  Points& pts = mesh_.property(points_, fh);
305  std::copy(pts.begin(), pts.end(), std::back_inserter(tmp_points_));
306  pts.clear();
307  }
308  if (faces.empty()) return; // should not happen anyway...
309 
310 
311  // collect points of the 2 deleted faces
312  if ((fh=_ci.fl).is_valid()) {
313  Points& pts = mesh_.property(points_, fh);
314  std::copy(pts.begin(), pts.end(), std::back_inserter(tmp_points_));
315  pts.clear();
316  }
317  if ((fh=_ci.fr).is_valid()) {
318  Points& pts = mesh_.property(points_, fh);
319  std::copy(pts.begin(), pts.end(), std::back_inserter(tmp_points_));
320  pts.clear();
321  }
322 
323  // add the deleted point
324  tmp_points_.push_back(_ci.p0);
325 
326 
327  // setup iterators
328  typename std::vector<FaceHandle>::iterator fh_it, fh_end(faces.end());
329  typename Points::const_iterator p_it, p_end(tmp_points_.end());
330 
331  // re-distribute points
332  Scalar emin, e;
333  typename Mesh::CFVIter fv_it;
334 
335  for (p_it=tmp_points_.begin(); p_it!=p_end; ++p_it) {
336  emin = FLT_MAX;
337 
338  for (fh_it=faces.begin(); fh_it!=fh_end; ++fh_it) {
339  fv_it=mesh_.cfv_iter(*fh_it);
340  const Point& p0 = mesh_.point(*fv_it);
341  const Point& p1 = mesh_.point(*(++fv_it));
342  const Point& p2 = mesh_.point(*(++fv_it));
343 
344  e = distPointTriangleSquared(*p_it, p0, p1, p2);
345  if (e < emin) {
346  emin = e;
347  fh = *fh_it;
348  }
349 
350  }
351 
352  mesh_.property(points_, fh).push_back(*p_it);
353  }
354 }
355 
356 
357 //-----------------------------------------------------------------------------
358 
359 
360 template <class MeshT>
361 typename ModHausdorffT<MeshT>::Scalar
363 compute_sqr_error(FaceHandle _fh, const Point& _p) const
364 {
365  typename Mesh::CFVIter fv_it = mesh_.cfv_iter(_fh);
366  const Point& p0 = mesh_.point(fv_it);
367  const Point& p1 = mesh_.point(++fv_it);
368  const Point& p2 = mesh_.point(++fv_it);
369 
370  const Points& points = mesh_.property(points_, _fh);
371  typename Points::const_iterator p_it = points.begin();
372  typename Points::const_iterator p_end = points.end();
373 
374  Point dummy;
375  Scalar e;
376  Scalar emax = distPointTriangleSquared(_p, p0, p1, p2);
377 
378 
379 
380  for (; p_it!=p_end; ++p_it) {
381  e = distPointTriangleSquared(*p_it, p0, p1, p2);
382  if (e > emax)
383  emax = e;
384  }
385 
386  return emax;
387 }
388 
389 
390 //=============================================================================
391 }
392 }
393 //=============================================================================
Mesh::VertexHandle v1
Remaining vertex.
Mesh::FaceHandle fl
Left face.
Mesh::FaceHandle fr
Right face.
Mesh::Point p1
Positions of remaining vertex.
Mesh::Point p0
Position of removed vertex.
virtual void postprocess_collapse(const CollapseInfo &_ci)
re-distribute points
Mesh::VertexHandle v0
Vertex to be removed.
virtual float collapse_priority(const CollapseInfo &_ci)
compute Hausdorff error for one-ring
Scalar compute_sqr_error(FaceHandle _fh, const Point &_p) const
compute max error for face _fh w.r.t. its point list and _p
void set_error_tolerance_factor(double _factor)
set the percentage of tolerance
Kernel::VertexFaceIter VertexFaceIter
Circulator.
Definition: PolyMeshT.hh:169
Scalar distPointTriangleSquared(const Point &_p, const Point &_v0, const Point &_v1, const Point &_v2)
squared distance from point _p to triangle (_v0, _v1, _v2)
Kernel::Scalar Scalar
Scalar type.
Definition: PolyMeshT.hh:113
virtual void initialize()
reset per-face point lists