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