Loading [MathJax]/extensions/TeX/AMSmath.js
Radium Engine  1.6.0
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LoopSubdivider.cpp
1#include <Core/Geometry/LoopSubdivider.hpp>
2#include <Core/Types.hpp>
3#include <Core/Utils/Index.hpp>
4#include <Eigen/Core>
5#include <OpenMesh/Core/Mesh/ArrayKernel.hh>
6#include <OpenMesh/Core/Mesh/AttribKernelT.hh>
7#include <OpenMesh/Core/Mesh/CirculatorsT.hh>
8#include <OpenMesh/Core/Mesh/Handles.hh>
9#include <OpenMesh/Core/Mesh/IteratorsT.hh>
10#include <OpenMesh/Core/Mesh/PolyConnectivity.hh>
11#include <OpenMesh/Core/Mesh/SmartHandles.hh>
12#include <OpenMesh/Tools/Utils/MeshCheckerT.hh>
13#include <memory>
14
15namespace Ra {
16namespace Core {
17namespace Geometry {
18
19bool LoopSubdivider::prepare( deprecated::TopologicalMesh& mesh ) {
20 uint maxValence = 0;
21 for ( auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it ) {
22 if ( mesh.valence( *v_it ) > maxValence ) { maxValence = mesh.valence( *v_it ); }
23 }
24 init_weights( maxValence + 1 );
25 mesh.add_property( m_vpPos );
26 mesh.add_property( m_epPos );
27 mesh.add_property( m_hV );
28 for ( uint i = 0; i < mesh.n_halfedges(); ++i ) {
29 auto h = mesh.halfedge_handle( i );
30 if ( !mesh.is_boundary( h ) ) { mesh.property( m_hV, h ) = mesh.to_vertex_handle( h ); }
31 }
32 return true;
33}
34
35bool LoopSubdivider::cleanup( deprecated::TopologicalMesh& mesh ) {
36 mesh.remove_property( m_vpPos );
37 mesh.remove_property( m_epPos );
38 mesh.remove_property( m_hV );
39 return true;
40}
41
42bool LoopSubdivider::subdivide( deprecated::TopologicalMesh& mesh,
43 size_t n,
44 const bool updatePoints ) {
45 m_oldVertexOps.clear();
46 m_newVertexOps.clear();
47 m_newEdgePropOps.clear();
48 m_newFacePropOps.clear();
49 m_oldVertexOps.resize( n );
50 m_newVertexOps.resize( n );
51 m_newEdgePropOps.resize( n );
52 m_newFacePropOps.resize( n );
53
54 deprecated::TopologicalMesh::FaceIter fit, f_end;
55 deprecated::TopologicalMesh::EdgeIter eit, e_end;
56
57 // Do n subdivisions
58 for ( uint iter = 0; iter < n; ++iter ) {
59 const size_t NV = mesh.n_vertices();
60 if ( updatePoints ) {
61 // compute new positions for old vertices
62 m_oldVertexOps[iter].reserve( NV );
63#pragma omp parallel for
64 for ( int i = 0; i < int( NV ); ++i ) {
65 const auto& vh = mesh.vertex_handle( i );
66 smooth( mesh, vh, iter );
67 }
68 }
69
70 // Compute position for new vertices and store them in the edge property
71 m_newVertexOps[iter].reserve( mesh.n_edges() );
72#pragma omp parallel for
73 for ( int i = 0; i < int( mesh.n_edges() ); ++i ) {
74 const auto& eh = mesh.edge_handle( i );
75 compute_midpoint( mesh, eh, iter );
76 }
77
78 // Split each edge at midpoint and store precomputed positions (stored in
79 // edge property ep_pos_) in the vertex property vp_pos_;
80 // Attention! Creating new edges, hence make sure the loop ends correctly.
81 e_end = mesh.edges_end();
82 m_newEdgePropOps[iter].reserve( 3 * mesh.n_edges() );
83 for ( eit = mesh.edges_begin(); eit != e_end; ++eit ) {
84 split_edge( mesh, *eit, iter );
85 }
86 m_newEdgePropOps[iter].shrink_to_fit();
87
88 // Commit changes in topology and reconsitute consistency
89 // Attention! Creating new faces, hence make sure the loop ends correctly.
90 f_end = mesh.faces_end();
91 m_newFacePropOps[iter].reserve( 6 * mesh.n_faces() );
92 for ( fit = mesh.faces_begin(); fit != f_end; ++fit ) {
93 split_face( mesh, *fit, iter );
94 }
95 m_newFacePropOps[iter].shrink_to_fit();
96
97 if ( updatePoints ) {
98 // Commit changes in geometry
99#pragma omp parallel for
100 for ( int i = 0; i < int( NV ); ++i ) {
101 const auto& vh = mesh.vertex_handle( i );
102 mesh.set_point( vh, mesh.property( m_vpPos, vh ) );
103 }
104 }
105
106 // Now we have an consistent mesh!
107 CORE_ASSERT( OpenMesh::Utils::MeshCheckerT<deprecated::TopologicalMesh>( mesh ).check(),
108 "LoopSubdivision ended with a bad topology." );
109 }
110
111 return true;
112}
113
114void LoopSubdivider::split_face( deprecated::TopologicalMesh& mesh,
115 const deprecated::TopologicalMesh::FaceHandle& fh,
116 size_t iter ) {
117 using HeHandle = deprecated::TopologicalMesh::HalfedgeHandle;
118 // get where to cut
119 HeHandle heh1( mesh.halfedge_handle( fh ) );
120 HeHandle heh2( mesh.next_halfedge_handle( mesh.next_halfedge_handle( heh1 ) ) );
121 HeHandle heh3( mesh.next_halfedge_handle( mesh.next_halfedge_handle( heh2 ) ) );
122
123 // Cutting off every corner of the 6_gon
124 corner_cutting( mesh, heh1, iter );
125 corner_cutting( mesh, heh2, iter );
126 corner_cutting( mesh, heh3, iter );
127}
128
129void LoopSubdivider::corner_cutting( deprecated::TopologicalMesh& mesh,
130 const deprecated::TopologicalMesh::HalfedgeHandle& he,
131 size_t iter ) {
132 using HeHandle = deprecated::TopologicalMesh::HalfedgeHandle;
133 using VHandle = deprecated::TopologicalMesh::VertexHandle;
134 using FHandle = deprecated::TopologicalMesh::FaceHandle;
135
136 // Define Halfedge Handles
137 HeHandle heh1( he );
138 HeHandle heh5( heh1 );
139 HeHandle heh6( mesh.next_halfedge_handle( heh1 ) );
140
141 // Cycle around the polygon to find correct Halfedge
142 for ( ; mesh.next_halfedge_handle( mesh.next_halfedge_handle( heh5 ) ) != heh1;
143 heh5 = mesh.next_halfedge_handle( heh5 ) )
144 ;
145
146 VHandle vh1 = mesh.to_vertex_handle( heh1 );
147 VHandle vh2 = mesh.to_vertex_handle( heh5 );
148
149 HeHandle heh2( mesh.next_halfedge_handle( heh5 ) );
150 HeHandle heh3( mesh.new_edge( vh1, vh2 ) );
151 HeHandle heh4( mesh.opposite_halfedge_handle( heh3 ) );
152
153 /* Intermediate result
154 *
155 * *
156 * 5 /|\
157 * /_ \
158 * vh2> * *
159 * /|\3 |\
160 * /_ \|4 \
161 * *----\*----\*
162 * 1 ^ 6
163 * vh1 (adjust_outgoing halfedge!)
164 */
165
166 // Old and new Face
167 FHandle fh_old( mesh.face_handle( heh6 ) );
168 FHandle fh_new( mesh.new_face() );
169
170 // Re-Set Handles around old Face
171 mesh.set_next_halfedge_handle( heh4, heh6 );
172 mesh.set_next_halfedge_handle( heh5, heh4 );
173
174 mesh.set_face_handle( heh4, fh_old );
175 mesh.set_face_handle( heh5, fh_old );
176 mesh.set_face_handle( heh6, fh_old );
177 mesh.set_halfedge_handle( fh_old, heh4 );
178
179 // Re-Set Handles around new Face
180 mesh.set_next_halfedge_handle( heh1, heh3 );
181 mesh.set_next_halfedge_handle( heh3, heh2 );
182
183 mesh.set_face_handle( heh1, fh_new );
184 mesh.set_face_handle( heh2, fh_new );
185 mesh.set_face_handle( heh3, fh_new );
186
187 mesh.set_halfedge_handle( fh_new, heh1 );
188
189 // deal with custom properties
190 mesh.copyAllProps( heh1, heh4 );
191 mesh.copyAllProps( heh5, heh3 );
192 m_newFacePropOps[iter].push_back( { heh4, { { 1, heh1 } } } );
193 m_newFacePropOps[iter].push_back( { heh3, { { 1, heh5 } } } );
194}
195
196void LoopSubdivider::split_edge( deprecated::TopologicalMesh& mesh,
197 const deprecated::TopologicalMesh::EdgeHandle& eh,
198 size_t iter ) {
199 using HeHandle = deprecated::TopologicalMesh::HalfedgeHandle;
200 using VHandle = deprecated::TopologicalMesh::VertexHandle;
201
202 // prepare data
203 HeHandle heh = mesh.halfedge_handle( eh, 0 );
204 HeHandle opp_heh = mesh.halfedge_handle( eh, 1 );
205 HeHandle new_heh;
206 HeHandle opp_new_heh;
207 HeHandle t_heh;
208 VHandle vh;
209 VHandle vh1( mesh.to_vertex_handle( heh ) );
210
211 // new vertex
212 vh = mesh.property( m_epPos, eh );
213
214 // Re-link mesh entities
215 if ( mesh.is_boundary( eh ) ) {
216 for ( t_heh = heh; mesh.next_halfedge_handle( t_heh ) != opp_heh;
217 t_heh = mesh.opposite_halfedge_handle( mesh.next_halfedge_handle( t_heh ) ) )
218 ;
219 }
220 else {
221 for ( t_heh = mesh.next_halfedge_handle( opp_heh );
222 mesh.next_halfedge_handle( t_heh ) != opp_heh;
223 t_heh = mesh.next_halfedge_handle( t_heh ) )
224 ;
225 }
226
227 new_heh = mesh.new_edge( vh, vh1 );
228 opp_new_heh = mesh.opposite_halfedge_handle( new_heh );
229 mesh.set_vertex_handle( heh, vh );
230
231 mesh.set_next_halfedge_handle( t_heh, opp_new_heh );
232 mesh.set_next_halfedge_handle( new_heh, mesh.next_halfedge_handle( heh ) );
233 mesh.set_next_halfedge_handle( heh, new_heh );
234 mesh.set_next_halfedge_handle( opp_new_heh, opp_heh );
235
236 if ( mesh.face_handle( opp_heh ).is_valid() ) {
237 mesh.set_face_handle( opp_new_heh, mesh.face_handle( opp_heh ) );
238 mesh.set_halfedge_handle( mesh.face_handle( opp_new_heh ), opp_new_heh );
239
240 // deal with custom properties
241 mesh.interpolateAllProps( t_heh, opp_heh, opp_new_heh, 0.5 );
242 m_newEdgePropOps[iter].push_back( { opp_new_heh, { { 0.5, t_heh }, { 0.5, opp_heh } } } );
243 }
244
245 if ( mesh.face_handle( heh ).is_valid() ) {
246 mesh.set_face_handle( new_heh, mesh.face_handle( heh ) );
247 mesh.set_halfedge_handle( mesh.face_handle( heh ), heh );
248
249 // deal with custom properties
250 mesh.copyAllProps( heh, new_heh );
251 m_newEdgePropOps[iter].push_back( { new_heh, { { 1, heh } } } );
252 HeHandle heh_prev = mesh.prev_halfedge_handle( heh );
253 mesh.interpolateAllProps( heh_prev, new_heh, heh, 0.5 );
254 m_newEdgePropOps[iter].push_back( { heh, { { 0.5, heh_prev }, { 0.5, new_heh } } } );
255 }
256
257 mesh.set_halfedge_handle( vh, new_heh );
258 mesh.set_halfedge_handle( vh1, opp_new_heh );
259
260 // Never forget this, when playing with the topology
261 mesh.adjust_outgoing_halfedge( vh );
262 mesh.adjust_outgoing_halfedge( vh1 );
263}
264
265void LoopSubdivider::compute_midpoint( deprecated::TopologicalMesh& mesh,
266 const deprecated::TopologicalMesh::EdgeHandle& eh,
267 size_t iter ) {
268 deprecated::TopologicalMesh::HalfedgeHandle heh = mesh.halfedge_handle( eh, 0 );
269 deprecated::TopologicalMesh::HalfedgeHandle opp_heh = mesh.halfedge_handle( eh, 1 );
270
271 deprecated::TopologicalMesh::Point pos = mesh.point( mesh.to_vertex_handle( heh ) );
272 pos += mesh.point( mesh.to_vertex_handle( opp_heh ) );
273
275
276 // boundary edge: just average vertex positions
277 if ( mesh.is_boundary( eh ) ) {
278 pos *= 0.5;
279 ops.resize( 2 );
280 ops[0] = V_OP( 0.5, mesh.to_vertex_handle( heh ) );
281 ops[1] = V_OP( 0.5, mesh.to_vertex_handle( opp_heh ) );
282 }
283 else // inner edge: add neighbouring Vertices to sum
284 {
285 pos *= 3.0;
286 pos += mesh.point( mesh.to_vertex_handle( mesh.next_halfedge_handle( heh ) ) );
287 pos += mesh.point( mesh.to_vertex_handle( mesh.next_halfedge_handle( opp_heh ) ) );
288 pos *= 1.0 / 8.0;
289 ops.resize( 4 );
290 ops[0] = V_OP( 3.f / 8.f, mesh.to_vertex_handle( heh ) );
291 ops[1] = V_OP( 3.f / 8.f, mesh.to_vertex_handle( opp_heh ) );
292 ops[2] = V_OP( 1.f / 8.f, mesh.to_vertex_handle( mesh.next_halfedge_handle( heh ) ) );
293 ops[3] = V_OP( 1.f / 8.f, mesh.to_vertex_handle( mesh.next_halfedge_handle( opp_heh ) ) );
294 }
295
296#pragma omp critical
297 {
298 auto vh = mesh.add_vertex( pos );
299 mesh.property( m_epPos, eh ) = vh;
300 m_newVertexOps[iter].push_back( V_OPS( vh, ops ) );
301 }
302}
303
304void LoopSubdivider::smooth( deprecated::TopologicalMesh& mesh,
305 const deprecated::TopologicalMesh::VertexHandle& vh,
306 size_t iter ) {
307 using VHandle = deprecated::TopologicalMesh::VertexHandle;
308
309 deprecated::TopologicalMesh::Point pos( 0.0, 0.0, 0.0 );
310
312
313 if ( mesh.is_boundary( vh ) ) // if boundary: Point 1-6-1
314 {
315 deprecated::TopologicalMesh::HalfedgeHandle heh;
316 deprecated::TopologicalMesh::HalfedgeHandle prev_heh;
317 heh = mesh.halfedge_handle( vh );
318
319 assert( mesh.is_boundary( mesh.edge_handle( heh ) ) );
320
321 prev_heh = mesh.prev_halfedge_handle( heh );
322
323 VHandle to_vh = mesh.to_vertex_handle( heh );
324 VHandle from_vh = mesh.from_vertex_handle( prev_heh );
325
326 // ( v_l + 6 v + v_r ) / 8
327 pos = mesh.point( vh );
328 pos *= 6.0;
329 pos += mesh.point( to_vh );
330 pos += mesh.point( from_vh );
331 pos *= 1.0 / 8.0;
332
333 ops.resize( 3 );
334 ops[2] = V_OP( 6.f / 8.f, vh );
335 ops[0] = V_OP( 1.f / 8.f, to_vh );
336 ops[1] = V_OP( 1.f / 8.f, from_vh );
337 }
338 else // inner vertex: (1-a) * p + a/n * Sum q, q in one-ring of p
339 {
340 deprecated::TopologicalMesh::VertexVertexIter vvit;
341 const uint valence = mesh.valence( vh );
342 ops.resize( valence + 1 );
343 int i = 0;
344
345 // Calculate Valence and sum up neighbour points
346 for ( vvit = mesh.cvv_iter( vh ); vvit.is_valid(); ++vvit ) {
347 pos += mesh.point( *vvit );
348 ops[i++] = V_OP( m_weights[valence].second, *vvit );
349 }
350 pos *= m_weights[valence].second; // alpha(n)/n * Sum q, q in one-ring of p
351 pos += m_weights[valence].first * mesh.point( vh ); // + (1-a)*p
352 ops[i] = V_OP( m_weights[valence].first, vh );
353 }
354
355 mesh.property( m_vpPos, vh ) = pos;
356
357#pragma omp critical
358 { m_oldVertexOps[iter].push_back( V_OPS( vh, ops ) ); }
359}
360
361void LoopSubdivider::recompute( const Vector3Array& newCoarseVertices,
362 const Vector3Array& newCoarseNormals,
363 Vector3Array& newSubdivVertices,
364 Vector3Array& newSubdivNormals,
366 // update vertices
367 auto inTriIndexProp = mesh.getInputTriangleMeshIndexPropHandle();
368 auto hNormalProp = mesh.halfedge_normals_pph();
369#pragma omp parallel for
370 for ( int i = 0; i < int( mesh.n_halfedges() ); ++i ) {
371 auto h = mesh.halfedge_handle( i );
372 // set position on coarse mesh vertices
373 auto vh = mesh.property( m_hV, h );
374 if ( vh.idx() != -1 ) // avoid both boundary and non-coarse halfedges
375 {
376 auto idx = mesh.property( inTriIndexProp, h );
377 mesh.set_point( vh, newCoarseVertices[idx] );
378 mesh.property( hNormalProp, h ) = newCoarseNormals[idx];
379 }
380 }
381 // for each subdiv step
382 for ( int i = 0; i < int( m_oldVertexOps.size() ); ++i ) {
383 // first update new vertices
384#pragma omp parallel for schedule( static )
385 for ( int j = 0; j < int( m_newVertexOps[i].size() ); ++j ) {
386 Ra::Core::Vector3 pos( 0, 0, 0 );
387 const auto& ops = m_newVertexOps[i][j];
388 for ( const auto& op : ops.second ) {
389 pos += op.first * mesh.point( op.second );
390 }
391 mesh.set_point( ops.first, pos );
392 }
393 // then compute old vertices
394 std::vector<Ra::Core::Vector3> pos( m_oldVertexOps[i].size() );
395#pragma omp parallel for
396 for ( int j = 0; j < int( m_oldVertexOps[i].size() ); ++j ) {
397 pos[j] = Ra::Core::Vector3( 0, 0, 0 );
398 const auto& ops = m_oldVertexOps[i][j];
399 for ( const auto& op : ops.second ) {
400 pos[j] += op.first * mesh.point( op.second );
401 }
402 }
403 // then commit pos for old vertices
404#pragma omp parallel for schedule( static )
405 for ( int j = 0; j < int( m_oldVertexOps[i].size() ); ++j ) {
406 mesh.set_point( m_oldVertexOps[i][j].first, pos[j] );
407 }
408 // deal with normal on edge centers (other non-static properties can be updated the same
409 // way)
410 // This loop should not be parallelized!
411 for ( int j = 0; j < int( m_newEdgePropOps[i].size() ); ++j ) {
412 Ra::Core::Vector3 nor( 0, 0, 0 );
413 const auto& ops = m_newEdgePropOps[i][j];
414 for ( const auto& op : ops.second ) {
415 nor += op.first * mesh.property( hNormalProp, op.second );
416 }
417 mesh.property( hNormalProp, ops.first ) = nor.normalized();
418 }
419 // deal with normal on faces centers (other non-static properties can be updated the same
420 // way)
421#pragma omp parallel for
422 for ( int j = 0; j < int( m_newFacePropOps[i].size() ); ++j ) {
423 Ra::Core::Vector3 nor( 0, 0, 0 );
424 const auto& ops = m_newFacePropOps[i][j];
425 for ( const auto& op : ops.second ) {
426 nor += op.first * mesh.property( hNormalProp, op.second );
427 }
428 mesh.property( hNormalProp, ops.first ) = nor.normalized();
429 }
430 }
431 // update subdivided TriangleMesh vertices and normals
432 auto outTriIndexProp = mesh.getOutputTriangleMeshIndexPropHandle();
433#pragma omp parallel for
434 for ( int i = 0; i < int( mesh.n_halfedges() ); ++i ) {
435 auto h = mesh.halfedge_handle( i );
436 if ( !mesh.is_boundary( h ) ) {
437 auto idx = mesh.property( outTriIndexProp, h );
438 newSubdivVertices[idx] = mesh.point( mesh.to_vertex_handle( h ) );
439 newSubdivNormals[idx] = mesh.property( hNormalProp, h );
440 }
441 }
442}
443
444} // namespace Geometry
445} // namespace Core
446} // namespace Ra
void recompute(const Vector3Array &newCoarseVertices, const Vector3Array &newCoarseNormals, Vector3Array &newSubdivVertices, Vector3Array &newSubdivNormals, deprecated::TopologicalMesh &mesh)
void init_weights(size_t max_valence)
Pre-compute weights.
const OpenMesh::HPropHandleT< Index > & getInputTriangleMeshIndexPropHandle() const
const OpenMesh::HPropHandleT< Index > & getOutputTriangleMeshIndexPropHandle() const
HalfedgeHandle halfedge_handle(VertexHandle vh, FaceHandle fh) const
@ Geometry
"Geometry" render objects are those loaded using Radium::IO and generated by GeometrySystem
hepler function to manage enum as underlying types in VariableSet
Definition Cage.cpp:4
T resize(T... args)