1#include <Core/Geometry/CatmullClarkSubdivider.hpp>
2#include <Core/Utils/Index.hpp>
4#include <OpenMesh/Core/Mesh/ArrayKernel.hh>
5#include <OpenMesh/Core/Mesh/AttribKernelT.hh>
6#include <OpenMesh/Core/Mesh/CirculatorsT.hh>
7#include <OpenMesh/Core/Mesh/Handles.hh>
8#include <OpenMesh/Core/Mesh/IteratorsT.hh>
9#include <OpenMesh/Core/Mesh/PolyConnectivity.hh>
10#include <OpenMesh/Core/Mesh/SmartHandles.hh>
11#include <OpenMesh/Core/System/config.h>
12#include <OpenMesh/Tools/Utils/MeshCheckerT.hh>
19bool CatmullClarkSubdivider::prepare( deprecated::TopologicalMesh& mesh ) {
20 mesh.add_property( m_vpPos );
21 mesh.add_property( m_epH );
22 mesh.add_property( m_fpH );
23 mesh.add_property( m_creaseWeights );
24 mesh.createAllPropsOnFaces(
25 m_normalPropF, m_floatPropsF, m_vec2PropsF, m_vec3PropsF, m_vec4PropsF );
26 mesh.add_property( m_hV );
27 for ( uint i = 0; i < mesh.n_halfedges(); ++i ) {
28 auto h = mesh.halfedge_handle( i );
29 if ( !mesh.is_boundary( h ) ) { mesh.property( m_hV, h ) = mesh.to_vertex_handle( h ); }
33 for (
auto e_it = mesh.edges_begin(); e_it != mesh.edges_end(); ++e_it )
34 mesh.property( m_creaseWeights, *e_it ) = 0.0;
39bool CatmullClarkSubdivider::cleanup( deprecated::TopologicalMesh& mesh ) {
40 mesh.remove_property( m_vpPos );
41 mesh.remove_property( m_epH );
42 mesh.remove_property( m_fpH );
43 mesh.remove_property( m_creaseWeights );
44 mesh.clearAllProps( m_normalPropF, m_floatPropsF, m_vec2PropsF, m_vec3PropsF, m_vec4PropsF );
45 mesh.remove_property( m_hV );
49bool CatmullClarkSubdivider::subdivide( deprecated::TopologicalMesh& mesh,
51 const bool update_points ) {
52 m_oldVertexOps.clear();
53 m_newFaceVertexOps.clear();
54 m_newEdgeVertexOps.clear();
55 m_newEdgePropOps.clear();
56 m_newFacePropOps.clear();
57 m_oldVertexOps.resize( n );
58 m_newFaceVertexOps.resize( n );
59 m_newEdgeVertexOps.resize( n );
60 m_newEdgePropOps.resize( n );
61 m_newFacePropOps.resize( n );
63 for (
size_t iter = 0; iter < n; ++iter ) {
65 const size_t NV = mesh.n_vertices();
66 m_newFaceVertexOps[iter].reserve( mesh.n_faces() );
67#pragma omp parallel for
68 for (
int i = 0; i < int( mesh.n_faces() ); ++i ) {
69 const auto& fh = mesh.face_handle( i );
71 deprecated::TopologicalMesh::Point centroid;
72 mesh.calc_face_centroid( fh, centroid );
73 deprecated::TopologicalMesh::VertexHandle vh;
75 { vh = mesh.new_vertex( centroid ); }
76 mesh.property( m_fpH, fh ) = vh;
78 const uint v = mesh.valence( fh );
79 const Scalar inv_v = 1.f / v;
81 auto heh = mesh.halfedge_handle( fh );
82 for ( uint j = 0; j < v; ++j ) {
83 ops[j] = V_OP( inv_v, mesh.to_vertex_handle( heh ) );
84 heh = mesh.next_halfedge_handle( heh );
87 { m_newFaceVertexOps[iter].push_back( V_OPS( vh, ops ) ); }
89 mesh.interpolateAllPropsOnFaces(
90 fh, m_normalPropF, m_floatPropsF, m_vec2PropsF, m_vec3PropsF, m_vec4PropsF );
94 m_newEdgeVertexOps[iter].reserve( mesh.n_edges() );
95#pragma omp parallel for
96 for (
int i = 0; i < int( mesh.n_edges() ); ++i ) {
97 const auto& eh = mesh.edge_handle( i );
98 compute_midpoint( mesh, eh, update_points, iter );
102 if ( update_points ) {
104 m_oldVertexOps[iter].reserve( NV );
105#pragma omp parallel for
106 for (
int i = 0; i < int( NV ); ++i ) {
107 const auto& vh = mesh.vertex_handle( i );
108 update_vertex( mesh, vh, iter );
112#pragma omp parallel for
113 for (
int i = 0; i < int( NV ); ++i ) {
114 const auto& vh = mesh.vertex_handle( i );
115 mesh.set_point( vh, mesh.property( m_vpPos, vh ) );
121 auto e_end = mesh.edges_end();
122 m_newEdgePropOps[iter].reserve( 3 * mesh.n_edges() );
123 for (
auto e_itr = mesh.edges_begin(); e_itr != e_end; ++e_itr ) {
124 split_edge( mesh, *e_itr, iter );
126 m_newEdgePropOps[iter].shrink_to_fit();
130 auto f_end = mesh.faces_end();
131 m_newFacePropOps[iter].reserve( 2 * 4 * mesh.n_faces() );
132 for (
auto f_itr = mesh.faces_begin(); f_itr != f_end; ++f_itr ) {
133 split_face( mesh, *f_itr, iter );
135 m_newFacePropOps[iter].shrink_to_fit();
137 CORE_ASSERT( OpenMesh::Utils::MeshCheckerT<deprecated::TopologicalMesh>( mesh ).check(),
138 "CatmullClarkSubdivision ended with a bad topology." );
144 m_triangulationPropOps.
clear();
145 m_triangulationPropOps.
reserve( 2 * mesh.n_faces() );
147 auto f_end = mesh.faces_end();
148 for (
auto f_itr = mesh.faces_begin(); f_itr != f_end; ++f_itr ) {
150 auto heh1 = mesh.halfedge_handle( *f_itr );
151 auto heh2 = mesh.next_halfedge_handle( heh1 );
152 auto heh3 = mesh.next_halfedge_handle( heh2 );
153 auto heh4 = mesh.next_halfedge_handle( heh3 );
154 auto heh5 = mesh.new_edge( mesh.to_vertex_handle( heh1 ), mesh.to_vertex_handle( heh3 ) );
155 auto heh6 = mesh.opposite_halfedge_handle( heh5 );
158 auto fnew = mesh.new_face();
159 mesh.set_halfedge_handle( fnew, heh2 );
160 mesh.set_face_handle( heh5, *f_itr );
161 mesh.set_face_handle( heh2, fnew );
162 mesh.set_face_handle( heh3, fnew );
163 mesh.set_face_handle( heh6, fnew );
164 mesh.set_next_halfedge_handle( heh1, heh5 );
165 mesh.set_next_halfedge_handle( heh5, heh4 );
166 mesh.set_next_halfedge_handle( heh3, heh6 );
167 mesh.set_next_halfedge_handle( heh6, heh2 );
170 mesh.copyAllProps( heh3, heh5 );
171 mesh.copyAllProps( heh1, heh6 );
172 m_triangulationPropOps.
push_back( { heh5, { { 1, heh3 } } } );
173 m_triangulationPropOps.
push_back( { heh6, { { 1, heh1 } } } );
179void CatmullClarkSubdivider::split_face( deprecated::TopologicalMesh& mesh,
180 const deprecated::TopologicalMesh::FaceHandle& fh,
191 size_t valence = mesh.valence( fh ) / 2;
194 auto vh = mesh.property( m_fpH, fh );
197 auto hend = mesh.halfedge_handle( fh );
198 auto hh = mesh.next_halfedge_handle( hend );
199 auto hold = mesh.new_edge( mesh.to_vertex_handle( hend ), vh );
200 mesh.set_next_halfedge_handle( hend, hold );
201 mesh.set_face_handle( hold, fh );
204 mesh.copyAllPropsFromFace(
205 fh, hold, m_normalPropF, m_floatPropsF, m_vec2PropsF, m_vec3PropsF, m_vec4PropsF );
208 hold = mesh.opposite_halfedge_handle( hold );
211 mesh.copyAllProps( hend, hold );
212 m_newFacePropOps[iter].push_back( { hold, { { 1, hend } } } );
214 const Scalar inv_val = Scalar( 1 ) / valence;
216 p_ops[0] = { inv_val, hh };
218 for (
size_t i = 1; i < valence; i++ ) {
220 auto hnext = mesh.next_halfedge_handle( hh );
221 auto hnew = mesh.new_edge( mesh.to_vertex_handle( hnext ), vh );
224 auto fnew = mesh.new_face();
225 mesh.set_halfedge_handle( fnew, hh );
226 mesh.set_face_handle( hnew, fnew );
227 mesh.set_face_handle( hold, fnew );
228 mesh.set_face_handle( hh, fnew );
229 mesh.set_face_handle( hnext, fnew );
230 mesh.set_next_halfedge_handle( hnew, hold );
231 mesh.set_next_halfedge_handle( hold, hh );
234 mesh.copyAllPropsFromFace(
235 fh, hnew, m_normalPropF, m_floatPropsF, m_vec2PropsF, m_vec3PropsF, m_vec4PropsF );
238 hh = mesh.next_halfedge_handle( hnext );
239 hold = mesh.opposite_halfedge_handle( hnew );
240 mesh.set_next_halfedge_handle( hnext, hnew );
243 mesh.copyAllProps( hnext, hold );
244 m_newFacePropOps[iter].push_back( { hold, { { 1, hnext } } } );
245 p_ops[i] = { inv_val, hh };
249 mesh.set_next_halfedge_handle( hold, hh );
250 mesh.set_next_halfedge_handle( hh, hend );
251 hh = mesh.next_halfedge_handle( hend );
252 mesh.set_next_halfedge_handle( hh, hold );
253 mesh.set_face_handle( hold, fh );
254 mesh.set_halfedge_handle( vh, hold );
257 for (
auto vh_iter = mesh.vih_iter( vh ); vh_iter.is_valid(); ++vh_iter ) {
258 m_newFacePropOps[iter].push_back( { *vh_iter, p_ops } );
262void CatmullClarkSubdivider::split_edge( deprecated::TopologicalMesh& mesh,
263 const deprecated::TopologicalMesh::EdgeHandle& eh,
265 using HeHandle = deprecated::TopologicalMesh::HalfedgeHandle;
266 using VHandle = deprecated::TopologicalMesh::VertexHandle;
269 HeHandle heh = mesh.halfedge_handle( eh, 0 );
270 HeHandle opp_heh = mesh.halfedge_handle( eh, 1 );
272 HeHandle opp_new_heh;
275 VHandle vh1( mesh.to_vertex_handle( heh ) );
278 vh = mesh.property( m_epH, eh );
281 if ( mesh.is_boundary( eh ) ) {
283 for ( t_heh = heh; mesh.next_halfedge_handle( t_heh ) != opp_heh;
284 t_heh = mesh.opposite_halfedge_handle( mesh.next_halfedge_handle( t_heh ) ) )
289 for ( t_heh = mesh.next_halfedge_handle( opp_heh );
290 mesh.next_halfedge_handle( t_heh ) != opp_heh;
291 t_heh = mesh.next_halfedge_handle( t_heh ) )
295 new_heh = mesh.new_edge( vh, vh1 );
296 opp_new_heh = mesh.opposite_halfedge_handle( new_heh );
297 mesh.set_vertex_handle( heh, vh );
299 mesh.set_next_halfedge_handle( t_heh, opp_new_heh );
300 mesh.set_next_halfedge_handle( new_heh, mesh.next_halfedge_handle( heh ) );
301 mesh.set_next_halfedge_handle( heh, new_heh );
302 mesh.set_next_halfedge_handle( opp_new_heh, opp_heh );
304 if ( mesh.face_handle( opp_heh ).is_valid() ) {
305 mesh.set_face_handle( opp_new_heh, mesh.face_handle( opp_heh ) );
306 mesh.set_halfedge_handle( mesh.face_handle( opp_new_heh ), opp_new_heh );
309 mesh.interpolateAllProps( t_heh, opp_heh, opp_new_heh, 0.5 );
310 m_newEdgePropOps[iter].push_back( { opp_new_heh, { { 0.5, t_heh }, { 0.5, opp_heh } } } );
313 if ( mesh.face_handle( heh ).is_valid() ) {
314 mesh.set_face_handle( new_heh, mesh.face_handle( heh ) );
315 mesh.set_halfedge_handle( mesh.face_handle( heh ), heh );
318 mesh.copyAllProps( heh, new_heh );
319 m_newEdgePropOps[iter].push_back( { new_heh, { { 1, heh } } } );
320 HeHandle heh_prev = mesh.prev_halfedge_handle( heh );
321 mesh.interpolateAllProps( heh_prev, new_heh, heh, 0.5 );
322 m_newEdgePropOps[iter].push_back( { heh, { { 0.5, heh_prev }, { 0.5, new_heh } } } );
325 mesh.set_halfedge_handle( vh, new_heh );
326 mesh.set_halfedge_handle( vh1, opp_new_heh );
329 mesh.adjust_outgoing_halfedge( vh );
330 mesh.adjust_outgoing_halfedge( vh1 );
333void CatmullClarkSubdivider::compute_midpoint( deprecated::TopologicalMesh& mesh,
334 const deprecated::TopologicalMesh::EdgeHandle& eh,
335 const bool update_points,
337 deprecated::TopologicalMesh::HalfedgeHandle heh = mesh.halfedge_handle( eh, 0 );
338 deprecated::TopologicalMesh::HalfedgeHandle opp_heh = mesh.halfedge_handle( eh, 1 );
340 deprecated::TopologicalMesh::Point pos = mesh.point( mesh.to_vertex_handle( heh ) );
341 pos += mesh.point( mesh.to_vertex_handle( opp_heh ) );
348 if ( mesh.is_boundary( eh ) || !update_points ) {
352 ops[0] = V_OP( 0.5, mesh.to_vertex_handle( heh ) );
353 ops[1] = V_OP( 0.5, mesh.to_vertex_handle( opp_heh ) );
358 pos += mesh.point( mesh.property( m_fpH, mesh.face_handle( heh ) ) );
359 pos += mesh.point( mesh.property( m_fpH, mesh.face_handle( opp_heh ) ) );
363 ops[0] = V_OP( 0.25, mesh.to_vertex_handle( heh ) );
364 ops[1] = V_OP( 0.25, mesh.to_vertex_handle( opp_heh ) );
365 ops[2] = V_OP( 0.25, mesh.property( m_fpH, mesh.face_handle( heh ) ) );
366 ops[3] = V_OP( 0.25, mesh.property( m_fpH, mesh.face_handle( opp_heh ) ) );
372 auto vh = mesh.new_vertex( pos );
373 mesh.property( m_epH, eh ) = vh;
375 m_newEdgeVertexOps[iter].push_back( V_OPS( vh, ops ) );
379void CatmullClarkSubdivider::update_vertex( deprecated::TopologicalMesh& mesh,
380 const deprecated::TopologicalMesh::VertexHandle& vh,
382 deprecated::TopologicalMesh::Point pos( 0.0, 0.0, 0.0 );
391 if ( mesh.is_boundary( vh ) ) {
392 pos = mesh.point( vh );
394 ops[0] = V_OP( 1.f / 3.f, vh );
396 for (
auto ve_itr = mesh.ve_iter( vh ); ve_itr.is_valid(); ++ve_itr ) {
397 if ( mesh.is_boundary( *ve_itr ) ) {
398 pos += mesh.point( mesh.property( m_epH, *ve_itr ) );
399 ops[i++] = V_OP( 1.f / 3.f, mesh.property( m_epH, *ve_itr ) );
413 const uint valence = mesh.valence( vh );
414 const Scalar inv_v2 = 1.f / ( valence * valence );
415 ops.
resize( valence * 2 + 1 );
418 for (
auto voh_it = mesh.voh_iter( vh ); voh_it.is_valid(); ++voh_it ) {
419 pos += mesh.point( mesh.to_vertex_handle( *voh_it ) );
420 ops[i++] = V_OP( inv_v2, mesh.to_vertex_handle( *voh_it ) );
424 deprecated::TopologicalMesh::Point Q( 0, 0, 0 );
425 for (
auto vf_itr = mesh.vf_iter( vh ); vf_itr.is_valid(); ++vf_itr ) {
426 Q += mesh.point( mesh.property( m_fpH, *vf_itr ) );
427 ops[i++] = V_OP( inv_v2, mesh.property( m_fpH, *vf_itr ) );
432 pos += mesh.point( vh ) * ( valence - 2.0 ) / valence + Q;
433 ops[i] = V_OP( ( valence - 2.0 ) / valence, vh );
437 mesh.property( m_vpPos, vh ) = pos;
440 { m_oldVertexOps[iter].push_back( V_OPS( vh, ops ) ); }
444 const Vector3Array& newCoarseNormals,
445 Vector3Array& newSubdivVertices,
446 Vector3Array& newSubdivNormals,
450 auto hNormalProp = mesh.halfedge_normals_pph();
451#pragma omp parallel for
452 for (
int i = 0; i < int( mesh.n_halfedges() ); ++i ) {
455 auto vh = mesh.property( m_hV, h );
456 if ( vh.idx() != -1 )
458 auto idx = mesh.property( inTriIndexProp, h );
459 mesh.set_point( vh, newCoarseVertices[idx] );
460 mesh.property( hNormalProp, h ) = newCoarseNormals[idx];
464 for (
int i = 0; i < int( m_oldVertexOps.size() ); ++i ) {
466#pragma omp parallel for
467 for (
int j = 0; j < int( m_newFaceVertexOps[i].size() ); ++j ) {
468 Ra::Core::Vector3 pos( 0, 0, 0 );
469 const auto& ops = m_newFaceVertexOps[i][j];
470 for (
const auto& op : ops.second ) {
471 pos += op.first * mesh.point( op.second );
473 mesh.set_point( ops.first, pos );
476#pragma omp parallel for
477 for (
int j = 0; j < int( m_newEdgeVertexOps[i].size() ); ++j ) {
478 Ra::Core::Vector3 pos( 0, 0, 0 );
479 const auto& ops = m_newEdgeVertexOps[i][j];
480 for (
const auto& op : ops.second ) {
481 pos += op.first * mesh.point( op.second );
483 mesh.set_point( ops.first, pos );
487#pragma omp parallel for
488 for (
int j = 0; j < int( m_oldVertexOps[i].size() ); ++j ) {
489 pos[j] = Ra::Core::Vector3( 0, 0, 0 );
490 const auto& ops = m_oldVertexOps[i][j];
491 for (
const auto& op : ops.second ) {
492 pos[j] += op.first * mesh.point( op.second );
496#pragma omp parallel for schedule( static )
497 for (
int j = 0; j < int( m_oldVertexOps[i].size() ); ++j ) {
498 mesh.set_point( m_oldVertexOps[i][j].first, pos[j] );
503 for (
int j = 0; j < int( m_newEdgePropOps[i].size() ); ++j ) {
504 Ra::Core::Vector3 nor( 0, 0, 0 );
505 const auto& ops = m_newEdgePropOps[i][j];
506 for (
const auto& op : ops.second ) {
507 nor += op.first * mesh.property( hNormalProp, op.second );
509 mesh.property( hNormalProp, ops.first ) = nor.normalized();
513#pragma omp parallel for
514 for (
int j = 0; j < int( m_newFacePropOps[i].size() ); ++j ) {
515 Ra::Core::Vector3 nor( 0, 0, 0 );
516 const auto& ops = m_newFacePropOps[i][j];
517 for (
const auto& op : ops.second ) {
518 nor += op.first * mesh.property( hNormalProp, op.second );
520 mesh.property( hNormalProp, ops.first ) = nor.normalized();
525#pragma omp parallel for
526 for (
int j = 0; j < int( m_triangulationPropOps.
size() ); ++j ) {
527 Ra::Core::Vector3 nor( 0, 0, 0 );
528 const auto& ops = m_triangulationPropOps[j];
529 for (
const auto& op : ops.second ) {
530 nor += op.first * mesh.property( hNormalProp, op.second );
532 mesh.property( hNormalProp, ops.first ) = nor.normalized();
536#pragma omp parallel for
537 for (
int i = 0; i < int( mesh.n_halfedges() ); ++i ) {
539 if ( !mesh.is_boundary( h ) ) {
540 auto idx = mesh.property( outTriIndexProp, h );
541 newSubdivVertices[idx] = mesh.point( mesh.to_vertex_handle( h ) );
542 newSubdivNormals[idx] = mesh.property( hNormalProp, h );
void recompute(const Vector3Array &newCoarseVertices, const Vector3Array &newCoarseNormals, Vector3Array &newSubdivVertices, Vector3Array &newSubdivNormals, deprecated::TopologicalMesh &mesh)
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