Developer Documentation
CatmullClarkT_impl.hh
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 // CLASS CatmullClarkT - IMPLEMENTATION
45 //
46 //=============================================================================
47 
48 #define OPENMESH_SUBDIVIDER_UNIFORM_CATMULLCLARK_CC
49 
50 //== INCLUDES =================================================================
51 
52 #include "CatmullClarkT.hh"
53 #include <OpenMesh/Tools/Utils/MeshCheckerT.hh>
54 
55 //== NAMESPACES ===============================================================
56 
57 namespace OpenMesh { // BEGIN_NS_OPENMESH
58 namespace Subdivider { // BEGIN_NS_SUBVIDER
59 namespace Uniform { // BEGIN_NS_UNIFORM
60 
61 //== IMPLEMENTATION ==========================================================
62 
63 template <typename MeshType, typename RealType>
64 bool
66 {
67  _m.add_property( vp_pos_ );
68  _m.add_property( ep_pos_ );
69  _m.add_property( fp_pos_ );
70  _m.add_property( creaseWeights_ );
71 
72  // initialize all weights to 0 (= smooth edge)
73  for( EdgeIter e_it = _m.edges_begin(); e_it != _m.edges_end(); ++e_it)
74  _m.property(creaseWeights_, *e_it ) = 0.0;
75 
76  return true;
77 }
78 
79 //-----------------------------------------------------------------------------
80 
81 template <typename MeshType, typename RealType>
82 bool
84 {
85  _m.remove_property( vp_pos_ );
86  _m.remove_property( ep_pos_ );
87  _m.remove_property( fp_pos_ );
88  _m.remove_property( creaseWeights_ );
89  return true;
90 }
91 
92 //-----------------------------------------------------------------------------
93 
94 template <typename MeshType, typename RealType>
95 bool
96 CatmullClarkT<MeshType,RealType>::subdivide( MeshType& _m , size_t _n , const bool _update_points)
97 {
98  // Do _n subdivisions
99  for ( size_t i = 0; i < _n; ++i)
100  {
101 
102  // Compute face centroid
103  for ( auto fh : _m.faces())
104  {
105  Point centroid;
106  _m.calc_face_centroid( fh, centroid);
107  _m.property( fp_pos_, fh ) = centroid;
108  }
109 
110  // Compute position for new (edge-) vertices and store them in the edge property
111  for ( auto eh : _m.edges())
112  compute_midpoint( _m, eh, _update_points );
113 
114  // position updates activated?
115  if(_update_points)
116  {
117  // compute new positions for old vertices
118  for ( auto vh : _m.vertices())
119  update_vertex( _m, vh );
120 
121  // Commit changes in geometry
122  for ( auto vh : _m.vertices())
123  _m.set_point(vh, _m.property( vp_pos_, vh ) );
124  }
125 
126  // Split each edge at midpoint stored in edge property ep_pos_;
127  // Attention! Creating new edges, hence make sure the loop ends correctly.
128  for ( auto eh : _m.edges())
129  split_edge( _m, eh );
130 
131  // Commit changes in topology and reconsitute consistency
132  // Attention! Creating new faces, hence make sure the loop ends correctly.
133  for ( auto fh : _m.faces())
134  split_face( _m, fh);
135 
136 
137 #if defined(_DEBUG) || defined(DEBUG)
138  // Now we have an consistent mesh!
139  assert( OpenMesh::Utils::MeshCheckerT<MeshType>(_m).check() );
140 #endif
141  }
142 
143  _m.update_normals();
144 
145  return true;
146 }
147 
148 //-----------------------------------------------------------------------------
149 
150 template <typename MeshType, typename RealType>
151 void
152 CatmullClarkT<MeshType,RealType>::split_face( MeshType& _m, const FaceHandle& _fh)
153 {
154  /*
155  Split an n-gon into n quads by connecting
156  each vertex of fh to vh.
157 
158  - _fh will remain valid (it will become one of the quads)
159  - the halfedge handles of the new quads will
160  point to the old halfedges
161  */
162 
163  // Since edges already refined (valence*2)
164  size_t valence = _m.valence(_fh)/2;
165 
166  // new mesh vertex from face centroid
167  VertexHandle vh = _m.add_vertex(_m.property( fp_pos_, _fh ));
168 
169  HalfedgeHandle hend = _m.halfedge_handle(_fh);
170  HalfedgeHandle hh = _m.next_halfedge_handle(hend);
171 
172  HalfedgeHandle hold = _m.new_edge(_m.to_vertex_handle(hend), vh);
173 
174  _m.set_next_halfedge_handle(hend, hold);
175  _m.set_face_handle(hold, _fh);
176 
177  hold = _m.opposite_halfedge_handle(hold);
178 
179  for(size_t i = 1; i < valence; i++)
180  {
181  HalfedgeHandle hnext = _m.next_halfedge_handle(hh);
182 
183  FaceHandle fnew = _m.new_face();
184 
185  _m.set_halfedge_handle(fnew, hh);
186 
187  HalfedgeHandle hnew = _m.new_edge(_m.to_vertex_handle(hnext), vh);
188 
189  _m.set_face_handle(hnew, fnew);
190  _m.set_face_handle(hold, fnew);
191  _m.set_face_handle(hh, fnew);
192  _m.set_face_handle(hnext, fnew);
193 
194  _m.set_next_halfedge_handle(hnew, hold);
195  _m.set_next_halfedge_handle(hold, hh);
196  _m.set_next_halfedge_handle(hh, hnext);
197  hh = _m.next_halfedge_handle(hnext);
198  _m.set_next_halfedge_handle(hnext, hnew);
199 
200  hold = _m.opposite_halfedge_handle(hnew);
201  }
202 
203  _m.set_next_halfedge_handle(hold, hh);
204  _m.set_next_halfedge_handle(hh, hend);
205  hh = _m.next_halfedge_handle(hend);
206  _m.set_next_halfedge_handle(hend, hh);
207  _m.set_next_halfedge_handle(hh, hold);
208 
209  _m.set_face_handle(hold, _fh);
210 
211  _m.set_halfedge_handle(vh, hold);
212 }
213 
214 //-----------------------------------------------------------------------------
215 
216 template <typename MeshType, typename RealType>
217 void
218 CatmullClarkT<MeshType,RealType>::split_edge( MeshType& _m, const EdgeHandle& _eh)
219 {
220  HalfedgeHandle heh = _m.halfedge_handle(_eh, 0);
221  HalfedgeHandle opp_heh = _m.halfedge_handle(_eh, 1);
222 
223  HalfedgeHandle new_heh, opp_new_heh, t_heh;
224  VertexHandle vh;
225  VertexHandle vh1( _m.to_vertex_handle(heh));
226  Point zero(0,0,0);
227 
228  // new vertex
229  vh = _m.new_vertex( zero );
230  _m.set_point( vh, _m.property( ep_pos_, _eh ) );
231 
232  // Re-link mesh entities
233  if (_m.is_boundary(_eh))
234  {
235  for (t_heh = heh;
236  _m.next_halfedge_handle(t_heh) != opp_heh;
237  t_heh = _m.opposite_halfedge_handle(_m.next_halfedge_handle(t_heh)))
238  {}
239  }
240  else
241  {
242  for (t_heh = _m.next_halfedge_handle(opp_heh);
243  _m.next_halfedge_handle(t_heh) != opp_heh;
244  t_heh = _m.next_halfedge_handle(t_heh) )
245  {}
246  }
247 
248  new_heh = _m.new_edge(vh, vh1);
249  opp_new_heh = _m.opposite_halfedge_handle(new_heh);
250  _m.set_vertex_handle( heh, vh );
251 
252  _m.set_next_halfedge_handle(t_heh, opp_new_heh);
253  _m.set_next_halfedge_handle(new_heh, _m.next_halfedge_handle(heh));
254  _m.set_next_halfedge_handle(heh, new_heh);
255  _m.set_next_halfedge_handle(opp_new_heh, opp_heh);
256 
257  if (_m.face_handle(opp_heh).is_valid())
258  {
259  _m.set_face_handle(opp_new_heh, _m.face_handle(opp_heh));
260  _m.set_halfedge_handle(_m.face_handle(opp_new_heh), opp_new_heh);
261  }
262 
263  if( _m.face_handle(heh).is_valid())
264  {
265  _m.set_face_handle( new_heh, _m.face_handle(heh) );
266  _m.set_halfedge_handle( _m.face_handle(heh), heh );
267  }
268 
269  _m.set_halfedge_handle( vh, new_heh);
270  _m.set_halfedge_handle( vh1, opp_new_heh );
271 
272  // Never forget this, when playing with the topology
273  _m.adjust_outgoing_halfedge( vh );
274  _m.adjust_outgoing_halfedge( vh1 );
275 }
276 
277 //-----------------------------------------------------------------------------
278 
279 template <typename MeshType, typename RealType>
280 void
281 CatmullClarkT<MeshType,RealType>::compute_midpoint( MeshType& _m, const EdgeHandle& _eh, const bool _update_points)
282 {
283  HalfedgeHandle heh, opp_heh;
284 
285  heh = _m.halfedge_handle( _eh, 0);
286  opp_heh = _m.halfedge_handle( _eh, 1);
287 
288  Point pos( _m.point( _m.to_vertex_handle( heh)));
289 
290  pos += _m.point( _m.to_vertex_handle( opp_heh));
291 
292  // boundary edge: just average vertex positions
293  // this yields the [1/2 1/2] mask
294  if (_m.is_boundary(_eh) || !_update_points)
295  {
296  pos *= static_cast<RealType>(0.5);
297  }
298 // else if (_m.status(_eh).selected() )
299 // {
300 // pos *= 0.5; // change this
301 // }
302 
303  else // inner edge: add neighbouring Vertices to sum
304  // this yields the [1/16 1/16; 3/8 3/8; 1/16 1/16] mask
305  {
306  pos += _m.property(fp_pos_, _m.face_handle(heh));
307  pos += _m.property(fp_pos_, _m.face_handle(opp_heh));
308  pos *= static_cast<RealType>(0.25);
309  }
310  _m.property( ep_pos_, _eh ) = pos;
311 }
312 
313 //-----------------------------------------------------------------------------
314 
315 template <typename MeshType, typename RealType>
316 void
317 CatmullClarkT<MeshType,RealType>::update_vertex( MeshType& _m, const VertexHandle& _vh)
318 {
319  Point pos(0.0,0.0,0.0);
320 
321  // TODO boundary, Extraordinary Vertex and Creased Surfaces
322  // see "A Factored Approach to Subdivision Surfaces"
323  // http://faculty.cs.tamu.edu/schaefer/research/tutorial.pdf
324  // and http://www.cs.utah.edu/~lacewell/subdeval
325  if ( _m.is_boundary( _vh))
326  {
327  pos = _m.point(_vh);
328  VertexEdgeIter ve_itr;
329  for ( ve_itr = _m.ve_iter( _vh); ve_itr.is_valid(); ++ve_itr)
330  if ( _m.is_boundary( *ve_itr))
331  pos += _m.property( ep_pos_, *ve_itr);
332  pos /= static_cast<typename vector_traits<typename MeshType::Point>::value_type>(3.0);
333  }
334  else // inner vertex
335  {
336  /* For each (non boundary) vertex V, introduce a new vertex whose
337  position is F/n + 2E/n + (n-3)V/n where F is the average of
338  the new face vertices of all faces adjacent to the old vertex
339  V, E is the average of the midpoints of all edges incident
340  on the old vertex V, and n is the number of edges incident on
341  the vertex.
342  */
343 
344  /*
345  Normal Vec;
346  VertexEdgeIter ve_itr;
347  double valence(0.0);
348 
349  // R = Calculate Valence and sum of edge midpoints
350  for ( ve_itr = _m.ve_iter( _vh); ve_itr; ++ve_itr)
351  {
352  valence+=1.0;
353  pos += _m.property(ep_pos_, *ve_itr);
354  }
355  pos /= valence*valence;
356  */
357 
358  RealType valence(0.0);
359  VOHIter voh_it = _m.voh_iter( _vh );
360  for( ; voh_it.is_valid(); ++voh_it )
361  {
362  pos += _m.point( _m.to_vertex_handle( *voh_it ) );
363  valence+=1.0;
364  }
365  pos /= valence*valence;
366 
367  VertexFaceIter vf_itr;
368  Point Q(0, 0, 0);
369 
370  for ( vf_itr = _m.vf_iter( _vh); vf_itr.is_valid(); ++vf_itr) //, neigboring_faces += 1.0 )
371  {
372  Q += _m.property(fp_pos_, *vf_itr);
373  }
374 
375  Q /= valence*valence;//neigboring_faces;
376 
377  pos += _m.point(_vh) * (valence - RealType(2.0) )/valence + Q;
378  // pos = vector_cast<Vec>(_m.point(_vh));
379  }
380 
381  _m.property( vp_pos_, _vh ) = pos;
382 }
383 
384 //-----------------------------------------------------------------------------
385 
386 //=============================================================================
387 } // END_NS_UNIFORM
388 } // END_NS_SUBDIVIDER
389 } // END_NS_OPENMESH
390 //=============================================================================
virtual bool prepare(MeshType &_m) override
Initialize properties and weights.
virtual bool subdivide(MeshType &_m, size_t _n, const bool _update_points=true) override
Execute n subdivision steps.
virtual bool cleanup(MeshType &_m) override
Remove properties and weights.
T::value_type value_type
Type of the scalar value.