Developer Documentation

VoronoiAreaTriMeshT.hpp
1 /*
2 * VoronoiAreaTriMeshTriMeshT.hpp
3 *
4 * Created on: Jul 10, 2014
5 * Author: kremer
6 */
7
8 #ifndef VORONOIAREATRIMESHT_HPP_
9 #define VORONOIAREATRIMESHT_HPP_
10
11 #include <ObjectTypes/TriangleMesh/TriangleMesh.hh>
12
19 template <class MeshT>
20 class VoronoiAreaTriMeshT {
21 public:
22 VoronoiAreaTriMeshT(const MeshT& _mesh, bool _obtuse_handling = true) :
23 mesh_(_mesh),
24 obtuse_handling_(_obtuse_handling) {
25
26 const_cast<MeshT&>(mesh_).add_property(cotan_weights_, "VoronoiAreaTriMeshT::CotanWeights");
27
28 if(obtuse_handling_) {
29
30 const_cast<MeshT&>(mesh_).add_property(obtuse_tag_, "VoronoiAreaTriMeshT::ObtuseTag");
31
32 tag_obtuse_triangles();
33 }
34
35 compute_cotan_weights();
36 }
37
38 virtual ~VoronoiAreaTriMeshT() {
39
40 if(obtuse_handling_) {
41 const_cast<MeshT&>(mesh_).remove_property(obtuse_tag_);
42 }
43
44 const_cast<MeshT&>(mesh_).remove_property(cotan_weights_);
45 }
46
47 const OpenMesh::EPropHandleT<typename MeshT::Scalar>& cotan_weight_prop() const {
48 return cotan_weights_;
49 }
50
51 typename MeshT::Scalar get_voronoi_area(const typename MeshT::VertexHandle& _vh) const {
52
53 typename MeshT::Scalar w = 0.0;
54
55 bool obtuse = false;
56 if(obtuse_handling_) {
57 for(typename MeshT::ConstVertexFaceIter cvf_it = mesh_.cvf_iter(_vh),
58 cvf_end = mesh_.cvf_end(_vh); cvf_it != cvf_end; ++cvf_it) {
59 if(is_obtuse(*cvf_it)) {
60 obtuse = true;
61 break;
62 }
63 }
64 }
65
66 for(typename MeshT::ConstVertexOHalfedgeIter cvoh_it = mesh_.cvoh_iter(_vh),
67 cvoh_end = mesh_.cvoh_end(_vh); cvoh_it != cvoh_end; ++cvoh_it) {
68
69 if(obtuse) {
70
71 // At least one face in one-ring is obtuse
72
73 const typename MeshT::FaceHandle fh = mesh_.face_handle(*cvoh_it);
74
75 // Voronoi inappropriate
76 if(is_obtuse(fh)) {
77
78 w += area(fh) * 0.5;
79
80 } else {
81
82 w += area(fh) * 0.25;
83 }
84
85 } else {
86
87 // Appropriate computation of area
88 w += voronoi_area(*cvoh_it);
89 }
90 }
91
92 return w;
93 }
94
95 private:
96
97 inline typename MeshT::Scalar area(const typename MeshT::FaceHandle& _fh) const {
98
99 typename MeshT::ConstFaceVertexIter cfv_it = mesh_.cfv_iter(_fh);
100
101 const typename MeshT::Point& p0 = mesh_.point(*cfv_it); ++cfv_it;
102 const typename MeshT::Point& p1 = mesh_.point(*cfv_it); ++cfv_it;
103 const typename MeshT::Point& p2 = mesh_.point(*cfv_it);
104
105 return ((p1 - p0) % (p2 - p0)).norm() * 0.5;
106 }
107
108 inline typename MeshT::Scalar voronoi_area(const typename MeshT::HalfedgeHandle& _heh) const {
109
110 if(mesh_.is_boundary(_heh)) return 0.0;
111
112 const typename MeshT::Normal e = mesh_.point(mesh_.to_vertex_handle(_heh)) -
113 mesh_.point(mesh_.from_vertex_handle(_heh));
114
115 return (1.0/8.0) * mesh_.property(cotan_weights_, mesh_.edge_handle(_heh)) * (e | e);
116 }
117
118 void tag_obtuse_triangles() const {
119
120 for(typename MeshT::ConstFaceIter f_it = mesh_.faces_begin(),
121 f_end = mesh_.faces_end(); f_it != f_end; ++f_it) {
122
123 typename MeshT::ConstFaceVertexIter cfv_it = mesh_.cfv_iter(*f_it);
124
125 const typename MeshT::Point& p0 = mesh_.point(*cfv_it); ++cfv_it;
126 const typename MeshT::Point& p1 = mesh_.point(*cfv_it); ++cfv_it;
127 const typename MeshT::Point& p2 = mesh_.point(*cfv_it);
128
129 const_cast<MeshT&>(mesh_).property(obtuse_tag_, *f_it) = (
130 (((p1 - p0)|(p2 - p0)) < 0) ||
131 (((p2 - p1)|(p0 - p1)) < 0) ||
132 (((p0 - p2)|(p1 - p0)) < 0));
133 }
134 }
135
136 // Compute cotan weights
137 void compute_cotan_weights() const {
138
139 for(typename MeshT::ConstEdgeIter e_it = mesh_.edges_begin(),
140 e_end = mesh_.edges_end(); e_it != e_end; ++e_it) {
141
142 const typename MeshT::Scalar w = cotan_weight(mesh_.halfedge_handle(*e_it, 0)) +
143 cotan_weight(mesh_.halfedge_handle(*e_it, 1));
144
145 const_cast<MeshT&>(mesh_).property(cotan_weights_, *e_it) = w;
146 }
147 }
148
149 typename MeshT::Scalar cotan_weight(const typename MeshT::HalfedgeHandle& _he) const {
150
151 /*
152 * Compute Voronoi area for triangle incident to _he
153 */
154 const typename MeshT::Normal e0 = (mesh_.point(mesh_.to_vertex_handle(_he)) -
155 mesh_.point(mesh_.to_vertex_handle(mesh_.next_halfedge_handle(_he)))).normalized();
156
157 const typename MeshT::Normal e1 = (mesh_.point(mesh_.from_vertex_handle(_he)) -
158 mesh_.point(mesh_.to_vertex_handle(mesh_.next_halfedge_handle(_he)))).normalized();
159
160 const typename MeshT::Scalar cos_a = (e0 | e1);
161 const typename MeshT::Scalar sin_a = (e0 % e1).norm();
162
163 return cos_a / sin_a;
164 }
165
166 inline bool is_obtuse(const typename MeshT::FaceHandle& _fh) const {
167 return mesh_.property(obtuse_tag_, _fh);
168 }
169
170 const MeshT& mesh_;
171 const bool obtuse_handling_;
172
173 OpenMesh::FPropHandleT<bool> obtuse_tag_;
174
175 OpenMesh::EPropHandleT<typename MeshT::Scalar> cotan_weights_;
176 };
177
178 #endif /* VORONOIAREATRIMESHT_HPP_ */
OpenMesh::EPropHandleT< typename MeshT::Scalar >
TriangleMesh.hh
VoronoiAreaTriMeshT
Implementation of Voronoi area computation as described in "Discrete Differential-Geometry Operators ...
Definition: VoronoiAreaTriMeshT.hpp:20
OpenMesh::FPropHandleT< bool >