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