Developer Documentation

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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
Definition: MeshItems.hh:64
OpenMesh::VDPM::ViewingParameters
Definition: ViewingParameters.hh:79
OpenMesh::VDPM::VHierarchyNode
Definition: VHierarchyNode.hh:95
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::TriMesh_ArrayKernelT
Definition: TriMesh_ArrayKernelT.hh:98
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)
OpenMesh::VDPM::Plane3d
Definition: Plane3d.hh:83
VertexTraits
#define VertexTraits
Macro for defining the vertex traits. See Specifying your MyMesh.
Definition: Traits.hh:96
config.h
OpenMesh::VDPM::VHierarchy
Definition: VHierarchy.hh:78
OpenMesh::VDPM::VFront
Definition: VFront.hh:80
OpenMesh::Endian::LSB
Little endian (Intel family and clones)
Definition: Endian.hh:83
Timer.hh
OpenMesh::VectorT
Definition: Vector11T.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::VDPM::VHierarchyNodeIndex
Definition: VHierarchyNodeIndex.hh:78
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::DefaultTraits
Definition: Traits.hh:126
OpenMesh::VDPM::VHierarchyNodeHandle
Definition: VHierarchyNode.hh:82
OpenMesh::PolyMeshT::VertexFaceIter
Kernel::VertexFaceIter VertexFaceIter
Circulator.
Definition: PolyMeshT.hh:169
OpenMesh::Utils::Timer
Definition: Timer.hh:87
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
AnalyzerTraits
Definition: vdpmanalyzer.cc:94
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