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 #if defined(__FreeBSD__)
55 #include <unistd.h>
56 #endif
57 // -------------------- OpenMesh
58 #include <OpenMesh/Core/IO/MeshIO.hh>
59 #include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
60 #include <OpenMesh/Core/Utils/Endian.hh>
61 #include <OpenMesh/Tools/Utils/Timer.hh>
62 #include <OpenMesh/Tools/Utils/getopt.h>
63 // -------------------- OpenMesh VDPM
64 #include <OpenMesh/Tools/VDPM/StreamingDef.hh>
65 #include <OpenMesh/Tools/VDPM/ViewingParameters.hh>
66 #include <OpenMesh/Tools/VDPM/VHierarchy.hh>
67 #include <OpenMesh/Tools/VDPM/VFront.hh>
68
69 // ----------------------------------------------------------------------------
70
71
72 // ----------------------------------------------------------------------------
73
74 using namespace OpenMesh;
75
76 // ----------------------------------------------------------------------------
77
78 // using view dependent progressive mesh
79 using VDPM::Plane3d;
80 using VDPM::VFront;
81 using VDPM::VHierarchy;
82 using VDPM::VHierarchyNode;
83 using VDPM::VHierarchyNodeIndex;
84 using VDPM::VHierarchyNodeHandle;
85 using VDPM::VHierarchyNodeHandleContainer;
86 using VDPM::ViewingParameters;
87
88
89 // ------------------------------------------------------------- mesh type ----
90 // Generate an OpenMesh suitable for the analyzing task
91
92 struct AnalyzerTraits : public DefaultTraits
93 {
94 VertexTraits
95 {
96 public:
97
98 VHierarchyNodeHandle vhierarchy_node_handle()
99 { return node_handle_; }
100
101 void set_vhierarchy_node_handle(VHierarchyNodeHandle _node_handle)
102 { node_handle_ = _node_handle; }
103
104 bool is_ancestor(const VHierarchyNodeIndex &_other)
105 { return false; }
106
107 private:
108
109 VHierarchyNodeHandle node_handle_;
110
111 };
112
113
114 HalfedgeTraits
115 {
116 public:
117
118 VHierarchyNodeHandle
119 vhierarchy_leaf_node_handle()
120 { return leaf_node_handle_; }
121
122 void
123 set_vhierarchy_leaf_node_handle(VHierarchyNodeHandle _leaf_node_handle)
124 { leaf_node_handle_ = _leaf_node_handle; }
125
126 private:
127
128 VHierarchyNodeHandle leaf_node_handle_;
129
130 };
131
132 VertexAttributes(Attributes::Status |
133 Attributes::Normal);
134 HalfedgeAttributes(Attributes::PrevHalfedge);
135 EdgeAttributes(Attributes::Status);
136 FaceAttributes(Attributes::Status |
137 Attributes::Normal);
138
139 };
140
141 typedef TriMesh_ArrayKernelT<AnalyzerTraits> Mesh;
142
143 // ----------------------------------------------------------------- types ----
144
145 struct PMInfo
146 {
147 Mesh::Point p0;
148 Mesh::VertexHandle v0, v1, vl, vr;
149 };
150
151 typedef std::vector<PMInfo> PMInfoContainer;
152 typedef PMInfoContainer::iterator PMInfoIter;
153 typedef std::vector<VertexHandle> VertexHandleContainer;
154 typedef std::vector<Vec3f> ResidualContainer;
155
156
157 // -------------------------------------------------------------- forwards ----
158
159 /// open progressive mesh
160 void open_prog_mesh(const std::string &_filename);
161
162 /// save view-dependent progressive mesh
163 void save_vd_prog_mesh(const std::string &_filename);
164
165
166 /// locate fundamental cut vertices
167 void locate_fund_cut_vertices();
168
169 void create_vertex_hierarchy();
170
171
172 /// refine mesh up to _n vertices
173 void refine(unsigned int _n);
174
175 /// coarsen mesh down to _n vertices
176 void coarsen(unsigned int _n);
177
178 void vdpm_analysis();
179
180 void get_leaf_node_handles(VHierarchyNodeHandle node_handle,
181 VHierarchyNodeHandleContainer &leaf_nodes);
182 void compute_bounding_box(VHierarchyNodeHandle node_handle,
183 VHierarchyNodeHandleContainer &leaf_nodes);
184 void compute_cone_of_normals(VHierarchyNodeHandle node_handle,
185 VHierarchyNodeHandleContainer &leaf_nodes);
186 void compute_screen_space_error(VHierarchyNodeHandle node_handle,
187 VHierarchyNodeHandleContainer &leaf_nodes);
188 void compute_mue_sigma(VHierarchyNodeHandle node_handle,
189 ResidualContainer &residuals);
190
191 Vec3f
192 point2triangle_residual(const Vec3f &p,
193 const Vec3f tri[3], float &s, float &t);
194
195
196 void PrintOutFundCuts();
197 void PrintVertexNormals();
198
199 // --------------------------------------------------------------- globals ----
200 // mesh data
201
202 Mesh mesh_;
203 PMInfoContainer pminfos_;
204 PMInfoIter pmiter_;
205
206 VHierarchy vhierarchy_;
207
208 unsigned int n_base_vertices_, n_base_faces_, n_details_;
209 unsigned int n_current_res_;
210 unsigned int n_max_res_;
211 bool verbose = false;
212
213
214 // ----------------------------------------------------------------------------
215
usage_and_exit(int xcode)216 void usage_and_exit(int xcode)
217 {
218 using namespace std;
219
220 cout << "Usage: vdpmanalyzer [-h] [-o output.spm] input.pm\n";
221
222 exit(xcode);
223 }
224
225 // ----------------------------------------------------------------------------
226
227 inline
228 std::string&
replace_extension(std::string & _s,const std::string & _e)229 replace_extension( std::string& _s, const std::string& _e )
230 {
231 std::string::size_type dot = _s.rfind(".");
232 if (dot == std::string::npos)
233 { _s += "." + _e; }
234 else
235 { _s = _s.substr(0,dot+1)+_e; }
236 return _s;
237 }
238
239 // ----------------------------------------------------------------------------
240
241 inline
242 std::string
basename(const std::string & _f)243 basename(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
255 typedef std::vector<OpenMesh::Vec3f> MyPoints;
256 typedef MyPoints::iterator MyPointsIter;
257
258 MyPoints projected_points;
259 MyPoints original_points;
260
261
262 // ------------------------------------------------------------------ main ----
263
264
main(int argc,char ** argv)265 int 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
324 void
open_prog_mesh(const std::string & _filename)325 open_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
455 void
save_vd_prog_mesh(const std::string & _filename)456 save_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
refine(unsigned int _n)591 void 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
coarsen(unsigned int _n)619 void 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
645 void
vdpm_analysis()646 vdpm_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
770 void
get_leaf_node_handles(VHierarchyNodeHandle node_handle,VHierarchyNodeHandleContainer & leaf_nodes)771 get_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
790 void
compute_bounding_box(VHierarchyNodeHandle node_handle,VHierarchyNodeHandleContainer & leaf_nodes)791 compute_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
812 void
compute_cone_of_normals(VHierarchyNodeHandle node_handle,VHierarchyNodeHandleContainer & leaf_nodes)813 compute_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
842 void
compute_screen_space_error(VHierarchyNodeHandle node_handle,VHierarchyNodeHandleContainer & leaf_nodes)843 compute_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
897 void
compute_mue_sigma(VHierarchyNodeHandle node_handle,ResidualContainer & residuals)898 compute_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
955 Vec3f
point2triangle_residual(const Vec3f & p,const Vec3f tri[3],float & s,float & t)956 point2triangle_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 // ============================================================================
1184