Developer Documentation
RulesT_impl.hh
Go to the documentation of this file.
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 
48 //=============================================================================
49 //
50 // Rules - IMPLEMENTATION
51 //
52 //=============================================================================
53 
54 
55 #define OPENMESH_SUBDIVIDER_ADAPTIVE_RULEST_CC
56 
57 
58 //== INCLUDES =================================================================
59 
60 #include <OpenMesh/Core/System/config.h>
61 #include <OpenMesh/Core/IO/MeshIO.hh>
62 #include "RulesT.hh"
63 // --------------------
64 #if defined(OM_CC_MIPS)
65 # include <math.h>
66 #else
67 # include <cmath>
68 #endif
69 
70 #if defined(OM_CC_MSVC)
71 # pragma warning(disable:4244)
72 #endif
73 
74 //== NAMESPACE ================================================================
75 
76 namespace OpenMesh { // BEGIN_NS_OPENMESH
77 namespace Subdivider { // BEGIN_NS_DECIMATER
78 namespace Adaptive { // BEGIN_NS_ADAPTIVE
79 
80 
81 //== IMPLEMENTATION ==========================================================
82 
83 #define MOBJ Base::mesh_.data
84 #define FH face_handle
85 #define VH vertex_handle
86 #define EH edge_handle
87 #define HEH halfedge_handle
88 #define NHEH next_halfedge_handle
89 #define PHEH prev_halfedge_handle
90 #define OHEH opposite_halfedge_handle
91 #define TVH to_vertex_handle
92 #define FVH from_vertex_handle
93 
94 // ------------------------------------------------------------------ Tvv3 ----
95 
96 
97 template<class M>
98 void
99 Tvv3<M>::raise(typename M::FaceHandle& _fh, state_t _target_state)
100 {
101  if (MOBJ(_fh).state() < _target_state)
102  {
103  this->update(_fh, _target_state);
104 
105  typename M::VertexVertexIter vv_it;
106  typename M::FaceVertexIter fv_it;
107  typename M::VertexHandle vh;
108  typename M::Point position(0.0, 0.0, 0.0);
109  typename M::Point face_position;
110  const typename M::Point zero_point(0.0, 0.0, 0.0);
111  std::vector<typename M::VertexHandle> vertex_vector;
112 
113 
114  // raise all adjacent vertices to level x-1
115  for (fv_it = Base::mesh_.fv_iter(_fh); fv_it.is_valid(); ++fv_it) {
116 
117  vertex_vector.push_back(*fv_it);
118  }
119 
120  while(!vertex_vector.empty()) {
121 
122  vh = vertex_vector.back();
123  vertex_vector.pop_back();
124 
125  if (_target_state > 1)
126  Base::prev_rule()->raise(vh, _target_state - 1);
127  }
128 
129  face_position = MOBJ(_fh).position(_target_state - 1);
130 
131  typename M::EdgeHandle eh;
132  std::vector<typename M::EdgeHandle> edge_vector;
133 
134  // interior face
135  if (!Base::mesh_.is_boundary(_fh) || MOBJ(_fh).final()) {
136 
137  // insert new vertex
138  vh = Base::mesh_.new_vertex();
139 
140  Base::mesh_.split(_fh, vh);
141 
142  typename M::Scalar valence(0.0);
143 
144  // calculate display position for new vertex
145  for (vv_it = Base::mesh_.vv_iter(vh); vv_it.is_valid(); ++vv_it)
146  {
147  position += Base::mesh_.point(*vv_it);
148  valence += 1.0;
149  }
150 
151  position /= valence;
152 
153  // set attributes for new vertex
154  Base::mesh_.set_point(vh, position);
155  MOBJ(vh).set_position(_target_state, zero_point);
156  MOBJ(vh).set_state(_target_state);
157  MOBJ(vh).set_not_final();
158 
159  typename M::VertexOHalfedgeIter voh_it;
160  // check for edge flipping
161  for (voh_it = Base::mesh_.voh_iter(vh); voh_it.is_valid(); ++voh_it) {
162 
163  if (Base::mesh_.FH(*voh_it).is_valid()) {
164 
165  MOBJ(Base::mesh_.FH(*voh_it)).set_state(_target_state);
166  MOBJ(Base::mesh_.FH(*voh_it)).set_not_final();
167  MOBJ(Base::mesh_.FH(*voh_it)).set_position(_target_state - 1, face_position);
168 
169 
170  for (state_t j = 0; j < _target_state; ++j) {
171  MOBJ(Base::mesh_.FH(*voh_it)).set_position(j, MOBJ(_fh).position(j));
172  }
173 
174  if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it))).is_valid()) {
175 
176  if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it)))).state() == _target_state) {
177 
178  if (Base::mesh_.is_flip_ok(Base::mesh_.EH(Base::mesh_.NHEH(*voh_it)))) {
179 
180  edge_vector.push_back(Base::mesh_.EH(Base::mesh_.NHEH(*voh_it)));
181  }
182  }
183  }
184  }
185  }
186  }
187 
188  // boundary face
189  else {
190 
191  typename M::VertexHandle vh1 = Base::mesh_.new_vertex(),
192  vh2 = Base::mesh_.new_vertex();
193 
194  typename M::HalfedgeHandle hh2 = Base::mesh_.HEH(_fh),
195  hh1, hh3;
196 
197  while (!Base::mesh_.is_boundary(Base::mesh_.OHEH(hh2)))
198  hh2 = Base::mesh_.NHEH(hh2);
199 
200  eh = Base::mesh_.EH(hh2);
201 
202  hh2 = Base::mesh_.NHEH(hh2);
203  hh1 = Base::mesh_.NHEH(hh2);
204 
205  assert(Base::mesh_.is_boundary(eh));
206 
207  Base::mesh_.split(eh, vh1);
208 
209  eh = Base::mesh_.EH(Base::mesh_.PHEH(hh2));
210 
211  assert(Base::mesh_.is_boundary(eh));
212 
213  Base::mesh_.split(eh, vh2);
214 
215  hh3 = Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(hh1)));
216 
217  typename M::VertexHandle vh0(Base::mesh_.TVH(hh1)),
218  vh3(Base::mesh_.FVH(hh2));
219 
220  // set display position and attributes for new vertices
221  Base::mesh_.set_point(vh1, (Base::mesh_.point(vh0) * static_cast<typename M::Point::value_type>(2.0) + Base::mesh_.point(vh3)) / static_cast<typename M::Point::value_type>(3.0) );
222 
223  MOBJ(vh1).set_position(_target_state, zero_point);
224  MOBJ(vh1).set_state(_target_state);
225  MOBJ(vh1).set_not_final();
226 
227  MOBJ(vh0).set_position(_target_state, MOBJ(vh0).position(_target_state - 1) * static_cast<typename M::Point::value_type>(3.0));
228  MOBJ(vh0).set_state(_target_state);
229  MOBJ(vh0).set_not_final();
230 
231  // set display position and attributes for old vertices
232  Base::mesh_.set_point(vh2, (Base::mesh_.point(vh3) * static_cast<typename M::Point::value_type>(2.0) + Base::mesh_.point(vh0)) / static_cast<typename M::Point::value_type>(3.0) );
233  MOBJ(vh2).set_position(_target_state, zero_point);
234  MOBJ(vh2).set_state(_target_state);
235  MOBJ(vh2).set_not_final();
236 
237  MOBJ(vh3).set_position(_target_state, MOBJ(vh3).position(_target_state - 1) * static_cast<typename M::Point::value_type>(3.0) );
238  MOBJ(vh3).set_state(_target_state);
239  MOBJ(vh3).set_not_final();
240 
241  // init 3 faces
242  MOBJ(Base::mesh_.FH(hh1)).set_state(_target_state);
243  MOBJ(Base::mesh_.FH(hh1)).set_not_final();
244  MOBJ(Base::mesh_.FH(hh1)).set_position(_target_state - 1, face_position);
245 
246  MOBJ(Base::mesh_.FH(hh2)).set_state(_target_state);
247  MOBJ(Base::mesh_.FH(hh2)).set_not_final();
248  MOBJ(Base::mesh_.FH(hh2)).set_position(_target_state - 1, face_position);
249 
250  MOBJ(Base::mesh_.FH(hh3)).set_state(_target_state);
251  MOBJ(Base::mesh_.FH(hh3)).set_final();
252  MOBJ(Base::mesh_.FH(hh3)).set_position(_target_state - 1, face_position);
253 
254 
255  for (state_t j = 0; j < _target_state; ++j) {
256  MOBJ(Base::mesh_.FH(hh1)).set_position(j, MOBJ(_fh).position(j));
257  }
258 
259  for (state_t j = 0; j < _target_state; ++j) {
260 
261  MOBJ(Base::mesh_.FH(hh2)).set_position(j, MOBJ(_fh).position(j));
262  }
263 
264  for (state_t j = 0; j < _target_state; ++j) {
265 
266  MOBJ(Base::mesh_.FH(hh3)).set_position(j, MOBJ(_fh).position(j));
267  }
268 
269  // check for edge flipping
270  if (Base::mesh_.FH(Base::mesh_.OHEH(hh1)).is_valid()) {
271 
272  if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(hh1))).state() == _target_state) {
273 
274  if (Base::mesh_.is_flip_ok(Base::mesh_.EH(hh1))) {
275 
276  edge_vector.push_back(Base::mesh_.EH(hh1));
277  }
278  }
279  }
280 
281  if (Base::mesh_.FH(Base::mesh_.OHEH(hh2)).is_valid()) {
282 
283  if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(hh2))).state() == _target_state) {
284 
285  if (Base::mesh_.is_flip_ok(Base::mesh_.EH(hh2))) {
286 
287  edge_vector.push_back(Base::mesh_.EH(hh2));
288  }
289  }
290  }
291  }
292 
293  // flip edges
294  while (!edge_vector.empty()) {
295 
296  eh = edge_vector.back();
297  edge_vector.pop_back();
298 
299  assert(Base::mesh_.is_flip_ok(eh));
300 
301  Base::mesh_.flip(eh);
302 
303  MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 0))).set_final();
304  MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 1))).set_final();
305 
306  MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 0))).set_state(_target_state);
307  MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 1))).set_state(_target_state);
308 
309  MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 0))).set_position(_target_state, face_position);
310  MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 1))).set_position(_target_state, face_position);
311  }
312  }
313 }
314 
315 
316 template<class M>
317 void Tvv3<M>::raise(typename M::VertexHandle& _vh, state_t _target_state) {
318 
319  if (MOBJ(_vh).state() < _target_state) {
320 
321  this->update(_vh, _target_state);
322 
323  // multiply old position by 3
324  MOBJ(_vh).set_position(_target_state, MOBJ(_vh).position(_target_state - 1) * static_cast<typename M::Point::value_type>(3.0) );
325 
326  MOBJ(_vh).inc_state();
327 
328  assert(MOBJ(_vh).state() == _target_state);
329  }
330 }
331 
332 
333 // ------------------------------------------------------------------ Tvv4 ----
334 
335 
336 template<class M>
337 void
338 Tvv4<M>::raise(typename M::FaceHandle& _fh, state_t _target_state)
339 {
340 
341  if (MOBJ(_fh).state() < _target_state) {
342 
343  this->update(_fh, _target_state);
344 
345  typename M::FaceVertexIter fv_it;
346  typename M::VertexHandle temp_vh;
347  typename M::Point face_position;
348  const typename M::Point zero_point(0.0, 0.0, 0.0);
349  std::vector<typename M::VertexHandle> vertex_vector;
350  std::vector<typename M::HalfedgeHandle> halfedge_vector;
351 
352  // raise all adjacent vertices to level x-1
353  for (fv_it = Base::mesh_.fv_iter(_fh); fv_it.is_valid(); ++fv_it) {
354 
355  vertex_vector.push_back(*fv_it);
356  }
357 
358  while(!vertex_vector.empty()) {
359 
360  temp_vh = vertex_vector.back();
361  vertex_vector.pop_back();
362 
363  if (_target_state > 1) {
364  Base::prev_rule()->raise(temp_vh, _target_state - 1);
365  }
366  }
367 
368  face_position = MOBJ(_fh).position(_target_state - 1);
369 
370  typename M::HalfedgeHandle hh[3];
371  typename M::VertexHandle vh[3];
372  typename M::VertexHandle new_vh[3];
373  typename M::FaceHandle fh[4];
374  typename M::EdgeHandle eh;
375  typename M::HalfedgeHandle temp_hh;
376 
377  // normal (final) face
378  if (MOBJ(_fh).final()) {
379 
380  // define three halfedge handles around the face
381  hh[0] = Base::mesh_.HEH(_fh);
382  hh[1] = Base::mesh_.NHEH(hh[0]);
383  hh[2] = Base::mesh_.NHEH(hh[1]);
384 
385  assert(hh[0] == Base::mesh_.NHEH(hh[2]));
386 
387  vh[0] = Base::mesh_.TVH(hh[0]);
388  vh[1] = Base::mesh_.TVH(hh[1]);
389  vh[2] = Base::mesh_.TVH(hh[2]);
390 
391  new_vh[0] = Base::mesh_.add_vertex(zero_point);
392  new_vh[1] = Base::mesh_.add_vertex(zero_point);
393  new_vh[2] = Base::mesh_.add_vertex(zero_point);
394 
395  // split three edges
396  split_edge(hh[0], new_vh[0], _target_state);
397  eh = Base::mesh_.EH(Base::mesh_.PHEH(hh[2]));
398  split_edge(hh[1], new_vh[1], _target_state);
399  split_edge(hh[2], new_vh[2], _target_state);
400 
401  assert(Base::mesh_.FVH(hh[2]) == vh[1]);
402  assert(Base::mesh_.FVH(hh[1]) == vh[0]);
403  assert(Base::mesh_.FVH(hh[0]) == vh[2]);
404 
405  if (Base::mesh_.FH(Base::mesh_.OHEH(hh[0])).is_valid())
406  {
407  temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[0])));
408  if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
409  halfedge_vector.push_back(temp_hh);
410  }
411 
412  if (Base::mesh_.FH(Base::mesh_.OHEH(hh[1])).is_valid()) {
413 
414  temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[1])));
415  if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
416  halfedge_vector.push_back(temp_hh);
417  }
418 
419  if (Base::mesh_.FH(Base::mesh_.OHEH(hh[2])).is_valid()) {
420 
421  temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[2])));
422  if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
423  halfedge_vector.push_back(temp_hh);
424  }
425  }
426 
427  // splitted face, check for type
428  else {
429 
430  // define red halfedge handle
431  typename M::HalfedgeHandle red_hh(MOBJ(_fh).red_halfedge());
432 
433  if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh))).is_valid()
434  && Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)))).is_valid()
435  && MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)))).red_halfedge() == red_hh
436  && MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh))))).red_halfedge() == red_hh)
437  {
438 
439  // three times divided face
440  vh[0] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)))));
441  vh[1] = Base::mesh_.TVH(red_hh);
442  vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh))));
443 
444  new_vh[0] = Base::mesh_.FVH(red_hh);
445  new_vh[1] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)));
446  new_vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(red_hh));
447 
448  hh[0] = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)));
449  hh[1] = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh))));
450  hh[2] = Base::mesh_.NHEH(red_hh);
451 
452  eh = Base::mesh_.EH(red_hh);
453  }
454 
455  else
456  {
457 
458  if ((Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh))).is_valid() &&
459  MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)))).red_halfedge()
460  == red_hh )
461  || (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)))).is_valid()
462  && MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh))))).red_halfedge() == red_hh))
463  {
464 
465  // double divided face
466  if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)))).red_halfedge() == red_hh)
467  {
468  // first case
469  vh[0] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)));
470  vh[1] = Base::mesh_.TVH(red_hh);
471  vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh))));
472 
473  new_vh[0] = Base::mesh_.FVH(red_hh);
474  new_vh[1] = Base::mesh_.add_vertex(zero_point);
475  new_vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(red_hh));
476 
477  hh[0] = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)));
478  hh[1] = Base::mesh_.PHEH(Base::mesh_.OHEH(red_hh));
479  hh[2] = Base::mesh_.NHEH(red_hh);
480 
481  // split one edge
482  eh = Base::mesh_.EH(red_hh);
483 
484  split_edge(hh[1], new_vh[1], _target_state);
485 
486  assert(Base::mesh_.FVH(hh[2]) == vh[1]);
487  assert(Base::mesh_.FVH(hh[1]) == vh[0]);
488  assert(Base::mesh_.FVH(hh[0]) == vh[2]);
489 
490  if (Base::mesh_.FH(Base::mesh_.OHEH(hh[1])).is_valid())
491  {
492  temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[1])));
493  if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
494  halfedge_vector.push_back(temp_hh);
495  }
496  }
497  else
498  {
499 
500  // second case
501  vh[0] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)))));
502  vh[1] = Base::mesh_.TVH(red_hh);
503  vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(red_hh));
504 
505  new_vh[0] = Base::mesh_.FVH(red_hh);
506  new_vh[1] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)));
507  new_vh[2] = Base::mesh_.add_vertex(zero_point);
508 
509  hh[0] = Base::mesh_.PHEH(red_hh);
510  hh[1] = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh))));
511  hh[2] = Base::mesh_.NHEH(red_hh);
512 
513  // split one edge
514  eh = Base::mesh_.EH(red_hh);
515 
516  split_edge(hh[2], new_vh[2], _target_state);
517 
518  assert(Base::mesh_.FVH(hh[2]) == vh[1]);
519  assert(Base::mesh_.FVH(hh[1]) == vh[0]);
520  assert(Base::mesh_.FVH(hh[0]) == vh[2]);
521 
522  if (Base::mesh_.FH(Base::mesh_.OHEH(hh[2])).is_valid()) {
523 
524  temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[2])));
525  if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
526  halfedge_vector.push_back(temp_hh);
527  }
528  }
529  }
530 
531  else {
532 
533  // one time divided face
534  vh[0] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)));
535  vh[1] = Base::mesh_.TVH(red_hh);
536  vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(red_hh));
537 
538  new_vh[0] = Base::mesh_.FVH(red_hh);
539  new_vh[1] = Base::mesh_.add_vertex(zero_point);
540  new_vh[2] = Base::mesh_.add_vertex(zero_point);
541 
542  hh[0] = Base::mesh_.PHEH(red_hh);
543  hh[1] = Base::mesh_.PHEH(Base::mesh_.OHEH(red_hh));
544  hh[2] = Base::mesh_.NHEH(red_hh);
545 
546  // split two edges
547  eh = Base::mesh_.EH(red_hh);
548 
549  split_edge(hh[1], new_vh[1], _target_state);
550  split_edge(hh[2], new_vh[2], _target_state);
551 
552  assert(Base::mesh_.FVH(hh[2]) == vh[1]);
553  assert(Base::mesh_.FVH(hh[1]) == vh[0]);
554  assert(Base::mesh_.FVH(hh[0]) == vh[2]);
555 
556  if (Base::mesh_.FH(Base::mesh_.OHEH(hh[1])).is_valid()) {
557 
558  temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[1])));
559  if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
560  halfedge_vector.push_back(temp_hh);
561  }
562 
563  if (Base::mesh_.FH(Base::mesh_.OHEH(hh[2])).is_valid()) {
564 
565  temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[2])));
566  if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
567  halfedge_vector.push_back(temp_hh);
568  }
569  }
570  }
571  }
572 
573  // continue here for all cases
574 
575  // flip edge
576  if (Base::mesh_.is_flip_ok(eh)) {
577 
578  Base::mesh_.flip(eh);
579  }
580 
581  // search new faces
582  fh[0] = Base::mesh_.FH(hh[0]);
583  fh[1] = Base::mesh_.FH(hh[1]);
584  fh[2] = Base::mesh_.FH(hh[2]);
585  fh[3] = Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(hh[0])));
586 
587  assert(_fh == fh[0] || _fh == fh[1] || _fh == fh[2] || _fh == fh[3]);
588 
589  // init new faces
590  for (int i = 0; i <= 3; ++i) {
591 
592  MOBJ(fh[i]).set_state(_target_state);
593  MOBJ(fh[i]).set_final();
594  MOBJ(fh[i]).set_position(_target_state, face_position);
595  MOBJ(fh[i]).set_red_halfedge(Base::mesh_.InvalidHalfedgeHandle);
596  }
597 
598  // init new vertices and edges
599  for (int i = 0; i <= 2; ++i) {
600 
601  MOBJ(new_vh[i]).set_position(_target_state, zero_point);
602  MOBJ(new_vh[i]).set_state(_target_state);
603  MOBJ(new_vh[i]).set_not_final();
604 
605  Base::mesh_.set_point(new_vh[i], (Base::mesh_.point(vh[i]) + Base::mesh_.point(vh[(i + 2) % 3])) * 0.5);
606 
607  MOBJ(Base::mesh_.EH(hh[i])).set_state(_target_state);
608  MOBJ(Base::mesh_.EH(hh[i])).set_position(_target_state, zero_point);
609  MOBJ(Base::mesh_.EH(hh[i])).set_final();
610 
611  MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh[i]))).set_state(_target_state);
612  MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh[i]))).set_position(_target_state, zero_point);
613  MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh[i]))).set_final();
614 
615  MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh[i]))).set_state(_target_state);
616  MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh[i]))).set_position(_target_state, zero_point);
617  MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh[i]))).set_final();
618  }
619 
620  // check, if opposite triangle needs splitting
621  while (!halfedge_vector.empty()) {
622 
623  temp_hh = halfedge_vector.back();
624  halfedge_vector.pop_back();
625 
626  check_edge(temp_hh, _target_state);
627  }
628 
629  assert(MOBJ(fh[0]).state() == _target_state);
630  assert(MOBJ(fh[1]).state() == _target_state);
631  assert(MOBJ(fh[2]).state() == _target_state);
632  assert(MOBJ(fh[3]).state() == _target_state);
633  }
634 }
635 
636 
637 template<class M>
638 void
639 Tvv4<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
640 {
641 
642  if (MOBJ(_vh).state() < _target_state)
643  {
644 
645  this->update(_vh, _target_state);
646 
647  // multiply old position by 4
648  MOBJ(_vh).set_position(_target_state, MOBJ(_vh).position(_target_state - 1) * static_cast<typename M::Point::value_type>(4.0));
649 
650  MOBJ(_vh).inc_state();
651  }
652 }
653 
654 
655 template<class M>
656 void
657 Tvv4<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state)
658 {
659  if (MOBJ(_eh).state() < _target_state)
660  {
661  this->update(_eh, _target_state);
662 
663  typename M::FaceHandle fh(Base::mesh_.FH(Base::mesh_.HEH(_eh, 0)));
664 
665  if (!fh.is_valid())
666  fh=Base::mesh_.FH(Base::mesh_.HEH(_eh, 1));
667 
668  raise(fh, _target_state);
669 
670  assert(MOBJ(_eh).state() == _target_state);
671  }
672 }
673 
674 #ifndef DOXY_IGNORE_THIS
675 
676 template<class M>
677 void
678 Tvv4<M>::split_edge(typename M::HalfedgeHandle &_hh,
679  typename M::VertexHandle &_vh,
680  state_t _target_state)
681  {
682  typename M::HalfedgeHandle temp_hh;
683 
684  if (Base::mesh_.FH(Base::mesh_.OHEH(_hh)).is_valid())
685  {
686  if (!MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).final())
687  {
688  if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).red_halfedge().is_valid())
689  {
690  temp_hh = MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).red_halfedge();
691  }
692  else
693  {
694  // two cases for divided, but not visited face
695  if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh))))).state()
696  == MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).state())
697  {
698  temp_hh = Base::mesh_.PHEH(Base::mesh_.OHEH(_hh));
699  }
700 
701  else if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(_hh))))).state()
702  == MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).state())
703  {
704  temp_hh = Base::mesh_.NHEH(Base::mesh_.OHEH(_hh));
705  }
706  }
707  }
708  else
709  temp_hh = Base::mesh_.InvalidHalfedgeHandle;
710  }
711  else
712  temp_hh = Base::mesh_.InvalidHalfedgeHandle;
713 
714  // split edge
715  Base::mesh_.split(Base::mesh_.EH(_hh), _vh);
716 
717  if (Base::mesh_.FVH(_hh) == _vh)
718  {
719  MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(_hh))))).set_state(MOBJ(Base::mesh_.EH(_hh)).state());
720  _hh = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(_hh)));
721  }
722 
723  if (Base::mesh_.FH(Base::mesh_.OHEH(_hh)).is_valid()) {
724 
725  MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh)))).set_not_final();
726  MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).set_state(_target_state-1);
727  MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh)))))).set_state(_target_state-1);
728 
729  MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).set_not_final();
730  MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh)))))).set_not_final();
731 
732  MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh)))).set_state(_target_state);
733 
734  if (temp_hh.is_valid()) {
735 
736  MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).set_red_halfedge(temp_hh);
737  MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh)))))).set_red_halfedge(temp_hh);
738  }
739  else {
740 
741  typename M::FaceHandle
742  fh1(Base::mesh_.FH(Base::mesh_.OHEH(_hh))),
743  fh2(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh))))));
744 
745  MOBJ(fh1).set_red_halfedge(Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh))));
746  MOBJ(fh2).set_red_halfedge(Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh))));
747 
748  const typename M::Point zero_point(0.0, 0.0, 0.0);
749 
750  MOBJ(fh1).set_position(_target_state - 1, zero_point);
751  MOBJ(fh2).set_position(_target_state - 1, zero_point);
752  }
753  }
754 
755  // init edges
756  MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh))))).set_state(_target_state - 1);
757  MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh))))).set_final();
758 
759  MOBJ(Base::mesh_.EH(_hh)).set_state(_target_state - 1);
760  MOBJ(Base::mesh_.EH(_hh)).set_final();
761 }
762 
763 
764 template<class M>
765 void Tvv4<M>::check_edge(const typename M::HalfedgeHandle& _hh,
766  state_t _target_state)
767 {
768  typename M::FaceHandle fh1(Base::mesh_.FH(_hh)),
769  fh2(Base::mesh_.FH(Base::mesh_.OHEH(_hh)));
770 
771  assert(fh1.is_valid());
772  assert(fh2.is_valid());
773 
774  typename M::HalfedgeHandle red_hh(MOBJ(fh1).red_halfedge());
775 
776  if (!MOBJ(fh1).final()) {
777 
778  assert (MOBJ(fh1).final() == MOBJ(fh2).final());
779  assert (!MOBJ(fh1).final());
780  assert (MOBJ(fh1).red_halfedge() == MOBJ(fh2).red_halfedge());
781 
782  const typename M::Point zero_point(0.0, 0.0, 0.0);
783 
784  MOBJ(fh1).set_position(_target_state - 1, zero_point);
785  MOBJ(fh2).set_position(_target_state - 1, zero_point);
786 
787  assert(red_hh.is_valid());
788 
789  if (!red_hh.is_valid()) {
790 
791  MOBJ(fh1).set_state(_target_state - 1);
792  MOBJ(fh2).set_state(_target_state - 1);
793 
794  MOBJ(fh1).set_red_halfedge(_hh);
795  MOBJ(fh2).set_red_halfedge(_hh);
796 
797  MOBJ(Base::mesh_.EH(_hh)).set_not_final();
798  MOBJ(Base::mesh_.EH(_hh)).set_state(_target_state - 1);
799  }
800 
801  else {
802 
803  MOBJ(Base::mesh_.EH(_hh)).set_not_final();
804  MOBJ(Base::mesh_.EH(_hh)).set_state(_target_state - 1);
805 
806  raise(fh1, _target_state);
807 
808  assert(MOBJ(fh1).state() == _target_state);
809  }
810  }
811 }
812 
813 
814 // -------------------------------------------------------------------- VF ----
815 
816 
817 template<class M>
818 void VF<M>::raise(typename M::FaceHandle& _fh, state_t _target_state)
819 {
820  if (MOBJ(_fh).state() < _target_state) {
821 
822  this->update(_fh, _target_state);
823 
824  // raise all neighbor vertices to level x-1
825  typename M::FaceVertexIter fv_it;
826  typename M::VertexHandle vh;
827  std::vector<typename M::VertexHandle> vertex_vector;
828 
829  if (_target_state > 1) {
830 
831  for (fv_it = Base::mesh_.fv_iter(_fh); fv_it.is_valid(); ++fv_it) {
832 
833  vertex_vector.push_back(*fv_it);
834  }
835 
836  while (!vertex_vector.empty()) {
837 
838  vh = vertex_vector.back();
839  vertex_vector.pop_back();
840 
841  Base::prev_rule()->raise(vh, _target_state - 1);
842  }
843  }
844 
845  // calculate new position
846  typename M::Point position(0.0, 0.0, 0.0);
847  typename M::Scalar valence(0.0);
848 
849  for (fv_it = Base::mesh_.fv_iter(_fh); fv_it.is_valid(); ++fv_it) {
850 
851  valence += 1.0;
852  position += Base::mesh_.data(*fv_it).position(_target_state - 1);
853  }
854 
855  position /= valence;
856 
857  // boundary rule
858  if (Base::number() == Base::subdiv_rule()->number() + 1 &&
859  Base::mesh_.is_boundary(_fh) &&
860  !MOBJ(_fh).final())
861  position *= static_cast<typename M::Scalar>(0.5);
862 
863  MOBJ(_fh).set_position(_target_state, position);
864  MOBJ(_fh).inc_state();
865 
866  assert(_target_state == MOBJ(_fh).state());
867  }
868 }
869 
870 
871 // -------------------------------------------------------------------- FF ----
872 
873 
874 template<class M>
875 void FF<M>::raise(typename M::FaceHandle& _fh, state_t _target_state) {
876 
877  if (MOBJ(_fh).state() < _target_state) {
878 
879  this->update(_fh, _target_state);
880 
881  // raise all neighbor faces to level x-1
882  typename M::FaceFaceIter ff_it;
883  typename M::FaceHandle fh;
884  std::vector<typename M::FaceHandle> face_vector;
885 
886  if (_target_state > 1) {
887 
888  for (ff_it = Base::mesh_.ff_iter(_fh); ff_it.is_valid(); ++ff_it) {
889 
890  face_vector.push_back(*ff_it);
891  }
892 
893  while (!face_vector.empty()) {
894 
895  fh = face_vector.back();
896  face_vector.pop_back();
897 
898  Base::prev_rule()->raise(fh, _target_state - 1);
899  }
900 
901  for (ff_it = Base::mesh_.ff_iter(_fh); ff_it.is_valid(); ++ff_it) {
902 
903  face_vector.push_back(*ff_it);
904  }
905 
906  while (!face_vector.empty()) {
907 
908  fh = face_vector.back();
909  face_vector.pop_back();
910 
911  while (MOBJ(fh).state() < _target_state - 1)
912  Base::prev_rule()->raise(fh, _target_state - 1);
913  }
914  }
915 
916  // calculate new position
917  typename M::Point position(0.0, 0.0, 0.0);
918  typename M::Scalar valence(0.0);
919 
920  for (ff_it = Base::mesh_.ff_iter(_fh); ff_it.is_valid(); ++ff_it) {
921 
922  valence += 1.0;
923 
924  position += Base::mesh_.data(*ff_it).position(_target_state - 1);
925  }
926 
927  position /= valence;
928 
929  MOBJ(_fh).set_position(_target_state, position);
930  MOBJ(_fh).inc_state();
931  }
932 }
933 
934 
935 // ------------------------------------------------------------------- FFc ----
936 
937 
938 template<class M>
939 void FFc<M>::raise(typename M::FaceHandle& _fh, state_t _target_state)
940 {
941  if (MOBJ(_fh).state() < _target_state) {
942 
943  this->update(_fh, _target_state);
944 
945  // raise all neighbor faces to level x-1
946  typename M::FaceFaceIter ff_it(Base::mesh_.ff_iter(_fh));
947  typename M::FaceHandle fh;
948  std::vector<typename M::FaceHandle> face_vector;
949 
950  if (_target_state > 1)
951  {
952  for (; ff_it.is_valid(); ++ff_it)
953  face_vector.push_back(*ff_it);
954 
955  while (!face_vector.empty())
956  {
957  fh = face_vector.back();
958  face_vector.pop_back();
959  Base::prev_rule()->raise(fh, _target_state - 1);
960  }
961 
962  for (ff_it = Base::mesh_.ff_iter(_fh); ff_it.is_valid(); ++ff_it)
963  face_vector.push_back(*ff_it);
964 
965  while (!face_vector.empty()) {
966 
967  fh = face_vector.back();
968  face_vector.pop_back();
969 
970  while (MOBJ(fh).state() < _target_state - 1)
971  Base::prev_rule()->raise(fh, _target_state - 1);
972  }
973  }
974 
975  // calculate new position
976  typename M::Point position(0.0, 0.0, 0.0);
977  typename M::Scalar valence(0.0);
978 
979  for (ff_it = Base::mesh_.ff_iter(_fh); ff_it.is_valid(); ++ff_it)
980  {
981  valence += 1.0;
982  position += Base::mesh_.data(*ff_it).position(_target_state - 1);
983  }
984 
985  position /= valence;
986 
987  // choose coefficient c
988  typename M::Scalar c = Base::coeff();
989 
990  position *= (static_cast<typename M::Scalar>(1.0) - c);
991  position += MOBJ(_fh).position(_target_state - 1) * c;
992 
993  MOBJ(_fh).set_position(_target_state, position);
994  MOBJ(_fh).inc_state();
995  }
996 }
997 
998 
999 // -------------------------------------------------------------------- FV ----
1000 
1001 
1002 template<class M>
1003 void FV<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
1004 {
1005 
1006  if (MOBJ(_vh).state() < _target_state) {
1007 
1008  this->update(_vh, _target_state);
1009 
1010  // raise all neighbor vertices to level x-1
1011  typename M::VertexFaceIter vf_it(Base::mesh_.vf_iter(_vh));
1012  typename M::FaceHandle fh;
1013  std::vector<typename M::FaceHandle> face_vector;
1014 
1015  if (_target_state > 1) {
1016 
1017  for (; vf_it.is_valid(); ++vf_it) {
1018 
1019  face_vector.push_back(*vf_it);
1020  }
1021 
1022  while (!face_vector.empty()) {
1023 
1024  fh = face_vector.back();
1025  face_vector.pop_back();
1026 
1027  Base::prev_rule()->raise(fh, _target_state - 1);
1028  }
1029 
1030  for (vf_it = Base::mesh_.vf_iter(_vh); vf_it.is_valid(); ++vf_it) {
1031 
1032  face_vector.push_back(*vf_it);
1033  }
1034 
1035  while (!face_vector.empty()) {
1036 
1037  fh = face_vector.back();
1038  face_vector.pop_back();
1039 
1040  while (MOBJ(fh).state() < _target_state - 1)
1041  Base::prev_rule()->raise(fh, _target_state - 1);
1042  }
1043  }
1044 
1045  // calculate new position
1046  typename M::Point position(0.0, 0.0, 0.0);
1047  typename M::Scalar valence(0.0);
1048 
1049  for (vf_it = Base::mesh_.vf_iter(_vh); vf_it.is_valid(); ++vf_it) {
1050 
1051  valence += 1.0;
1052  position += Base::mesh_.data(*vf_it).position(_target_state - 1);
1053  }
1054 
1055  position /= valence;
1056 
1057  MOBJ(_vh).set_position(_target_state, position);
1058  MOBJ(_vh).inc_state();
1059 
1060  if (Base::number() == Base::n_rules() - 1) {
1061 
1062  Base::mesh_.set_point(_vh, position);
1063  MOBJ(_vh).set_final();
1064  }
1065  }
1066 }
1067 
1068 
1069 // ------------------------------------------------------------------- FVc ----
1070 
1071 
1072 template<class M>
1073 void FVc<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
1074 {
1075  if (MOBJ(_vh).state() < _target_state) {
1076 
1077  this->update(_vh, _target_state);
1078 
1079  typename M::VertexOHalfedgeIter voh_it;
1080  typename M::FaceHandle fh;
1081  std::vector<typename M::FaceHandle> face_vector;
1082  int valence(0);
1083 
1084  face_vector.clear();
1085 
1086  // raise all neighbour faces to level x-1
1087  if (_target_state > 1) {
1088 
1089  for (voh_it = Base::mesh_.voh_iter(_vh); voh_it.is_valid(); ++voh_it) {
1090 
1091  if (Base::mesh_.FH(*voh_it).is_valid()) {
1092 
1093  face_vector.push_back(Base::mesh_.FH(*voh_it));
1094 
1095  if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it))).is_valid()) {
1096 
1097  face_vector.push_back(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it))));
1098  }
1099  }
1100  }
1101 
1102  while (!face_vector.empty()) {
1103 
1104  fh = face_vector.back();
1105  face_vector.pop_back();
1106 
1107  Base::prev_rule()->raise(fh, _target_state - 1);
1108  }
1109 
1110  for (voh_it = Base::mesh_.voh_iter(_vh); voh_it.is_valid(); ++voh_it) {
1111 
1112  if (Base::mesh_.FH(*voh_it).is_valid()) {
1113 
1114  face_vector.push_back(Base::mesh_.FH(*voh_it));
1115 
1116  if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it))).is_valid()) {
1117 
1118  face_vector.push_back(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it))));
1119  }
1120  }
1121  }
1122 
1123  while (!face_vector.empty()) {
1124 
1125  fh = face_vector.back();
1126  face_vector.pop_back();
1127 
1128  while (MOBJ(fh).state() < _target_state - 1)
1129  Base::prev_rule()->raise(fh, _target_state - 1);
1130  }
1131 
1132  for (voh_it = Base::mesh_.voh_iter(_vh); voh_it.is_valid(); ++voh_it) {
1133 
1134  if (Base::mesh_.FH(*voh_it).is_valid()) {
1135 
1136  face_vector.push_back(Base::mesh_.FH(*voh_it));
1137 
1138  if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it))).is_valid()) {
1139 
1140  face_vector.push_back(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it))));
1141  }
1142  }
1143  }
1144 
1145  while (!face_vector.empty()) {
1146 
1147  fh = face_vector.back();
1148  face_vector.pop_back();
1149 
1150  while (MOBJ(fh).state() < _target_state - 1)
1151  Base::prev_rule()->raise(fh, _target_state - 1);
1152  }
1153  }
1154 
1155  // calculate new position
1156  typename M::Point position(0.0, 0.0, 0.0);
1157  typename M::Scalar c;
1158 #if 0
1159  const typename M::Scalar _2pi(2.0*M_PI);
1160  const typename M::Scalar _2over3(2.0/3.0);
1161 
1162  for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it)
1163  {
1164  ++valence;
1165  }
1166 
1167  // choose coefficient c
1168  c = _2over3 * ( cos( _2pi / valence) + 1.0);
1169 #else
1170  valence = Base::mesh_.valence(_vh);
1171  c = typename M::Scalar(coeff(valence));
1172 #endif
1173 
1174 
1175  for (voh_it = Base::mesh_.voh_iter(_vh); voh_it.is_valid(); ++voh_it) {
1176 
1177  fh = Base::mesh_.FH(*voh_it);
1178  if (fh.is_valid())
1179  Base::prev_rule()->raise(fh, _target_state - 1);
1180 
1181  fh = Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it)));
1182  if (fh.is_valid())
1183  Base::prev_rule()->raise(fh, _target_state - 1);
1184 
1185  if (Base::mesh_.FH(*voh_it).is_valid()) {
1186 
1187  if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it))).is_valid()) {
1188 
1189  position += MOBJ(Base::mesh_.FH(*voh_it)).position(_target_state - 1) * c;
1190 
1191  position += MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(*voh_it)))).position(_target_state - 1) * ( typename M::Scalar(1.0) - c);
1192  }
1193  else {
1194 
1195  position += MOBJ(Base::mesh_.FH(*voh_it)).position(_target_state - 1);
1196  }
1197  }
1198 
1199  else {
1200 
1201  --valence;
1202  }
1203  }
1204 
1205  position /= typename M::Scalar(valence);
1206 
1207  MOBJ(_vh).set_position(_target_state, position);
1208  MOBJ(_vh).inc_state();
1209 
1210  assert(MOBJ(_vh).state() == _target_state);
1211 
1212  // check if last rule
1213  if (Base::number() == Base::n_rules() - 1) {
1214 
1215  Base::mesh_.set_point(_vh, position);
1216  MOBJ(_vh).set_final();
1217  }
1218  }
1219 }
1220 
1221 template<class M>
1222 std::vector<double> FVc<M>::coeffs_;
1223 
1224 template <class M>
1225 void FVc<M>::init_coeffs(size_t _max_valence)
1226 {
1227  if ( coeffs_.size() == _max_valence+1 )
1228  return;
1229 
1230  if ( coeffs_.size() < _max_valence+1 )
1231  {
1232  const double _2pi(2.0*M_PI);
1233  const double _2over3(2.0/3.0);
1234 
1235  if (coeffs_.empty())
1236  coeffs_.push_back(0.0); // dummy for valence 0
1237 
1238  for(size_t v=coeffs_.size(); v <= _max_valence; ++v)
1239  coeffs_.push_back(_2over3 * ( cos( _2pi / v) + 1.0));
1240  }
1241 }
1242 
1243 
1244 // -------------------------------------------------------------------- VV ----
1245 
1246 
1247 template<class M>
1248 void VV<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
1249 {
1250  if (MOBJ(_vh).state() < _target_state)
1251  {
1252  this->update(_vh, _target_state);
1253 
1254  // raise all neighbor vertices to level x-1
1255  typename M::VertexVertexIter vv_it(Base::mesh_.vv_iter(_vh));
1256  typename M::VertexHandle vh;
1257  std::vector<typename M::VertexHandle> vertex_vector;
1258 
1259  if (_target_state > 1) {
1260 
1261  for (; vv_it.is_valid(); ++vv_it) {
1262 
1263  vertex_vector.push_back(*vv_it);
1264  }
1265 
1266  while (!vertex_vector.empty()) {
1267 
1268  vh = vertex_vector.back();
1269  vertex_vector.pop_back();
1270 
1271  Base::prev_rule()->raise(vh, _target_state - 1);
1272  }
1273 
1274  for (; vv_it.is_valid(); ++vv_it) {
1275 
1276  vertex_vector.push_back(*vv_it);
1277  }
1278 
1279  while (!vertex_vector.empty()) {
1280 
1281  vh = vertex_vector.back();
1282  vertex_vector.pop_back();
1283 
1284  Base::prev_rule()->raise(vh, _target_state - 1);
1285  }
1286  }
1287 
1288  // calculate new position
1289  typename M::Point position(0.0, 0.0, 0.0);
1290  typename M::Scalar valence(0.0);
1291 
1292  for (vv_it = Base::mesh_.vv_iter(_vh); vv_it.is_valid(); ++vv_it) {
1293 
1294  valence += 1.0;
1295  position += Base::mesh_.data(*vv_it).position(_target_state - 1);
1296  }
1297 
1298  position /= valence;
1299 
1300  MOBJ(_vh).set_position(_target_state, position);
1301  MOBJ(_vh).inc_state();
1302 
1303  // check if last rule
1304  if (Base::number() == Base::n_rules() - 1) {
1305 
1306  Base::mesh_.set_point(_vh, position);
1307  MOBJ(_vh).set_final();
1308  }
1309  }
1310 }
1311 
1312 
1313 // ------------------------------------------------------------------- VVc ----
1314 
1315 
1316 template<class M>
1317 void VVc<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
1318 {
1319  if (MOBJ(_vh).state() < _target_state) {
1320 
1321  this->update(_vh, _target_state);
1322 
1323  // raise all neighbor vertices to level x-1
1324  typename M::VertexVertexIter vv_it(Base::mesh_.vv_iter(_vh));
1325  typename M::VertexHandle vh;
1326  std::vector<typename M::VertexHandle> vertex_vector;
1327 
1328  if (_target_state > 1) {
1329 
1330  for (; vv_it.is_valid(); ++vv_it) {
1331 
1332  vertex_vector.push_back(*vv_it);
1333  }
1334 
1335  while (!vertex_vector.empty()) {
1336 
1337  vh = vertex_vector.back();
1338  vertex_vector.pop_back();
1339 
1340  Base::prev_rule()->raise(vh, _target_state - 1);
1341  }
1342 
1343  for (; vv_it.is_valid(); ++vv_it) {
1344 
1345  vertex_vector.push_back(*vv_it);
1346  }
1347 
1348  while (!vertex_vector.empty()) {
1349 
1350  vh = vertex_vector.back();
1351  vertex_vector.pop_back();
1352 
1353  Base::prev_rule()->raise(vh, _target_state - 1);
1354  }
1355  }
1356 
1357  // calculate new position
1358  typename M::Point position(0.0, 0.0, 0.0);
1359  typename M::Scalar valence(0.0);
1360  typename M::Scalar c;
1361 
1362  for (vv_it = Base::mesh_.vv_iter(_vh); vv_it.is_valid(); ++vv_it)
1363  {
1364  valence += 1.0;
1365  position += Base::mesh_.data(*vv_it).position(_target_state - 1);
1366  }
1367 
1368  position /= valence;
1369 
1370  // choose coefficient c
1371  c = Base::coeff();
1372 
1373  position *= (static_cast<typename M::Scalar>(1.0) - c);
1374  position += MOBJ(_vh).position(_target_state - 1) * c;
1375 
1376  MOBJ(_vh).set_position(_target_state, position);
1377  MOBJ(_vh).inc_state();
1378 
1379  if (Base::number() == Base::n_rules() - 1)
1380  {
1381  Base::mesh_.set_point(_vh, position);
1382  MOBJ(_vh).set_final();
1383  }
1384  }
1385 }
1386 
1387 
1388 // -------------------------------------------------------------------- VE ----
1389 
1390 
1391 template<class M>
1392 void VE<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state)
1393 {
1394  if (MOBJ(_eh).state() < _target_state) {
1395 
1396  this->update(_eh, _target_state);
1397 
1398  // raise all neighbour vertices to level x-1
1399  typename M::VertexHandle vh;
1400  typename M::HalfedgeHandle hh1(Base::mesh_.HEH(_eh, 0)),
1401  hh2(Base::mesh_.HEH(_eh, 1));
1402 
1403  if (_target_state > 1) {
1404 
1405  vh = Base::mesh_.TVH(hh1);
1406 
1407  Base::prev_rule()->raise(vh, _target_state - 1);
1408 
1409  vh = Base::mesh_.TVH(hh2);
1410 
1411  Base::prev_rule()->raise(vh, _target_state - 1);
1412  }
1413 
1414  // calculate new position
1415  typename M::Point position(0.0, 0.0, 0.0);
1416  const typename M::Scalar valence(2.0);
1417 
1418  position += MOBJ(Base::mesh_.TVH(hh1)).position(_target_state - 1);
1419  position += MOBJ(Base::mesh_.TVH(hh2)).position(_target_state - 1);
1420 
1421  position /= valence;
1422 
1423  MOBJ(_eh).set_position(_target_state, position);
1424  MOBJ(_eh).inc_state();
1425  }
1426 }
1427 
1428 
1429 // ------------------------------------------------------------------- VdE ----
1430 
1431 
1432 template<class M>
1433 void VdE<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state)
1434 {
1435  if (MOBJ(_eh).state() < _target_state)
1436  {
1437  this->update(_eh, _target_state);
1438 
1439  // raise all neighbor vertices to level x-1
1440  typename M::VertexHandle vh;
1441  typename M::HalfedgeHandle hh1(Base::mesh_.HEH(_eh, 0)),
1442  hh2(Base::mesh_.HEH(_eh, 1));
1443  typename M::FaceHandle fh1, fh2;
1444 
1445  if (_target_state > 1) {
1446 
1447  fh1 = Base::mesh_.FH(hh1);
1448  fh2 = Base::mesh_.FH(hh2);
1449 
1450  if (fh1.is_valid()) {
1451 
1452  Base::prev_rule()->raise(fh1, _target_state - 1);
1453 
1454  vh = Base::mesh_.TVH(Base::mesh_.NHEH(hh1));
1455  Base::prev_rule()->raise(vh, _target_state - 1);
1456  }
1457 
1458  if (fh2.is_valid()) {
1459 
1460  Base::prev_rule()->raise(fh2, _target_state - 1);
1461 
1462  vh = Base::mesh_.TVH(Base::mesh_.NHEH(hh2));
1463  Base::prev_rule()->raise(vh, _target_state - 1);
1464  }
1465 
1466  vh = Base::mesh_.TVH(hh1);
1467  Base::prev_rule()->raise(vh, _target_state - 1);
1468 
1469  vh = Base::mesh_.TVH(hh2);
1470  Base::prev_rule()->raise(vh, _target_state - 1);
1471  }
1472 
1473  // calculate new position
1474  typename M::Point position(0.0, 0.0, 0.0);
1475  typename M::Scalar valence(2.0);
1476 
1477  position += MOBJ(Base::mesh_.TVH(hh1)).position(_target_state - 1);
1478  position += MOBJ(Base::mesh_.TVH(hh2)).position(_target_state - 1);
1479 
1480  if (fh1.is_valid()) {
1481 
1482  position += MOBJ(Base::mesh_.TVH(Base::mesh_.NHEH(hh1))).position(_target_state - 1);
1483  valence += 1.0;
1484  }
1485 
1486  if (fh2.is_valid()) {
1487 
1488  position += MOBJ(Base::mesh_.TVH(Base::mesh_.NHEH(hh2))).position(_target_state - 1);
1489  valence += 1.0;
1490  }
1491 
1492  if (Base::number() == Base::subdiv_rule()->Base::number() + 1)
1493  valence = 4.0;
1494 
1495  position /= valence;
1496 
1497  MOBJ(_eh).set_position(_target_state, position);
1498  MOBJ(_eh).inc_state();
1499  }
1500 }
1501 
1502 
1503 // ------------------------------------------------------------------ VdEc ----
1504 
1505 
1506 template<class M>
1507 void
1508 VdEc<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state)
1509 {
1510  if (MOBJ(_eh).state() < _target_state)
1511  {
1512  this->update(_eh, _target_state);
1513 
1514  // raise all neighbor vertices to level x-1
1515  typename M::VertexHandle vh;
1516  typename M::HalfedgeHandle hh1(Base::mesh_.HEH(_eh, 0)),
1517  hh2(Base::mesh_.HEH(_eh, 1));
1518  std::vector<typename M::VertexHandle> vertex_vector;
1519  typename M::FaceHandle fh1, fh2;
1520 
1521  if (_target_state > 1) {
1522 
1523  fh1 = Base::mesh_.FH(Base::mesh_.HEH(_eh, 0));
1524  fh2 = Base::mesh_.FH(Base::mesh_.HEH(_eh, 1));
1525 
1526  Base::prev_rule()->raise(fh1, _target_state - 1);
1527  Base::prev_rule()->raise(fh2, _target_state - 1);
1528 
1529  vertex_vector.push_back(Base::mesh_.TVH(hh1));
1530  vertex_vector.push_back(Base::mesh_.TVH(hh2));
1531 
1532  vertex_vector.push_back(Base::mesh_.TVH(Base::mesh_.NHEH(hh1)));
1533  vertex_vector.push_back(Base::mesh_.TVH(Base::mesh_.NHEH(hh2)));
1534 
1535  while (!vertex_vector.empty()) {
1536 
1537  vh = vertex_vector.back();
1538  vertex_vector.pop_back();
1539 
1540  Base::prev_rule()->raise(vh, _target_state - 1);
1541  }
1542 
1543  vertex_vector.push_back(Base::mesh_.TVH(hh1));
1544  vertex_vector.push_back(Base::mesh_.TVH(hh2));
1545 
1546  vertex_vector.push_back(Base::mesh_.TVH(Base::mesh_.NHEH(hh1)));
1547  vertex_vector.push_back(Base::mesh_.TVH(Base::mesh_.NHEH(hh2)));
1548 
1549  while (!vertex_vector.empty()) {
1550 
1551  vh = vertex_vector.back();
1552  vertex_vector.pop_back();
1553 
1554  Base::prev_rule()->raise(vh, _target_state - 1);
1555  }
1556  }
1557 
1558  // calculate new position
1559  typename M::Point position(0.0, 0.0, 0.0);
1560  const typename M::Scalar valence(4.0);
1561  typename M::Scalar c;
1562 
1563  // choose coefficient c
1564  c = Base::coeff();
1565 
1566  position += MOBJ(Base::mesh_.TVH(hh1)).position(_target_state - 1) * c;
1567  position += MOBJ(Base::mesh_.TVH(hh2)).position(_target_state - 1) * c;
1568  position += MOBJ(Base::mesh_.TVH(Base::mesh_.NHEH(hh1))).position(_target_state - 1) * (0.5 - c);
1569  position += MOBJ(Base::mesh_.TVH(Base::mesh_.NHEH(hh2))).position(_target_state - 1) * (0.5 - c);
1570 
1571  position /= valence;
1572 
1573  MOBJ(_eh).set_position(_target_state, position);
1574  MOBJ(_eh).inc_state();
1575  }
1576 }
1577 
1578 
1579 // -------------------------------------------------------------------- EV ----
1580 
1581 
1582 template<class M>
1583 void EV<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
1584 {
1585  if (MOBJ(_vh).state() < _target_state) {
1586 
1587  this->update(_vh, _target_state);
1588 
1589  // raise all neighbor vertices to level x-1
1590  typename M::VertexEdgeIter ve_it(Base::mesh_.ve_iter(_vh));
1591  typename M::EdgeHandle eh;
1592  std::vector<typename M::EdgeHandle> edge_vector;
1593 
1594  if (_target_state > 1) {
1595 
1596  for (; ve_it.is_valid(); ++ve_it) {
1597 
1598  edge_vector.push_back(*ve_it);
1599  }
1600 
1601  while (!edge_vector.empty()) {
1602 
1603  eh = edge_vector.back();
1604  edge_vector.pop_back();
1605 
1606  Base::prev_rule()->raise(eh, _target_state - 1);
1607  }
1608 
1609  for (ve_it = Base::mesh_.ve_iter(_vh); ve_it.is_valid(); ++ve_it) {
1610 
1611  edge_vector.push_back(*ve_it);
1612  }
1613 
1614  while (!edge_vector.empty()) {
1615 
1616  eh = edge_vector.back();
1617  edge_vector.pop_back();
1618 
1619  while (MOBJ(eh).state() < _target_state - 1)
1620  Base::prev_rule()->raise(eh, _target_state - 1);
1621  }
1622  }
1623 
1624  // calculate new position
1625  typename M::Point position(0.0, 0.0, 0.0);
1626  typename M::Scalar valence(0.0);
1627 
1628  for (ve_it = Base::mesh_.ve_iter(_vh); ve_it.is_valid(); ++ve_it) {
1629 
1630  if (Base::mesh_.data(*ve_it).final()) {
1631 
1632  valence += 1.0;
1633 
1634  position += Base::mesh_.data(*ve_it).position(_target_state - 1);
1635  }
1636  }
1637 
1638  position /= valence;
1639 
1640  MOBJ(_vh).set_position(_target_state, position);
1641  MOBJ(_vh).inc_state();
1642 
1643  // check if last rule
1644  if (Base::number() == Base::n_rules() - 1) {
1645 
1646  Base::mesh_.set_point(_vh, position);
1647  MOBJ(_vh).set_final();
1648  }
1649  }
1650 }
1651 
1652 
1653 // ------------------------------------------------------------------- EVc ----
1654 
1655 template<class M>
1656 std::vector<double> EVc<M>::coeffs_;
1657 
1658 template<class M>
1659 void EVc<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
1660 {
1661  if (MOBJ(_vh).state() < _target_state)
1662  {
1663  this->update(_vh, _target_state);
1664 
1665  // raise all neighbour vertices to level x-1
1666  typename M::VertexOHalfedgeIter voh_it;
1667  typename M::EdgeHandle eh;
1668  typename M::FaceHandle fh;
1669  std::vector<typename M::EdgeHandle> edge_vector;
1670  std::vector<typename M::FaceHandle> face_vector;
1671 
1672  if (_target_state > 1) {
1673 
1674  for (voh_it = Base::mesh_.voh_iter(_vh); voh_it.is_valid(); ++voh_it) {
1675 
1676  face_vector.push_back(Base::mesh_.FH(*voh_it));
1677  }
1678 
1679  while (!face_vector.empty()) {
1680 
1681  fh = face_vector.back();
1682  face_vector.pop_back();
1683 
1684  if (fh.is_valid())
1685  Base::prev_rule()->raise(fh, _target_state - 1);
1686  }
1687 
1688  for (voh_it = Base::mesh_.voh_iter(_vh); voh_it.is_valid(); ++voh_it) {
1689 
1690  edge_vector.push_back(Base::mesh_.EH(*voh_it));
1691 
1692  edge_vector.push_back(Base::mesh_.EH(Base::mesh_.NHEH(*voh_it)));
1693  }
1694 
1695  while (!edge_vector.empty()) {
1696 
1697  eh = edge_vector.back();
1698  edge_vector.pop_back();
1699 
1700  while (MOBJ(eh).state() < _target_state - 1)
1701  Base::prev_rule()->raise(eh, _target_state - 1);
1702  }
1703  }
1704 
1705 
1706  // calculate new position
1707  typename M::Point position(0.0, 0.0, 0.0);
1708  typename M::Scalar c;
1709  typename M::Point zero_point(0.0, 0.0, 0.0);
1710  size_t valence(0);
1711 
1712  valence = Base::mesh_.valence(_vh);
1713  c = static_cast<typename M::Scalar>(coeff( valence ));
1714 
1715  for (voh_it = Base::mesh_.voh_iter(_vh); voh_it.is_valid(); ++voh_it)
1716  {
1717  if (MOBJ(Base::mesh_.EH(*voh_it)).final())
1718  {
1719  position += MOBJ(Base::mesh_.EH(*voh_it)).position(_target_state-1)*c;
1720 
1721  if ( Base::mesh_.FH(*voh_it).is_valid() &&
1722  MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(*voh_it))).final() &&
1723  MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(*voh_it))).position(_target_state - 1) != zero_point)
1724  {
1725  position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(*voh_it))).position(_target_state-1) * (1.0-c);
1726  }
1727  else {
1728  position += MOBJ(Base::mesh_.EH(*voh_it)).position(_target_state - 1) * (1.0 - c);
1729  }
1730  }
1731  else {
1732  --valence;
1733  }
1734  }
1735 
1736  position /= valence;
1737 
1738  MOBJ(_vh).set_position(_target_state, position);
1739  MOBJ(_vh).inc_state();
1740 
1741  // check if last rule
1742  if (Base::number() == Base::n_rules() - 1) {
1743 
1744  Base::mesh_.set_point(_vh, position);
1745  MOBJ(_vh).set_final();
1746  }
1747  }
1748 }
1749 
1750 
1751 template <class M>
1752 void
1753 EVc<M>::init_coeffs(size_t _max_valence)
1754 {
1755  if ( coeffs_.size() == _max_valence+1 ) // equal? do nothing
1756  return;
1757 
1758  if (coeffs_.size() < _max_valence+1) // less than? add additional valences
1759  {
1760  const double _2pi = 2.0*M_PI;
1761 
1762  if (coeffs_.empty())
1763  coeffs_.push_back(0.0); // dummy for invalid valences 0,1,2
1764 
1765  for(size_t v=coeffs_.size(); v <= _max_valence; ++v)
1766  {
1767  // ( 3/2 + cos ( 2 PI / valence ) )� / 2 - 1
1768  double c = 1.5 + cos( _2pi / v );
1769  c = c * c * 0.5 - 1.0;
1770  coeffs_.push_back(c);
1771  }
1772  }
1773 }
1774 
1775 
1776 // -------------------------------------------------------------------- EF ----
1777 
1778 template<class M>
1779 void
1780 EF<M>::raise(typename M::FaceHandle& _fh, state_t _target_state) {
1781 
1782  if (MOBJ(_fh).state() < _target_state) {
1783 
1784  this->update(_fh, _target_state);
1785 
1786  // raise all neighbour edges to level x-1
1787  typename M::FaceEdgeIter fe_it(Base::mesh_.fe_iter(_fh));
1788  typename M::EdgeHandle eh;
1789  std::vector<typename M::EdgeHandle> edge_vector;
1790 
1791  if (_target_state > 1) {
1792 
1793  for (; fe_it.is_valid(); ++fe_it) {
1794 
1795  edge_vector.push_back(*fe_it);
1796  }
1797 
1798  while (!edge_vector.empty()) {
1799 
1800  eh = edge_vector.back();
1801  edge_vector.pop_back();
1802 
1803  Base::prev_rule()->raise(eh, _target_state - 1);
1804  }
1805 
1806  for (fe_it = Base::mesh_.fe_iter(_fh); fe_it.is_valid(); ++fe_it) {
1807 
1808  edge_vector.push_back(*fe_it);
1809  }
1810 
1811  while (!edge_vector.empty()) {
1812 
1813  eh = edge_vector.back();
1814  edge_vector.pop_back();
1815 
1816  while (MOBJ(eh).state() < _target_state - 1)
1817  Base::prev_rule()->raise(eh, _target_state - 1);
1818  }
1819  }
1820 
1821  // calculate new position
1822  typename M::Point position(0.0, 0.0, 0.0);
1823  typename M::Scalar valence(0.0);
1824 
1825  for (fe_it = Base::mesh_.fe_iter(_fh); fe_it.is_valid(); ++fe_it) {
1826 
1827  if (Base::mesh_.data(*fe_it).final()) {
1828 
1829  valence += 1.0;
1830 
1831  position += Base::mesh_.data(*fe_it).position(_target_state - 1);
1832  }
1833  }
1834 
1835  assert (valence == 3.0);
1836 
1837  position /= valence;
1838 
1839  MOBJ(_fh).set_position(_target_state, position);
1840  MOBJ(_fh).inc_state();
1841  }
1842 }
1843 
1844 
1845 // -------------------------------------------------------------------- FE ----
1846 
1847 
1848 template<class M>
1849 void
1850 FE<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) {
1851 
1852  if (MOBJ(_eh).state() < _target_state) {
1853 
1854  this->update(_eh, _target_state);
1855 
1856  // raise all neighbor faces to level x-1
1857  typename M::FaceHandle fh;
1858 
1859  if (_target_state > 1) {
1860 
1861  fh = Base::mesh_.FH(Base::mesh_.HEH(_eh, 0));
1862  Base::prev_rule()->raise(fh, _target_state - 1);
1863 
1864  fh = Base::mesh_.FH(Base::mesh_.HEH(_eh, 1));
1865  Base::prev_rule()->raise(fh, _target_state - 1);
1866 
1867  fh = Base::mesh_.FH(Base::mesh_.HEH(_eh, 0));
1868  Base::prev_rule()->raise(fh, _target_state - 1);
1869 
1870  fh = Base::mesh_.FH(Base::mesh_.HEH(_eh, 1));
1871  Base::prev_rule()->raise(fh, _target_state - 1);
1872  }
1873 
1874  // calculate new position
1875  typename M::Point position(0.0, 0.0, 0.0);
1876  const typename M::Scalar valence(2.0);
1877 
1878  position += MOBJ(Base::mesh_.FH(Base::mesh_.HEH(_eh, 0))).position(_target_state - 1);
1879 
1880  position += MOBJ(Base::mesh_.FH(Base::mesh_.HEH(_eh, 1))).position(_target_state - 1);
1881 
1882  position /= valence;
1883 
1884  MOBJ(_eh).set_position(_target_state, position);
1885  MOBJ(_eh).inc_state();
1886  }
1887 }
1888 
1889 
1890 // ------------------------------------------------------------------- EdE ----
1891 
1892 
1893 template<class M>
1894 void
1895 EdE<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) {
1896 
1897  if (MOBJ(_eh).state() < _target_state) {
1898 
1899  this->update(_eh, _target_state);
1900 
1901  // raise all neighbor faces and edges to level x-1
1902  typename M::HalfedgeHandle hh1, hh2;
1903  typename M::FaceHandle fh;
1904  typename M::EdgeHandle eh;
1905 
1906  hh1 = Base::mesh_.HEH(_eh, 0);
1907  hh2 = Base::mesh_.HEH(_eh, 1);
1908 
1909  if (_target_state > 1) {
1910 
1911  fh = Base::mesh_.FH(hh1);
1912  Base::prev_rule()->raise(fh, _target_state - 1);
1913 
1914  fh = Base::mesh_.FH(hh2);
1915  Base::prev_rule()->raise(fh, _target_state - 1);
1916 
1917  eh = Base::mesh_.EH(Base::mesh_.NHEH(hh1));
1918  Base::prev_rule()->raise(eh, _target_state - 1);
1919 
1920  eh = Base::mesh_.EH(Base::mesh_.PHEH(hh1));
1921  Base::prev_rule()->raise(eh, _target_state - 1);
1922 
1923  eh = Base::mesh_.EH(Base::mesh_.NHEH(hh2));
1924  Base::prev_rule()->raise(eh, _target_state - 1);
1925 
1926  eh = Base::mesh_.EH(Base::mesh_.PHEH(hh2));
1927  Base::prev_rule()->raise(eh, _target_state - 1);
1928  }
1929 
1930  // calculate new position
1931  typename M::Point position(0.0, 0.0, 0.0);
1932  const typename M::Scalar valence(4.0);
1933 
1934  position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh1))).position(_target_state - 1);
1935  position += MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh1))).position(_target_state - 1);
1936  position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh2))).position(_target_state - 1);
1937  position += MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh2))).position(_target_state - 1);
1938 
1939  position /= valence;
1940 
1941  MOBJ(_eh).set_position(_target_state, position);
1942  MOBJ(_eh).inc_state();
1943  }
1944 }
1945 
1946 
1947 // ------------------------------------------------------------------ EdEc ----
1948 
1949 
1950 template<class M>
1951 void
1952 EdEc<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state)
1953 {
1954  if (MOBJ(_eh).state() < _target_state) {
1955 
1956  this->update(_eh, _target_state);
1957 
1958  // raise all neighbor faces and edges to level x-1
1959  typename M::HalfedgeHandle hh1, hh2;
1960  typename M::FaceHandle fh;
1961  typename M::EdgeHandle eh;
1962 
1963  hh1 = Base::mesh_.HEH(_eh, 0);
1964  hh2 = Base::mesh_.HEH(_eh, 1);
1965 
1966  if (_target_state > 1) {
1967 
1968  fh = Base::mesh_.FH(hh1);
1969  Base::prev_rule()->raise(fh, _target_state - 1);
1970 
1971  fh = Base::mesh_.FH(hh2);
1972  Base::prev_rule()->raise(fh, _target_state - 1);
1973 
1974  eh = Base::mesh_.EH(Base::mesh_.NHEH(hh1));
1975  Base::prev_rule()->raise(eh, _target_state - 1);
1976 
1977  eh = Base::mesh_.EH(Base::mesh_.PHEH(hh1));
1978  Base::prev_rule()->raise(eh, _target_state - 1);
1979 
1980  eh = Base::mesh_.EH(Base::mesh_.NHEH(hh2));
1981  Base::prev_rule()->raise(eh, _target_state - 1);
1982 
1983  eh = Base::mesh_.EH(Base::mesh_.PHEH(hh2));
1984  Base::prev_rule()->raise(eh, _target_state - 1);
1985  }
1986 
1987  // calculate new position
1988  typename M::Point position(0.0, 0.0, 0.0);
1989  const typename M::Scalar valence(4.0);
1990  typename M::Scalar c;
1991 
1992  position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh1))).position(_target_state - 1);
1993  position += MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh1))).position(_target_state - 1);
1994  position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh2))).position(_target_state - 1);
1995  position += MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh2))).position(_target_state - 1);
1996 
1997  position /= valence;
1998 
1999  // choose coefficient c
2000  c = Base::coeff();
2001 
2002  position *= ( static_cast<typename M::Scalar>(1.0) - c);
2003 
2004  position += MOBJ(_eh).position(_target_state - 1) * c;
2005 
2006  MOBJ(_eh).set_position(_target_state, position);
2007  MOBJ(_eh).inc_state();
2008  }
2009 }
2010 
2011 #endif // DOXY_IGNORE_THIS
2012 
2013 #undef FH
2014 #undef VH
2015 #undef EH
2016 #undef HEH
2017 #undef M
2018 #undef MOBJ
2019 
2020 //=============================================================================
2021 } // END_NS_ADAPTIVE
2022 } // END_NS_SUBDIVIDER
2023 } // END_NS_OPENMESH
2024 //=============================================================================
OpenMesh::Subdivider::Adaptive::state_t
CompositeTraits::state_t state_t
Definition: CompositeTraits.hh:250
OpenMesh::Subdivider::Adaptive::Tvv3::raise
void raise(typename M::FaceHandle &_fh, state_t _target_state)
Raise item to target state _target_state.
Definition: RulesT_impl.hh:99
OpenMesh::Subdivider::Adaptive::EVc::raise
void raise(typename M::VertexHandle &_vh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::FV::raise
void raise(typename M::VertexHandle &_vh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::VF::raise
void raise(typename M::FaceHandle &_fh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::VVc::raise
void raise(typename M::VertexHandle &_vh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::FVc::raise
void raise(typename M::VertexHandle &_vh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::VdEc::raise
void raise(typename M::EdgeHandle &_eh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::FE::raise
void raise(typename M::EdgeHandle &_eh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::EF::raise
void raise(typename M::FaceHandle &_fh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::Tvv4::raise
void raise(typename M::FaceHandle &_fh, state_t _target_state)
Raise item to target state _target_state.
Definition: RulesT_impl.hh:338
OpenMesh::Subdivider::Adaptive::EdEc::raise
void raise(typename M::EdgeHandle &_eh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::VV::raise
void raise(typename M::VertexHandle &_vh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::VE::raise
void raise(typename M::EdgeHandle &_eh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::EV::raise
void raise(typename M::VertexHandle &_vh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::EdE::raise
void raise(typename M::EdgeHandle &_eh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::FF::raise
void raise(typename M::FaceHandle &_fh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::VdE::raise
void raise(typename M::EdgeHandle &_eh, state_t _target_state)
Raise item to target state _target_state.
OpenMesh::Subdivider::Adaptive::FFc::raise
void raise(typename M::FaceHandle &_fh, state_t _target_state)
Raise item to target state _target_state.