Developer Documentation
vdpmanalyzer.cc
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 // -------------------------------------------------------------- includes ----
45 
46 #include <OpenMesh/Core/System/config.h>
47 // -------------------- STL
48 #include <iostream>
49 #include <fstream>
50 #include <algorithm>
51 #include <limits>
52 #include <exception>
53 #include <cmath>
54 // -------------------- OpenMesh
55 #include <OpenMesh/Core/IO/MeshIO.hh>
56 #include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
57 #include <OpenMesh/Core/Utils/Endian.hh>
58 #include <OpenMesh/Tools/Utils/Timer.hh>
59 #include <OpenMesh/Tools/Utils/getopt.h>
60 // -------------------- OpenMesh VDPM
61 #include <OpenMesh/Tools/VDPM/StreamingDef.hh>
62 #include <OpenMesh/Tools/VDPM/ViewingParameters.hh>
63 #include <OpenMesh/Tools/VDPM/VHierarchy.hh>
64 #include <OpenMesh/Tools/VDPM/VFront.hh>
65 
66 // ----------------------------------------------------------------------------
67 
68 
69 // ----------------------------------------------------------------------------
70 
71 using namespace OpenMesh;
72 
73 // ----------------------------------------------------------------------------
74 
75 // using view dependent progressive mesh
76 using VDPM::Plane3d;
77 using VDPM::VFront;
78 using VDPM::VHierarchy;
79 using VDPM::VHierarchyNode;
80 using VDPM::VHierarchyNodeIndex;
81 using VDPM::VHierarchyNodeHandle;
82 using VDPM::VHierarchyNodeHandleContainer;
83 using VDPM::ViewingParameters;
84 
85 
86 // ------------------------------------------------------------- mesh type ----
87 // Generate an OpenMesh suitable for the analyzing task
88 
89 struct AnalyzerTraits : public DefaultTraits
90 {
91  VertexTraits
92  {
93  public:
94 
95  VHierarchyNodeHandle vhierarchy_node_handle()
96  { return node_handle_; }
97 
98  void set_vhierarchy_node_handle(VHierarchyNodeHandle _node_handle)
99  { node_handle_ = _node_handle; }
100 
101  bool is_ancestor(const VHierarchyNodeIndex &_other)
102  { return false; }
103 
104  private:
105 
106  VHierarchyNodeHandle node_handle_;
107 
108  };
109 
110 
111  HalfedgeTraits
112  {
113  public:
114 
115  VHierarchyNodeHandle
116  vhierarchy_leaf_node_handle()
117  { return leaf_node_handle_; }
118 
119  void
120  set_vhierarchy_leaf_node_handle(VHierarchyNodeHandle _leaf_node_handle)
121  { leaf_node_handle_ = _leaf_node_handle; }
122 
123  private:
124 
125  VHierarchyNodeHandle leaf_node_handle_;
126 
127  };
128 
129  VertexAttributes(Attributes::Status |
130  Attributes::Normal);
131  HalfedgeAttributes(Attributes::PrevHalfedge);
132  EdgeAttributes(Attributes::Status);
133  FaceAttributes(Attributes::Status |
134  Attributes::Normal);
135 
136 };
137 
138 typedef TriMesh_ArrayKernelT<AnalyzerTraits> Mesh;
139 
140 // ----------------------------------------------------------------- types ----
141 
142 struct PMInfo
143 {
144  Mesh::Point p0;
145  Mesh::VertexHandle v0, v1, vl, vr;
146 };
147 
148 typedef std::vector<PMInfo> PMInfoContainer;
149 typedef PMInfoContainer::iterator PMInfoIter;
150 typedef std::vector<VertexHandle> VertexHandleContainer;
151 typedef std::vector<Vec3f> ResidualContainer;
152 
153 
154 // -------------------------------------------------------------- forwards ----
155 
157 void open_prog_mesh(const std::string &_filename);
158 
160 void save_vd_prog_mesh(const std::string &_filename);
161 
162 
164 void locate_fund_cut_vertices();
165 
166 void create_vertex_hierarchy();
167 
168 
170 void refine(unsigned int _n);
171 
173 void coarsen(unsigned int _n);
174 
175 void vdpm_analysis();
176 
177 void get_leaf_node_handles(VHierarchyNodeHandle node_handle,
178  VHierarchyNodeHandleContainer &leaf_nodes);
179 void compute_bounding_box(VHierarchyNodeHandle node_handle,
180  VHierarchyNodeHandleContainer &leaf_nodes);
181 void compute_cone_of_normals(VHierarchyNodeHandle node_handle,
182  VHierarchyNodeHandleContainer &leaf_nodes);
183 void compute_screen_space_error(VHierarchyNodeHandle node_handle,
184  VHierarchyNodeHandleContainer &leaf_nodes);
185 void compute_mue_sigma(VHierarchyNodeHandle node_handle,
186  ResidualContainer &residuals);
187 
188 Vec3f
189 point2triangle_residual(const Vec3f &p,
190  const Vec3f tri[3], float &s, float &t);
191 
192 
193 void PrintOutFundCuts();
194 void PrintVertexNormals();
195 
196 // --------------------------------------------------------------- globals ----
197 // mesh data
198 
199 Mesh mesh_;
200 PMInfoContainer pminfos_;
201 PMInfoIter pmiter_;
202 
203 VHierarchy vhierarchy_;
204 
205 unsigned int n_base_vertices_, n_base_faces_, n_details_;
206 unsigned int n_current_res_;
207 unsigned int n_max_res_;
208 bool verbose = false;
209 
210 
211 // ----------------------------------------------------------------------------
212 
213 void usage_and_exit(int xcode)
214 {
215  using namespace std;
216 
217  cout << "Usage: vdpmanalyzer [-h] [-o output.spm] input.pm\n";
218 
219  exit(xcode);
220 }
221 
222 // ----------------------------------------------------------------------------
223 
224 inline
225 std::string&
226 replace_extension( std::string& _s, const std::string& _e )
227 {
228  std::string::size_type dot = _s.rfind(".");
229  if (dot == std::string::npos)
230  { _s += "." + _e; }
231  else
232  { _s = _s.substr(0,dot+1)+_e; }
233  return _s;
234 }
235 
236 // ----------------------------------------------------------------------------
237 
238 inline
239 std::string
240 basename(const std::string& _f)
241 {
242  std::string::size_type dot = _f.rfind("/");
243  if (dot == std::string::npos)
244  return _f;
245  return _f.substr(dot+1, _f.length()-(dot+1));
246 }
247 
248 
249 // ----------------------------------------------------------------------------
250 // just for debugging
251 
252 typedef std::vector<OpenMesh::Vec3f> MyPoints;
253 typedef MyPoints::iterator MyPointsIter;
254 
255 MyPoints projected_points;
256 MyPoints original_points;
257 
258 
259 // ------------------------------------------------------------------ main ----
260 
261 
262 int main(int argc, char **argv)
263 {
264  int c;
265  std::string ifname;
266  std::string ofname;
267 
268  while ( (c=getopt(argc, argv, "o:"))!=-1 )
269  {
270  switch(c)
271  {
272  case 'v': verbose = true; break;
273  case 'o': ofname = optarg; break;
274  case 'h': usage_and_exit(0); break;
275  default: usage_and_exit(1);
276  }
277  }
278 
279  if (optind >= argc)
280  usage_and_exit(1);
281 
282  ifname = argv[optind];
283 
284  if (ofname == "." || ofname == ".." )
285  ofname += "/" + basename(ifname);
286  std::string spmfname = ofname.empty() ? ifname : ofname;
287  replace_extension(spmfname, "spm");
288 
289  if ( ifname.empty() || spmfname.empty() )
290  {
291  usage_and_exit(1);
292  }
293 
294  try
295  {
296  open_prog_mesh(ifname);
297  vdpm_analysis();
298  save_vd_prog_mesh(spmfname);
299  }
300  catch( std::bad_alloc& )
301  {
302  std::cerr << "Error: out of memory!\n" << std::endl;
303  return 1;
304  }
305  catch( std::exception& x )
306  {
307  std::cerr << "Error: " << x.what() << std::endl;
308  return 1;
309  }
310  catch( ... )
311  {
312  std::cerr << "Fatal! Unknown error!\n";
313  return 1;
314  }
315  return 0;
316 }
317 
318 
319 // ----------------------------------------------------------------------------
320 
321 void
322 open_prog_mesh(const std::string& _filename)
323 {
324  Mesh::Point p;
325  unsigned int i, i0, i1, i2;
326  unsigned int v1, vl, vr;
327  char c[10];
328  VertexHandle vertex_handle;
329  VHierarchyNodeHandle node_handle, lchild_handle, rchild_handle;
330  VHierarchyNodeIndex node_index;
331 
332  std::ifstream ifs(_filename.c_str(), std::ios::binary);
333  if (!ifs)
334  {
335  std::cerr << "read error\n";
336  exit(1);
337  }
338 
339  //
340  bool swap = Endian::local() != Endian::LSB;
341 
342  // read header
343  ifs.read(c, 8); c[8] = '\0';
344  if (std::string(c) != std::string("ProgMesh"))
345  {
346  std::cerr << "Wrong file format.\n";
347  ifs.close();
348  exit(1);
349  }
350  IO::restore(ifs, n_base_vertices_, swap);
351  IO::restore(ifs, n_base_faces_, swap);
352  IO::restore(ifs, n_details_, swap);
353 
354  vhierarchy_.set_num_roots(n_base_vertices_);
355 
356  for (i=0; i<n_base_vertices_; ++i)
357  {
358  IO::restore(ifs, p, swap);
359 
360  vertex_handle = mesh_.add_vertex(p);
361  node_index = vhierarchy_.generate_node_index(i, 1);
362  node_handle = vhierarchy_.add_node();
363 
364  vhierarchy_.node(node_handle).set_index(node_index);
365  vhierarchy_.node(node_handle).set_vertex_handle(vertex_handle);
366  mesh_.data(vertex_handle).set_vhierarchy_node_handle(node_handle);
367  }
368 
369  for (i=0; i<n_base_faces_; ++i)
370  {
371  IO::restore(ifs, i0, swap);
372  IO::restore(ifs, i1, swap);
373  IO::restore(ifs, i2, swap);
374  mesh_.add_face(mesh_.vertex_handle(i0),
375  mesh_.vertex_handle(i1),
376  mesh_.vertex_handle(i2));
377  }
378 
379  // load progressive detail
380  for (i=0; i<n_details_; ++i)
381  {
382  IO::restore(ifs, p, swap);
383  IO::restore(ifs, v1, swap);
384  IO::restore(ifs, vl, swap);
385  IO::restore(ifs, vr, swap);
386 
387  PMInfo pminfo;
388  pminfo.p0 = p;
389  pminfo.v0 = mesh_.add_vertex(p);
390  pminfo.v1 = Mesh::VertexHandle(v1);
391  pminfo.vl = Mesh::VertexHandle(vl);
392  pminfo.vr = Mesh::VertexHandle(vr);
393  pminfos_.push_back(pminfo);
394 
395  node_handle = mesh_.data(pminfo.v1).vhierarchy_node_handle();
396 
397  vhierarchy_.make_children(node_handle);
398  lchild_handle = vhierarchy_.lchild_handle(node_handle);
399  rchild_handle = vhierarchy_.rchild_handle(node_handle);
400 
401  mesh_.data(pminfo.v0).set_vhierarchy_node_handle(lchild_handle);
402  mesh_.data(pminfo.v1).set_vhierarchy_node_handle(rchild_handle);
403  vhierarchy_.node(lchild_handle).set_vertex_handle(pminfo.v0);
404  vhierarchy_.node(rchild_handle).set_vertex_handle(pminfo.v1);
405  }
406 
407  ifs.close();
408 
409 
410  // recover mapping between basemesh vertices to roots of vertex hierarchy
411  for (i=0; i<n_base_vertices_; ++i)
412  {
413  node_handle = vhierarchy_.root_handle(i);
414  vertex_handle = vhierarchy_.node(node_handle).vertex_handle();
415 
416  mesh_.data(vertex_handle).set_vhierarchy_node_handle(node_handle);
417  }
418 
419  pmiter_ = pminfos_.begin();
420  n_current_res_ = 0;
421  n_max_res_ = n_details_;
422 
423 
424  // update face and vertex normals
425  mesh_.update_face_normals();
426  mesh_.update_vertex_normals();
427 
428  // bounding box
429  Mesh::ConstVertexIter
430  vIt(mesh_.vertices_begin()),
431  vEnd(mesh_.vertices_end());
432 
433  Mesh::Point bbMin, bbMax;
434 
435  bbMin = bbMax = mesh_.point(*vIt);
436  for (; vIt!=vEnd; ++vIt)
437  {
438  bbMin.minimize(mesh_.point(*vIt));
439  bbMax.maximize(mesh_.point(*vIt));
440  }
441 
442  // info
443  std::cerr << mesh_.n_vertices() << " vertices, "
444  << mesh_.n_edges() << " edge, "
445  << mesh_.n_faces() << " faces, "
446  << n_details_ << " detail vertices\n";
447 }
448 
449 
450 // ----------------------------------------------------------------------------
451 
452 void
453 save_vd_prog_mesh(const std::string &_filename)
454 {
455  unsigned i;
456  int fvi[3];
457  Mesh::Point p;
458  Vec3f normal;
459  float radius, sin_square, mue_square, sigma_square;
460  Mesh::FaceIter f_it;
461  Mesh::HalfedgeHandle hh;
462  Mesh::VertexHandle vh;
463  VHierarchyNodeIndex node_index;
464  VHierarchyNodeIndex fund_lcut_index, fund_rcut_index;
465  VHierarchyNodeHandle node_handle, lchild_handle, rchild_handle;
466  std::map<VertexHandle, unsigned int> handle2index_map;
467 
468  std::ofstream ofs(_filename.c_str(), std::ios::binary);
469  if (!ofs)
470  {
471  std::cerr << "write error\n";
472  exit(1);
473  }
474 
475  //
476  bool swap = Endian::local() != Endian::LSB;
477 
478  // write header
479  ofs << "VDProgMesh";
480 
481  IO::store(ofs, n_base_vertices_, swap);
482  IO::store(ofs, n_base_faces_, swap);
483  IO::store(ofs, n_details_, swap);
484 
485  // write base mesh
486  coarsen(0);
487  mesh_.garbage_collection( false, true, true );
488 
489  for (i=0; i<n_base_vertices_; ++i)
490  {
491  node_handle = vhierarchy_.root_handle(i);
492  vh = vhierarchy_.node(node_handle).vertex_handle();
493 
494  p = mesh_.point(vh);
495  radius = vhierarchy_.node(node_handle).radius();
496  normal = vhierarchy_.node(node_handle).normal();
497  sin_square = vhierarchy_.node(node_handle).sin_square();
498  mue_square = vhierarchy_.node(node_handle).mue_square();
499  sigma_square = vhierarchy_.node(node_handle).sigma_square();
500 
501  IO::store(ofs, p, swap);
502  IO::store(ofs, radius, swap);
503  IO::store(ofs, normal, swap);
504  IO::store(ofs, sin_square, swap);
505  IO::store(ofs, mue_square, swap);
506  IO::store(ofs, sigma_square, swap);
507 
508  handle2index_map[vh] = i;
509  }
510 
511 
512  for (f_it=mesh_.faces_begin(); f_it!=mesh_.faces_end(); ++f_it) {
513  hh = mesh_.halfedge_handle(*f_it);
514  vh = mesh_.to_vertex_handle(hh);
515  fvi[0] = handle2index_map[vh];
516 
517  hh = mesh_.next_halfedge_handle(hh);
518  vh = mesh_.to_vertex_handle(hh);
519  fvi[1] = handle2index_map[vh];
520 
521  hh = mesh_.next_halfedge_handle(hh);
522  vh = mesh_.to_vertex_handle(hh);
523  fvi[2] = handle2index_map[vh];
524 
525  IO::store(ofs, fvi[0], swap);
526  IO::store(ofs, fvi[1], swap);
527  IO::store(ofs, fvi[2], swap);
528  }
529 
530 
531  // write progressive detail (vertex hierarchy)
532 
533  for (i=0; i<n_details_; ++i)
534  {
535  PMInfo pminfo = *pmiter_;
536 
537  p = mesh_.point(pminfo.v0);
538 
539  IO::store(ofs, p, swap);
540 
541 
542  node_handle = mesh_.data(pminfo.v1).vhierarchy_node_handle();
543  lchild_handle = vhierarchy_.lchild_handle(node_handle);
544  rchild_handle = vhierarchy_.rchild_handle(node_handle);
545 
546  node_index = vhierarchy_.node(node_handle).node_index();
547  fund_lcut_index = vhierarchy_.node(node_handle).fund_lcut_index();
548  fund_rcut_index = vhierarchy_.node(node_handle).fund_rcut_index();
549 
550  IO::store(ofs, node_index.value(), swap);
551  IO::store(ofs, fund_lcut_index.value(), swap);
552  IO::store(ofs, fund_rcut_index.value(), swap);
553 
554  radius = vhierarchy_.node(lchild_handle).radius();
555  normal = vhierarchy_.node(lchild_handle).normal();
556  sin_square = vhierarchy_.node(lchild_handle).sin_square();
557  mue_square = vhierarchy_.node(lchild_handle).mue_square();
558  sigma_square = vhierarchy_.node(lchild_handle).sigma_square();
559 
560  IO::store(ofs, radius, swap);
561  IO::store(ofs, normal, swap);
562  IO::store(ofs, sin_square, swap);
563  IO::store(ofs, mue_square, swap);
564  IO::store(ofs, sigma_square, swap);
565 
566  radius = vhierarchy_.node(rchild_handle).radius();
567  normal = vhierarchy_.node(rchild_handle).normal();
568  sin_square = vhierarchy_.node(rchild_handle).sin_square();
569  mue_square = vhierarchy_.node(rchild_handle).mue_square();
570  sigma_square = vhierarchy_.node(rchild_handle).sigma_square();
571 
572  IO::store(ofs, radius, swap);
573  IO::store(ofs, normal, swap);
574  IO::store(ofs, sin_square, swap);
575  IO::store(ofs, mue_square, swap);
576  IO::store(ofs, sigma_square, swap);
577 
578  refine(i);
579  }
580 
581  ofs.close();
582 
583  std::cout << "save view-dependent progressive mesh" << std::endl;
584 }
585 
586 //-----------------------------------------------------------------------------
587 
588 void refine(unsigned int _n)
589 {
590  while (n_current_res_ < _n && pmiter_ != pminfos_.end())
591  {
592  mesh_.vertex_split(pmiter_->v0,
593  pmiter_->v1,
594  pmiter_->vl,
595  pmiter_->vr);
596 
597  VHierarchyNodeHandle
598  parent_handle = mesh_.data(pmiter_->v1).vhierarchy_node_handle();
599 
600  VHierarchyNodeHandle
601  lchild_handle = vhierarchy_.lchild_handle(parent_handle),
602  rchild_handle = vhierarchy_.rchild_handle(parent_handle);
603 
604  mesh_.data(pmiter_->v0).set_vhierarchy_node_handle(lchild_handle);
605  mesh_.data(pmiter_->v1).set_vhierarchy_node_handle(rchild_handle);
606 
607  ++pmiter_;
608  ++n_current_res_;
609  }
610 }
611 
612 
613 //-----------------------------------------------------------------------------
614 
615 
616 void coarsen(unsigned int _n)
617 {
618  while (n_current_res_ > _n && pmiter_ != pminfos_.begin())
619  {
620  --pmiter_;
621 
622  Mesh::HalfedgeHandle hh =
623  mesh_.find_halfedge(pmiter_->v0, pmiter_->v1);
624  mesh_.collapse(hh);
625 
626  VHierarchyNodeHandle
627  rchild_handle = mesh_.data(pmiter_->v1).vhierarchy_node_handle();
628 
629  VHierarchyNodeHandle
630  parent_handle = vhierarchy_.parent_handle(rchild_handle);
631 
632  mesh_.data(pmiter_->v1).set_vhierarchy_node_handle(parent_handle);
633 
634  --n_current_res_;
635  }
636 }
637 
638 
639 //-----------------------------------------------------------------------------
640 
641 
642 void
643 vdpm_analysis()
644 {
645  unsigned int i;
646  Mesh::VertexHandle vh;
647  Mesh::VertexIter v_it;
648  Mesh::HalfedgeIter h_it;
649  Mesh::HalfedgeHandle h, o, hn, op, hpo, on, ono;
650  VHierarchyNodeHandleContainer leaf_nodes;
651 
652  OpenMesh::Utils::Timer tana;
653  tana.start();
654 
655  refine(n_max_res_);
656 
657  mesh_.update_face_normals();
658  mesh_.update_vertex_normals();
659 
660  std::cout << "Init view-dependent PM analysis" << std::endl;
661 
662  // initialize
663  for (h_it=mesh_.halfedges_begin(); h_it!=mesh_.halfedges_end(); ++h_it)
664  {
665  vh = mesh_.to_vertex_handle(*h_it);
666  mesh_.data(*h_it).set_vhierarchy_leaf_node_handle(mesh_.data(vh).vhierarchy_node_handle());
667  }
668 
669  for (v_it=mesh_.vertices_begin(); v_it!=mesh_.vertices_end(); ++v_it)
670  {
671  VHierarchyNodeHandle
672  node_handle = mesh_.data(*v_it).vhierarchy_node_handle();
673 
674  vhierarchy_.node(node_handle).set_normal(mesh_.normal(*v_it));
675  }
676 
677  std::cout << "Start view-dependent PM analysis" << std::endl;
678 
679  // locate fundamental cut vertices in each edge collapse
680  OpenMesh::Utils::Timer t;
681 
682  for (i=n_max_res_; i>0; --i)
683  {
684  t.start();
685  PMInfo pminfo = pminfos_[i-1];
686 
687  if (verbose)
688  std::cout << "Analyzing " << i << "-th detail vertex" << std::endl;
689 
690  // maintain leaf node pointers & locate fundamental cut vertices
691  h = mesh_.find_halfedge(pminfo.v0, pminfo.v1);
692  o = mesh_.opposite_halfedge_handle(h);
693  hn = mesh_.next_halfedge_handle(h);
694  hpo = mesh_.opposite_halfedge_handle(mesh_.prev_halfedge_handle(h));
695  op = mesh_.prev_halfedge_handle(o);
696  on = mesh_.next_halfedge_handle(o);
697  ono = mesh_.opposite_halfedge_handle(on);
698 
699  VHierarchyNodeHandle
700  rchild_handle = mesh_.data(pminfo.v1).vhierarchy_node_handle();
701 
702  VHierarchyNodeHandle
703  parent_handle = vhierarchy_.parent_handle(rchild_handle);
704 
705  if (pminfo.vl != Mesh::InvalidVertexHandle)
706  {
707  VHierarchyNodeHandle
708  fund_lcut_handle = mesh_.data(hn).vhierarchy_leaf_node_handle();
709 
710  VHierarchyNodeHandle
711  left_leaf_handle = mesh_.data(hpo).vhierarchy_leaf_node_handle();
712 
713  mesh_.data(hn).set_vhierarchy_leaf_node_handle(left_leaf_handle);
714 
715  vhierarchy_.node(parent_handle).
716  set_fund_lcut(vhierarchy_.node_index(fund_lcut_handle));
717  }
718 
719  if (pminfo.vr != Mesh::InvalidVertexHandle)
720  {
721  VHierarchyNodeHandle
722  fund_rcut_handle = mesh_.data(on).vhierarchy_leaf_node_handle(),
723  right_leaf_handle = mesh_.data(ono).vhierarchy_leaf_node_handle();
724 
725  mesh_.data(op).set_vhierarchy_leaf_node_handle(right_leaf_handle);
726 
727  vhierarchy_.node(parent_handle).
728  set_fund_rcut(vhierarchy_.node_index(fund_rcut_handle));
729  }
730 
731  coarsen(i-1);
732 
733  leaf_nodes.clear();
734 
735  get_leaf_node_handles(parent_handle, leaf_nodes);
736  compute_bounding_box(parent_handle, leaf_nodes);
737  compute_cone_of_normals(parent_handle, leaf_nodes);
738  compute_screen_space_error(parent_handle, leaf_nodes);
739 
740  t.stop();
741 
742  if (verbose)
743  {
744  std::cout << " radius of bounding sphere: "
745  << vhierarchy_.node(parent_handle).radius() << std::endl;
746  std::cout << " direction of cone of normals: "
747  << vhierarchy_.node(parent_handle).normal() << std::endl;
748  std::cout << " sin(semi-angle of cone of normals) ^2: "
749  << vhierarchy_.node(parent_handle).sin_square() << std::endl;
750  std::cout << " (mue^2, sigma^2) : ("
751  << vhierarchy_.node(parent_handle).mue_square() << ", "
752  << vhierarchy_.node(parent_handle).sigma_square() << ")"
753  << std::endl;
754  std::cout << "- " << t.as_string() << std::endl;
755  }
756 
757  } // end for all collapses
758 
759  tana.stop();
760  std::cout << "Analyzing step completed in "
761  << tana.as_string() << std::endl;
762 }
763 
764 
765 // ----------------------------------------------------------------------------
766 
767 void
768 get_leaf_node_handles(VHierarchyNodeHandle node_handle,
769  VHierarchyNodeHandleContainer &leaf_nodes)
770 {
771  if (vhierarchy_.node(node_handle).is_leaf())
772  {
773  leaf_nodes.push_back(node_handle);
774  }
775  else
776  {
777  get_leaf_node_handles(vhierarchy_.node(node_handle).lchild_handle(),
778  leaf_nodes);
779  get_leaf_node_handles(vhierarchy_.node(node_handle).rchild_handle(),
780  leaf_nodes);
781  }
782 }
783 
784 
785 // ----------------------------------------------------------------------------
786 
787 void
788 compute_bounding_box(VHierarchyNodeHandle node_handle, VHierarchyNodeHandleContainer &leaf_nodes)
789 {
790  float max_distance;
791  Vec3f p, lp;
792  VHierarchyNodeHandleContainer::iterator n_it, n_end(leaf_nodes.end());
793 
794  max_distance = 0.0f;
795  VertexHandle vh = vhierarchy_.node(node_handle).vertex_handle();
796  p = mesh_.point(vh);
797  for ( n_it = leaf_nodes.begin(); n_it != n_end; ++n_it )
798  {
799  lp = mesh_.point(vhierarchy_.vertex_handle(*n_it));
800  max_distance = std::max(max_distance, (p - lp).length());
801  }
802 
803  vhierarchy_.node(node_handle).set_radius(max_distance);
804 }
805 
806 
807 // ----------------------------------------------------------------------------
808 
809 void
810 compute_cone_of_normals(VHierarchyNodeHandle node_handle,
811  VHierarchyNodeHandleContainer &leaf_nodes)
812 {
813  Vec3f n, ln;
814  VertexHandle vh = vhierarchy_.node(node_handle).vertex_handle();
815  VHierarchyNodeHandleContainer::iterator n_it, n_end(leaf_nodes.end());
816 
817  n = mesh_.calc_vertex_normal(vh);
818  float max_angle = 0.0f;
819 
820  n_it = leaf_nodes.begin();
821  while( n_it != n_end )
822  {
823  ln = vhierarchy_.node(*n_it).normal();
824  const float angle = acosf( dot(n,ln) );
825  max_angle = std::max(max_angle, angle );
826 
827  ++n_it;
828  }
829 
830  max_angle = std::min(max_angle, float(M_PI_2));
831  mesh_.set_normal(vh, n);
832  vhierarchy_.node(node_handle).set_normal(n);
833  vhierarchy_.node(node_handle).set_semi_angle(max_angle);
834 }
835 
836 
837 // ----------------------------------------------------------------------------
838 
839 void
840 compute_screen_space_error(VHierarchyNodeHandle node_handle, VHierarchyNodeHandleContainer &leaf_nodes)
841 {
842  std::vector<Vec3f> residuals;
843  Mesh::VertexFaceIter vf_it;
844  Mesh::HalfedgeHandle heh;
845  Mesh::VertexHandle vh;
846  Vec3f residual;
847  Vec3f res;
848  Vec3f lp;
849 #if ((defined(_MSC_VER) && (_MSC_VER >= 1800)) )
850  // Workaround for internal compiler error
851  Vec3f tri[3]{ {},{},{} };
852 #else
853  Vec3f tri[3];
854 #endif
855  float s, t;
856  VHierarchyNodeHandleContainer::iterator n_it, n_end(leaf_nodes.end());
857 
858  for ( n_it = leaf_nodes.begin(); n_it != n_end; ++n_it )
859  {
860  lp = mesh_.point(vhierarchy_.node(*n_it).vertex_handle());
861 
862  // compute residual of a leaf-vertex from the current mesh_
863  vh = vhierarchy_.node(node_handle).vertex_handle();
864  residual = lp - mesh_.point(vh);
865  float min_distance = residual.length();
866 
867  for (vf_it=mesh_.vf_iter(vh); vf_it.is_valid(); ++vf_it)
868  {
869  heh = mesh_.halfedge_handle(*vf_it);
870  tri[0] = mesh_.point(mesh_.to_vertex_handle(heh));
871  heh = mesh_.next_halfedge_handle(heh);
872  tri[1] = mesh_.point(mesh_.to_vertex_handle(heh));
873  heh = mesh_.next_halfedge_handle(heh);
874  tri[2] = mesh_.point(mesh_.to_vertex_handle(heh));
875 
876  res = point2triangle_residual(lp, tri, s, t);
877 
878  if (res.length() < min_distance)
879  {
880  residual = res;
881  min_distance = res.length();
882  }
883  }
884 
885  residuals.push_back(residual);
886  }
887 
888  compute_mue_sigma(node_handle, residuals);
889 }
890 
891 
892 // ----------------------------------------------------------------------------
893 
894 void
895 compute_mue_sigma(VHierarchyNodeHandle node_handle,
896  ResidualContainer &residuals)
897 {
898  Vec3f vn;
899  float max_inner, max_cross;
900  ResidualContainer::iterator r_it, r_end(residuals.end());
901 
902  max_inner = max_cross = 0.0f;
903  vn = mesh_.normal(vhierarchy_.node(node_handle).vertex_handle());
904  for (r_it = residuals.begin(); r_it != r_end; ++r_it)
905  {
906  float inner = fabsf(dot(*r_it, vn));
907  float cross = OpenMesh::cross(*r_it, vn).length();
908 
909  max_inner = std::max(max_inner, inner);
910  max_cross = std::max(max_cross, cross);
911  }
912 
913  if (max_cross < 1.0e-7)
914  {
915  vhierarchy_.node(node_handle).set_mue(max_cross);
916  vhierarchy_.node(node_handle).set_sigma(max_inner);
917  }
918  else {
919  float ratio = std::max(1.0f, max_inner/max_cross);
920  float whole_degree = acosf(1.0f/ratio);
921  float mue, max_mue;
922  Vec3f res;
923 
924  max_mue = 0.0f;
925  for (r_it = residuals.begin(); r_it != r_end; ++r_it)
926  {
927  res = *r_it;
928  float res_length = res.length();
929 
930  // TODO: take care when res.length() is too small
931  float degree = acosf(dot(vn,res) / res_length);
932 
933  if (degree < 0.0f) degree = -degree;
934  if (degree > float(M_PI_2)) degree = float(M_PI) - degree;
935 
936  if (degree < whole_degree)
937  mue = cosf(whole_degree - degree) * res_length;
938  else
939  mue = res_length;
940 
941  max_mue = std::max(max_mue, mue);
942  }
943 
944  vhierarchy_.node(node_handle).set_mue(max_mue);
945  vhierarchy_.node(node_handle).set_sigma(ratio*max_mue);
946  }
947 }
948 
949 // ----------------------------------------------------------------------------
950 
951 
952 Vec3f
953 point2triangle_residual(const Vec3f &p, const Vec3f tri[3], float &s, float &t)
954 {
955  OpenMesh::Vec3f B = tri[0]; // Tri.Origin();
956  OpenMesh::Vec3f E0 = tri[1] - tri[0]; // rkTri.Edge0()
957  OpenMesh::Vec3f E1 = tri[2] - tri[0]; // rkTri.Edge1()
958  OpenMesh::Vec3f D = tri[0] - p; // kDiff
959  float a = dot(E0, E0); // fA00
960  float b = dot(E0, E1); // fA01
961  float c = dot(E1, E1); // fA11
962  float d = dot(E0, D); // fB0
963  float e = dot(E1, D); // fB1
964  //float f = dot(D, D); // fC
965  float det = fabsf(a*c - b*b);
966  s = b*e-c*d;
967  t = b*d-a*e;
968 
969  OpenMesh::Vec3f residual;
970 
971 // float distance2;
972 
973  if ( s + t <= det )
974  {
975  if ( s < 0.0f )
976  {
977  if ( t < 0.0f ) // region 4
978  {
979  if ( d < 0.0f )
980  {
981  t = 0.0f;
982  if ( -d >= a )
983  {
984  s = 1.0f;
985 // distance2 = a+2.0f*d+f;
986  }
987  else
988  {
989  s = -d/a;
990 // distance2 = d*s+f;
991  }
992  }
993  else
994  {
995  s = 0.0f;
996  if ( e >= 0.0f )
997  {
998  t = 0.0f;
999 // distance2 = f;
1000  }
1001  else if ( -e >= c )
1002  {
1003  t = 1.0f;
1004 // distance2 = c+2.0f*e+f;
1005  }
1006  else
1007  {
1008  t = -e/c;
1009 // distance2 = e*t+f;
1010  }
1011  }
1012  }
1013  else // region 3
1014  {
1015  s = 0.0f;
1016  if ( e >= 0.0f )
1017  {
1018  t = 0.0f;
1019 // distance2 = f;
1020  }
1021  else if ( -e >= c )
1022  {
1023  t = 1.0f;
1024 // distance2 = c+2.0f*e+f;
1025  }
1026  else
1027  {
1028  t = -e/c;
1029 // distance2 = e*t+f;
1030  }
1031  }
1032  }
1033  else if ( t < 0.0f ) // region 5
1034  {
1035  t = 0.0f;
1036  if ( d >= 0.0f )
1037  {
1038  s = 0.0f;
1039 // distance2 = f;
1040  }
1041  else if ( -d >= a )
1042  {
1043  s = 1.0f;
1044 // distance2 = a+2.0f*d+f;
1045  }
1046  else
1047  {
1048  s = -d/a;
1049 // distance2 = d*s+f;
1050  }
1051  }
1052  else // region 0
1053  {
1054  // minimum at interior point
1055  float inv_det = 1.0f/det;
1056  s *= inv_det;
1057  t *= inv_det;
1058 // distance2 = s*(a*s+b*t+2.0f*d) + t*(b*s+c*t+2.0f*e)+f;
1059  }
1060  }
1061  else
1062  {
1063  float tmp0, tmp1, numer, denom;
1064 
1065  if ( s < 0.0f ) // region 2
1066  {
1067  tmp0 = b + d;
1068  tmp1 = c + e;
1069  if ( tmp1 > tmp0 )
1070  {
1071  numer = tmp1 - tmp0;
1072  denom = a-2.0f*b+c;
1073  if ( numer >= denom )
1074  {
1075  s = 1.0f;
1076  t = 0.0f;
1077 // distance2 = a+2.0f*d+f;
1078  }
1079  else
1080  {
1081  s = numer/denom;
1082  t = 1.0f - s;
1083 // distance2 = s*(a*s+b*t+2.0f*d) + t*(b*s+c*t+2.0f*e)+f;
1084  }
1085  }
1086  else
1087  {
1088  s = 0.0f;
1089  if ( tmp1 <= 0.0f )
1090  {
1091  t = 1.0f;
1092 // distance2 = c+2.0f*e+f;
1093  }
1094  else if ( e >= 0.0f )
1095  {
1096  t = 0.0f;
1097 // distance2 = f;
1098  }
1099  else
1100  {
1101  t = -e/c;
1102 // distance2 = e*t+f;
1103  }
1104  }
1105  }
1106  else if ( t < 0.0f ) // region 6
1107  {
1108  tmp0 = b + e;
1109  tmp1 = a + d;
1110  if ( tmp1 > tmp0 )
1111  {
1112  numer = tmp1 - tmp0;
1113  denom = a-2.0f*b+c;
1114  if ( numer >= denom )
1115  {
1116  t = 1.0f;
1117  s = 0.0f;
1118 // distance2 = c+2.0f*e+f;
1119  }
1120  else
1121  {
1122  t = numer/denom;
1123  s = 1.0f - t;
1124 // distance2 = s*(a*s+b*t+2.0f*d)+ t*(b*s+c*t+2.0f*e)+f;
1125  }
1126  }
1127  else
1128  {
1129  t = 0.0f;
1130  if ( tmp1 <= 0.0f )
1131  {
1132  s = 1.0f;
1133 // distance2 = a+2.0f*d+f;
1134  }
1135  else if ( d >= 0.0f )
1136  {
1137  s = 0.0f;
1138 // distance2 = f;
1139  }
1140  else
1141  {
1142  s = -d/a;
1143 // distance2 = d*s+f;
1144  }
1145  }
1146  }
1147  else // region 1
1148  {
1149  numer = c + e - b - d;
1150  if ( numer <= 0.0f )
1151  {
1152  s = 0.0f;
1153  t = 1.0f;
1154 // distance2 = c+2.0f*e+f;
1155  }
1156  else
1157  {
1158  denom = a-2.0f*b+c;
1159  if ( numer >= denom )
1160  {
1161  s = 1.0f;
1162  t = 0.0f;
1163 // distance2 = a+2.0f*d+f;
1164  }
1165  else
1166  {
1167  s = numer/denom;
1168  t = 1.0f - s;
1169 // distance2 = s*(a*s+b*t+2.0f*d) + t*(b*s+c*t+2.0f*e)+f;
1170  }
1171  }
1172  }
1173  }
1174 
1175  residual = p - (B + s*E0 + t*E1);
1176 
1177  return residual;
1178 }
1179 
1180 // ============================================================================
OpenMesh::VDPM::VHierarchyNode::rchild_handle
VHierarchyNodeHandle rchild_handle()
Returns handle to right child.
Definition: VHierarchyNode.hh:111
OpenMesh::cross
osg::Vec3f cross(const osg::Vec3f &_v1, const osg::Vec3f &_v2)
Adapter for osg vector member computing a scalar product.
Definition: VectorAdapter.hh:191
OpenMesh::PolyMeshT::VertexFaceIter
Kernel::VertexFaceIter VertexFaceIter
Circulator.
Definition: PolyMeshT.hh:166
OpenMesh::Attributes::PrevHalfedge
Add storage for previous halfedge (halfedges). The bit is set by default in the DefaultTraits.
Definition: Attributes.hh:84
OpenMesh::TriMesh_ArrayKernelT
Definition: PolyMesh_ArrayKernelT.hh:73
std
STL namespace.
OpenMesh::angle
T angle(T _cos_angle, T _sin_angle)
Definition: MathDefs.hh:140
OpenMesh::Utils::Timer::as_string
std::string as_string(Format format=Automatic)
OpenMesh::VertexHandle
Handle for a vertex entity.
Definition: Handles.hh:120
OpenMesh::Attributes::Status
Add status to mesh item (all items)
Definition: Attributes.hh:85
OpenMesh::TriMeshT::vertex_split
HalfedgeHandle vertex_split(Point _v0_point, VertexHandle _v1, VertexHandle _vl, VertexHandle _vr)
Vertex Split: inverse operation to collapse().
Definition: TriMeshT.hh:214
VertexTraits
#define VertexTraits
Macro for defining the vertex traits. See Specifying your MyMesh.
Definition: Traits.hh:91
OpenMesh::Attributes::Normal
Add normals to mesh item (vertices/faces)
Definition: Attributes.hh:82
OpenMesh::Utils::Timer::stop
void stop(void)
Stop measurement.
VertexAttributes
#define VertexAttributes(_i)
Macro for defining the vertex attributes. See Specifying your MyMesh.
Definition: Traits.hh:79
OpenMesh::Endian::LSB
Little endian (Intel family and clones)
Definition: Endian.hh:78
OpenMesh::VDPM::VHierarchyNode::lchild_handle
VHierarchyNodeHandle lchild_handle()
Returns handle to left child.
Definition: VHierarchyNode.hh:108
EdgeAttributes
#define EdgeAttributes(_i)
Macro for defining the edge attributes. See Specifying your MyMesh.
Definition: Traits.hh:85
OpenMesh::PolyMeshT::update_vertex_normals
void update_vertex_normals()
Update normal vectors for all vertices.
Definition: PolyMeshT_impl.hh:550
OpenMesh::PolyMeshT::add_vertex
SmartVertexHandle add_vertex(const Point &_p)
Alias for new_vertex(const Point&).
Definition: PolyMeshT.hh:235
OpenMesh::VDPM::VHierarchyNode::is_leaf
bool is_leaf() const
Returns true, if node is leaf else false.
Definition: VHierarchyNode.hh:101
OpenMesh::PolyMeshT::update_face_normals
void update_face_normals()
Update normal vectors for all faces.
Definition: PolyMeshT_impl.hh:335
OpenMesh::VectorT::length
auto length() const -> decltype(std::declval< VectorT< S, DIM >>().norm())
compute squared euclidean norm
Definition: Vector11T.hh:423
HalfedgeTraits
#define HalfedgeTraits
Macro for defining the halfedge traits. See Specifying your MyMesh.
Definition: Traits.hh:95
OpenMesh::Utils::Timer::start
void start(void)
Start measurement.
HalfedgeAttributes
#define HalfedgeAttributes(_i)
Macro for defining the halfedge attributes. See Specifying your MyMesh.
Definition: Traits.hh:82
OpenMesh::VDPM::VHierarchyNodeHandleContainer
std::vector< VHierarchyNodeHandle > VHierarchyNodeHandleContainer
Container for vertex hierarchy node handles.
Definition: VHierarchyNode.hh:183
FaceAttributes
#define FaceAttributes(_i)
Macro for defining the face attributes. See Specifying your MyMesh.
Definition: Traits.hh:88
OpenMesh::dot
osg::Vec3f::ValueType dot(const osg::Vec3f &_v1, const osg::Vec3f &_v2)
Adapter for osg vector member computing a scalar product.
Definition: VectorAdapter.hh:176
OpenMesh::Endian::local
static Type local()
Return endian type of host system.
Definition: Endian.hh:83
OpenMesh::PolyMeshT::calc_vertex_normal
Normal calc_vertex_normal(VertexHandle _vh) const
Calculate vertex normal for one specific vertex.
Definition: PolyMeshT_impl.hh:461