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