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::TriMesh_ArrayKernelT
Definition: PolyMesh_ArrayKernelT.hh:73
OpenMesh::VDPM::VHierarchyNode::lchild_handle
VHierarchyNodeHandle lchild_handle()
Returns handle to left child.
Definition: VHierarchyNode.hh:108
OpenMesh::Endian::LSB
Little endian (Intel family and clones)
Definition: Endian.hh:78
OpenMesh::VertexHandle
Handle for a vertex entity.
Definition: Handles.hh:120
OpenMesh::VDPM::VHierarchyNode::rchild_handle
VHierarchyNodeHandle rchild_handle()
Returns handle to right child.
Definition: VHierarchyNode.hh:111
OpenMesh::PolyMeshT::Point
Kernel::Point Point
Coordinate type.
Definition: PolyMeshT.hh:112
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::add_vertex
VertexHandle add_vertex(const Point &_p)
Alias for new_vertex(const Point&).
Definition: PolyMeshT.hh:235
OpenMesh::VDPM::VHierarchy
Definition: VHierarchy.hh:73
OpenMesh::Utils::Timer
Definition: Timer.hh:82
OpenMesh::Utils::Timer::as_string
std::string as_string(Format format=Automatic)
AnalyzerTraits
Definition: vdpmanalyzer.cc:89
OpenMesh::DefaultTraits
Definition: Traits.hh:121
OpenMesh::PolyMeshT::VertexFaceIter
Kernel::VertexFaceIter VertexFaceIter
Circulator.
Definition: PolyMeshT.hh:166
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
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::PolyMeshT::VertexHandle
Kernel::VertexHandle VertexHandle
Handle for referencing the corresponding item.
Definition: PolyMeshT.hh:136
OpenMesh::PolyMeshT::update_vertex_normals
void update_vertex_normals()
Update normal vectors for all vertices.
Definition: PolyMeshT_impl.hh:466
VertexTraits
#define VertexTraits
Macro for defining the vertex traits. See Specifying your MyMesh.
Definition: Traits.hh:91
config.h
OpenMesh::VDPM::VHierarchyNodeHandle
Definition: VHierarchyNode.hh:77
OpenMesh::angle
T angle(T _cos_angle, T _sin_angle)
Definition: MathDefs.hh:140
OpenMesh::Endian::local
static Type local()
Return endian type of host system.
Definition: Endian.hh:83
Timer.hh
HalfedgeAttributes
#define HalfedgeAttributes(_i)
Macro for defining the halfedge attributes. See Specifying your MyMesh.
Definition: Traits.hh:82
OpenMesh::VDPM::VHierarchyNode::is_leaf
bool is_leaf() const
Returns true, if node is leaf else false.
Definition: VHierarchyNode.hh:101
std
STL namespace.
OpenMesh::VDPM::Plane3d
Definition: Plane3d.hh:78
EdgeAttributes
#define EdgeAttributes(_i)
Macro for defining the edge attributes. See Specifying your MyMesh.
Definition: Traits.hh:85
OpenMesh::VectorT::length
auto length() const -> decltype(std::declval< VectorT< S, DIM >>().norm())
compute squared euclidean norm
Definition: Vector11T.hh:418
OpenMesh::VDPM::VHierarchyNodeIndex
Definition: VHierarchyNodeIndex.hh:73
OpenMesh::Attributes::Status
Add status to mesh item (all items)
Definition: Attributes.hh:85
OpenMesh::VDPM::VHierarchyNodeHandleContainer
std::vector< VHierarchyNodeHandle > VHierarchyNodeHandleContainer
Container for vertex hierarchy node handles.
Definition: VHierarchyNode.hh:183
OpenMesh::VectorT< float, 3 >
OpenMesh::PolyMeshT::update_face_normals
void update_face_normals()
Update normal vectors for all faces.
Definition: PolyMeshT_impl.hh:277
OpenMesh::VDPM::VFront
Definition: VFront.hh:75
OpenMesh::Attributes::PrevHalfedge
Add storage for previous halfedge (halfedges). The bit is set by default in the DefaultTraits.
Definition: Attributes.hh:84
OpenMesh::VDPM::ViewingParameters
Definition: ViewingParameters.hh:74
OpenMesh::Utils::Timer::start
void start(void)
Start measurement.
HalfedgeTraits
#define HalfedgeTraits
Macro for defining the halfedge traits. See Specifying your MyMesh.
Definition: Traits.hh:95
OpenMesh::VDPM::VHierarchyNode
Definition: VHierarchyNode.hh:90
OpenMesh::PolyMeshT::calc_vertex_normal
Normal calc_vertex_normal(VertexHandle _vh) const
Calculate vertex normal for one specific vertex.
Definition: PolyMeshT_impl.hh:391
OpenMesh::Attributes::Normal
Add normals to mesh item (vertices/faces)
Definition: Attributes.hh:82
OpenMesh
Definition: MeshItems.hh:59
VertexAttributes
#define VertexAttributes(_i)
Macro for defining the vertex attributes. See Specifying your MyMesh.
Definition: Traits.hh:79
OpenMesh::Utils::Timer::stop
void stop(void)
Stop measurement.
FaceAttributes
#define FaceAttributes(_i)
Macro for defining the face attributes. See Specifying your MyMesh.
Definition: Traits.hh:88