openmesh: RulesT.cc Source File

00001 /*===========================================================================*\
00002  *                                                                           *
00003  *                               OpenMesh                                    *
00004  *      Copyright (C) 2001-2011 by Computer Graphics Group, RWTH Aachen      *
00005  *                           www.openmesh.org                                *
00006  *                                                                           *
00007  *---------------------------------------------------------------------------* 
00008  *  This file is part of OpenMesh.                                           *
00009  *                                                                           *
00010  *  OpenMesh is free software: you can redistribute it and/or modify         * 
00011  *  it under the terms of the GNU Lesser General Public License as           *
00012  *  published by the Free Software Foundation, either version 3 of           *
00013  *  the License, or (at your option) any later version with the              *
00014  *  following exceptions:                                                    *
00015  *                                                                           *
00016  *  If other files instantiate templates or use macros                       *
00017  *  or inline functions from this file, or you compile this file and         *
00018  *  link it with other files to produce an executable, this file does        *
00019  *  not by itself cause the resulting executable to be covered by the        *
00020  *  GNU Lesser General Public License. This exception does not however       *
00021  *  invalidate any other reasons why the executable file might be            *
00022  *  covered by the GNU Lesser General Public License.                        *
00023  *                                                                           *
00024  *  OpenMesh is distributed in the hope that it will be useful,              *
00025  *  but WITHOUT ANY WARRANTY; without even the implied warranty of           *
00026  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            *
00027  *  GNU Lesser General Public License for more details.                      *
00028  *                                                                           *
00029  *  You should have received a copy of the GNU LesserGeneral Public          *
00030  *  License along with OpenMesh.  If not,                                    *
00031  *  see <http://www.gnu.org/licenses/>.                                      *
00032  *                                                                           *
00033 \*===========================================================================*/ 
00034 
00035 /*===========================================================================*\
00036  *                                                                           *             
00037  *   $Revision: 362 $                                                         *
00038  *   $Date: 2011-01-26 10:21:12 +0100 (Mi, 26 Jan 2011) $                   *
00039  *                                                                           *
00040 \*===========================================================================*/
00041 
00046 //=============================================================================
00047 //
00048 //  Rules - IMPLEMENTATION
00049 //
00050 //=============================================================================
00051 
00052 
00053 #define OPENMESH_SUBDIVIDER_ADAPTIVE_RULEST_CC
00054 
00055 
00056 //== INCLUDES =================================================================
00057 
00058 #include <OpenMesh/Core/System/config.h>
00059 #include <OpenMesh/Core/IO/MeshIO.hh>
00060 #include "RulesT.hh"
00061 // --------------------
00062 #if defined(OM_CC_MIPS)
00063 #  include <math.h>
00064 #else
00065 #  include <cmath>
00066 #endif
00067 
00068 #if defined(OM_CC_MSVC)
00069 #  pragma warning(disable:4244)
00070 #endif
00071 
00072 //== NAMESPACE ================================================================
00073 
00074 namespace OpenMesh   { // BEGIN_NS_OPENMESH
00075 namespace Subdivider { // BEGIN_NS_DECIMATER
00076 namespace Adaptive   { // BEGIN_NS_ADAPTIVE
00077 
00078 
00079 //== IMPLEMENTATION ========================================================== 
00080 
00081 #define MOBJ Base::mesh_.data
00082 #define FH face_handle
00083 #define VH vertex_handle
00084 #define EH edge_handle
00085 #define HEH halfedge_handle
00086 #define NHEH next_halfedge_handle
00087 #define PHEH prev_halfedge_handle
00088 #define OHEH opposite_halfedge_handle
00089 #define TVH  to_vertex_handle
00090 #define FVH  from_vertex_handle
00091 
00092 // ------------------------------------------------------------------ Tvv3 ----
00093 
00094 
00095 template<class M>
00096 void
00097 Tvv3<M>::raise(typename M::FaceHandle& _fh, state_t _target_state) 
00098 {
00099   if (MOBJ(_fh).state() < _target_state) 
00100   {
00101     update(_fh, _target_state);
00102 
00103     typename M::VertexVertexIter          vv_it;
00104     typename M::FaceVertexIter            fv_it;
00105     typename M::VertexHandle              vh;
00106     typename M::Point                     position(0.0, 0.0, 0.0);
00107     typename M::Point                     face_position;
00108     const typename M::Point               zero_point(0.0, 0.0, 0.0);
00109     std::vector<typename M::VertexHandle> vertex_vector;
00110     int                                      valence(0);
00111 
00112     // raise all adjacent vertices to level x-1
00113     for (fv_it = Base::mesh_.fv_iter(_fh); fv_it; ++fv_it) {
00114 
00115       vertex_vector.push_back(fv_it.handle());
00116     }
00117 
00118     while(!vertex_vector.empty()) {
00119 
00120       vh = vertex_vector.back();
00121       vertex_vector.pop_back();
00122 
00123       if (_target_state > 1)
00124         Base::prev_rule()->raise(vh, _target_state - 1);
00125     }
00126 
00127     face_position = MOBJ(_fh).position(_target_state - 1);
00128     
00129     typename M::EdgeHandle               eh;
00130     std::vector<typename M::EdgeHandle>  edge_vector;
00131 
00132     // interior face
00133     if (!Base::mesh_.is_boundary(_fh) || MOBJ(_fh).final()) { 
00134 
00135       // insert new vertex
00136       vh = Base::mesh_.new_vertex();
00137 
00138       Base::mesh_.split(_fh, vh);
00139 
00140       // calculate display position for new vertex
00141       for (vv_it = Base::mesh_.vv_iter(vh); vv_it; ++vv_it) 
00142       {
00143         position += Base::mesh_.point(vv_it.handle());
00144         ++valence;
00145       }
00146 
00147       position /= valence;
00148 
00149       // set attributes for new vertex
00150       Base::mesh_.set_point(vh, position);
00151       MOBJ(vh).set_position(_target_state, zero_point);
00152       MOBJ(vh).set_state(_target_state);
00153       MOBJ(vh).set_not_final();
00154 
00155       typename M::VertexOHalfedgeIter      voh_it;
00156       // check for edge flipping
00157       for (voh_it = Base::mesh_.voh_iter(vh); voh_it; ++voh_it) {
00158       
00159         if (Base::mesh_.FH(voh_it.handle()).is_valid()) {
00160 
00161           MOBJ(Base::mesh_.FH(voh_it.handle())).set_state(_target_state);
00162           MOBJ(Base::mesh_.FH(voh_it.handle())).set_not_final();
00163           MOBJ(Base::mesh_.FH(voh_it.handle())).set_position(_target_state - 1, face_position);
00164       
00165 
00166           for (state_t j = 0; j < _target_state; ++j) {
00167             MOBJ(Base::mesh_.FH(voh_it.handle())).set_position(j, MOBJ(_fh).position(j));
00168           }
00169       
00170           if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle()))).is_valid()) {
00171 
00172             if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle())))).state() == _target_state) {
00173 
00174               if (Base::mesh_.is_flip_ok(Base::mesh_.EH(Base::mesh_.NHEH(voh_it.handle())))) {
00175 
00176                 edge_vector.push_back(Base::mesh_.EH(Base::mesh_.NHEH(voh_it.handle())));
00177               }
00178             }
00179           }
00180         }
00181       }
00182     }
00183     
00184     // boundary face
00185     else { 
00186 
00187       typename M::VertexHandle vh1 = Base::mesh_.new_vertex(),
00188                                   vh2 = Base::mesh_.new_vertex();
00189       
00190       typename M::HalfedgeHandle hh2 = Base::mesh_.HEH(_fh),
00191                                     hh1, hh3;
00192       
00193       while (!Base::mesh_.is_boundary(Base::mesh_.OHEH(hh2)))
00194         hh2 = Base::mesh_.NHEH(hh2);
00195       
00196       eh = Base::mesh_.EH(hh2);
00197       
00198       hh2 = Base::mesh_.NHEH(hh2);
00199       hh1 = Base::mesh_.NHEH(hh2);
00200       
00201       assert(Base::mesh_.is_boundary(eh));
00202 
00203       Base::mesh_.split(eh, vh1);
00204 
00205       eh = Base::mesh_.EH(Base::mesh_.PHEH(hh2));
00206 
00207       assert(Base::mesh_.is_boundary(eh));
00208 
00209       Base::mesh_.split(eh, vh2);
00210 
00211       hh3 = Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(hh1)));
00212 
00213       typename M::VertexHandle   vh0(Base::mesh_.TVH(hh1)),
00214                                     vh3(Base::mesh_.FVH(hh2));
00215 
00216       // set display position and attributes for new vertices
00217       Base::mesh_.set_point(vh1, (Base::mesh_.point(vh0) * 2.0 + Base::mesh_.point(vh3)) / 3.0);
00218 
00219       MOBJ(vh1).set_position(_target_state, zero_point);
00220       MOBJ(vh1).set_state(_target_state);
00221       MOBJ(vh1).set_not_final();
00222 
00223       MOBJ(vh0).set_position(_target_state, MOBJ(vh0).position(_target_state - 1) * 3.0);
00224       MOBJ(vh0).set_state(_target_state);
00225       MOBJ(vh0).set_not_final();
00226 
00227       // set display position and attributes for old vertices
00228       Base::mesh_.set_point(vh2, (Base::mesh_.point(vh3) * 2.0 + Base::mesh_.point(vh0)) / 3.0);
00229       MOBJ(vh2).set_position(_target_state, zero_point);
00230       MOBJ(vh2).set_state(_target_state);
00231       MOBJ(vh2).set_not_final();
00232 
00233       MOBJ(vh3).set_position(_target_state, MOBJ(vh3).position(_target_state - 1) * 3.0);
00234       MOBJ(vh3).set_state(_target_state);
00235       MOBJ(vh3).set_not_final();
00236 
00237       // init 3 faces
00238       MOBJ(Base::mesh_.FH(hh1)).set_state(_target_state);
00239       MOBJ(Base::mesh_.FH(hh1)).set_not_final();
00240       MOBJ(Base::mesh_.FH(hh1)).set_position(_target_state - 1, face_position);
00241       
00242       MOBJ(Base::mesh_.FH(hh2)).set_state(_target_state);
00243       MOBJ(Base::mesh_.FH(hh2)).set_not_final();
00244       MOBJ(Base::mesh_.FH(hh2)).set_position(_target_state - 1, face_position);
00245 
00246       MOBJ(Base::mesh_.FH(hh3)).set_state(_target_state);
00247       MOBJ(Base::mesh_.FH(hh3)).set_final();
00248       MOBJ(Base::mesh_.FH(hh3)).set_position(_target_state - 1, face_position);
00249       
00250 
00251       for (state_t j = 0; j < _target_state; ++j) {
00252         MOBJ(Base::mesh_.FH(hh1)).set_position(j, MOBJ(_fh).position(j));
00253       }
00254       
00255       for (state_t j = 0; j < _target_state; ++j) {
00256 
00257         MOBJ(Base::mesh_.FH(hh2)).set_position(j, MOBJ(_fh).position(j));
00258       }
00259       
00260       for (state_t j = 0; j < _target_state; ++j) {
00261 
00262         MOBJ(Base::mesh_.FH(hh3)).set_position(j, MOBJ(_fh).position(j));
00263       }
00264       
00265       // check for edge flipping
00266       if (Base::mesh_.FH(Base::mesh_.OHEH(hh1)).is_valid()) {
00267 
00268         if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(hh1))).state() == _target_state) {
00269 
00270           if (Base::mesh_.is_flip_ok(Base::mesh_.EH(hh1))) {
00271 
00272             edge_vector.push_back(Base::mesh_.EH(hh1));
00273           }
00274         }
00275       }
00276 
00277       if (Base::mesh_.FH(Base::mesh_.OHEH(hh2)).is_valid()) {
00278 
00279         if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(hh2))).state() == _target_state) {
00280 
00281           if (Base::mesh_.is_flip_ok(Base::mesh_.EH(hh2))) {
00282 
00283             edge_vector.push_back(Base::mesh_.EH(hh2));
00284           }
00285         }
00286       }
00287     }
00288   
00289     // flip edges
00290     while (!edge_vector.empty()) {
00291 
00292       eh = edge_vector.back();
00293       edge_vector.pop_back();
00294       
00295       assert(Base::mesh_.is_flip_ok(eh));
00296 
00297       Base::mesh_.flip(eh);
00298 
00299       MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 0))).set_final();
00300       MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 1))).set_final();
00301 
00302       MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 0))).set_state(_target_state);
00303       MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 1))).set_state(_target_state);
00304 
00305       MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 0))).set_position(_target_state, face_position);
00306       MOBJ(Base::mesh_.FH(Base::mesh_.HEH(eh, 1))).set_position(_target_state, face_position);
00307     }
00308   }
00309 }
00310 
00311 
00312 template<class M>
00313 void Tvv3<M>::raise(typename M::VertexHandle& _vh, state_t _target_state) {
00314 
00315   if (MOBJ(_vh).state() < _target_state) {
00316 
00317     update(_vh, _target_state);
00318 
00319     // multiply old position by 3
00320     MOBJ(_vh).set_position(_target_state, MOBJ(_vh).position(_target_state - 1) * 3.0);
00321 
00322     MOBJ(_vh).inc_state();
00323 
00324     assert(MOBJ(_vh).state() == _target_state);
00325   }
00326 }
00327 
00328 
00329 // ------------------------------------------------------------------ Tvv4 ----
00330 
00331 
00332 template<class M>
00333 void
00334 Tvv4<M>::raise(typename M::FaceHandle& _fh, state_t _target_state) 
00335 {
00336 
00337   if (MOBJ(_fh).state() < _target_state) {
00338 
00339     update(_fh, _target_state);
00340 
00341     typename M::FaceVertexIter              fv_it;
00342     typename M::VertexHandle                temp_vh;
00343     typename M::Point                       face_position;
00344     const typename M::Point                 zero_point(0.0, 0.0, 0.0);
00345     std::vector<typename M::VertexHandle>   vertex_vector;
00346     std::vector<typename M::HalfedgeHandle> halfedge_vector;
00347 
00348     // raise all adjacent vertices to level x-1
00349     for (fv_it = Base::mesh_.fv_iter(_fh); fv_it; ++fv_it) {
00350 
00351       vertex_vector.push_back(fv_it.handle());
00352     }
00353 
00354     while(!vertex_vector.empty()) {
00355 
00356       temp_vh = vertex_vector.back();
00357       vertex_vector.pop_back();
00358 
00359       if (_target_state > 1) {
00360         Base::prev_rule()->raise(temp_vh, _target_state - 1);
00361       }
00362     }
00363 
00364     face_position = MOBJ(_fh).position(_target_state - 1);
00365 
00366     typename M::HalfedgeHandle hh[3];
00367     typename M::VertexHandle   vh[3];
00368     typename M::VertexHandle   new_vh[3];
00369     typename M::FaceHandle     fh[4];
00370     typename M::EdgeHandle     eh;
00371     typename M::HalfedgeHandle temp_hh;
00372 
00373     // normal (final) face
00374     if (MOBJ(_fh).final()) {
00375 
00376       // define three halfedge handles around the face
00377       hh[0] = Base::mesh_.HEH(_fh);
00378       hh[1] = Base::mesh_.NHEH(hh[0]);
00379       hh[2] = Base::mesh_.NHEH(hh[1]);
00380 
00381       assert(hh[0] == Base::mesh_.NHEH(hh[2]));
00382 
00383       vh[0] = Base::mesh_.TVH(hh[0]);
00384       vh[1] = Base::mesh_.TVH(hh[1]);
00385       vh[2] = Base::mesh_.TVH(hh[2]);
00386       
00387       new_vh[0] = Base::mesh_.add_vertex(zero_point);
00388       new_vh[1] = Base::mesh_.add_vertex(zero_point);
00389       new_vh[2] = Base::mesh_.add_vertex(zero_point);
00390 
00391       // split three edges
00392       split_edge(hh[0], new_vh[0], _target_state);
00393       eh = Base::mesh_.EH(Base::mesh_.PHEH(hh[2]));
00394       split_edge(hh[1], new_vh[1], _target_state);
00395       split_edge(hh[2], new_vh[2], _target_state);
00396 
00397       assert(Base::mesh_.FVH(hh[2]) == vh[1]);
00398       assert(Base::mesh_.FVH(hh[1]) == vh[0]);
00399       assert(Base::mesh_.FVH(hh[0]) == vh[2]);
00400 
00401       if (Base::mesh_.FH(Base::mesh_.OHEH(hh[0])).is_valid()) 
00402       {
00403         temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[0])));
00404         if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
00405           halfedge_vector.push_back(temp_hh);
00406       }
00407 
00408       if (Base::mesh_.FH(Base::mesh_.OHEH(hh[1])).is_valid()) {
00409 
00410         temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[1])));
00411         if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
00412           halfedge_vector.push_back(temp_hh);
00413       }
00414 
00415       if (Base::mesh_.FH(Base::mesh_.OHEH(hh[2])).is_valid()) {
00416 
00417         temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[2])));
00418         if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
00419           halfedge_vector.push_back(temp_hh);
00420       }
00421     }
00422     
00423     // splitted face, check for type
00424     else {
00425 
00426       // define red halfedge handle
00427       typename M::HalfedgeHandle red_hh(MOBJ(_fh).red_halfedge());
00428 
00429       if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh))).is_valid() 
00430           && Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)))).is_valid() 
00431           && MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)))).red_halfedge() == red_hh 
00432           && MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh))))).red_halfedge() == red_hh) 
00433       {
00434 
00435         // three times divided face
00436         vh[0] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)))));
00437         vh[1] = Base::mesh_.TVH(red_hh);
00438         vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh))));
00439           
00440         new_vh[0] = Base::mesh_.FVH(red_hh);
00441         new_vh[1] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)));
00442         new_vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(red_hh));
00443     
00444         hh[0] = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)));
00445         hh[1] = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh))));
00446         hh[2] = Base::mesh_.NHEH(red_hh);
00447     
00448         eh = Base::mesh_.EH(red_hh);
00449       }
00450       
00451       else 
00452       {
00453 
00454         if ((Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh))).is_valid() && 
00455              MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)))).red_halfedge() 
00456              == red_hh )
00457             || (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)))).is_valid() 
00458             && MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh))))).red_halfedge() == red_hh))
00459         {
00460 
00461           // double divided face
00462           if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)))).red_halfedge() == red_hh) 
00463           {
00464         // first case
00465         vh[0] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)));
00466         vh[1] = Base::mesh_.TVH(red_hh);
00467         vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh))));
00468       
00469         new_vh[0] = Base::mesh_.FVH(red_hh);
00470         new_vh[1] = Base::mesh_.add_vertex(zero_point);
00471         new_vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(red_hh));
00472   
00473         hh[0] = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(red_hh)));
00474         hh[1] = Base::mesh_.PHEH(Base::mesh_.OHEH(red_hh));
00475         hh[2] = Base::mesh_.NHEH(red_hh);
00476   
00477         // split one edge
00478         eh = Base::mesh_.EH(red_hh);
00479   
00480         split_edge(hh[1], new_vh[1], _target_state);
00481   
00482         assert(Base::mesh_.FVH(hh[2]) == vh[1]);
00483         assert(Base::mesh_.FVH(hh[1]) == vh[0]);
00484         assert(Base::mesh_.FVH(hh[0]) == vh[2]);
00485         
00486         if (Base::mesh_.FH(Base::mesh_.OHEH(hh[1])).is_valid()) 
00487         {
00488           temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[1])));
00489           if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
00490             halfedge_vector.push_back(temp_hh);
00491         }
00492           }
00493           else 
00494           {
00495 
00496             // second case
00497             vh[0] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)))));
00498             vh[1] = Base::mesh_.TVH(red_hh);
00499             vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(red_hh));
00500       
00501             new_vh[0] = Base::mesh_.FVH(red_hh);
00502             new_vh[1] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)));
00503             new_vh[2] = Base::mesh_.add_vertex(zero_point);
00504 
00505             hh[0] = Base::mesh_.PHEH(red_hh);
00506             hh[1] = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh))));
00507             hh[2] = Base::mesh_.NHEH(red_hh);
00508 
00509             // split one edge
00510             eh = Base::mesh_.EH(red_hh);
00511 
00512             split_edge(hh[2], new_vh[2], _target_state);
00513 
00514             assert(Base::mesh_.FVH(hh[2]) == vh[1]);
00515             assert(Base::mesh_.FVH(hh[1]) == vh[0]);
00516             assert(Base::mesh_.FVH(hh[0]) == vh[2]);
00517 
00518             if (Base::mesh_.FH(Base::mesh_.OHEH(hh[2])).is_valid()) {
00519             
00520               temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[2])));
00521               if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
00522                 halfedge_vector.push_back(temp_hh);
00523             }
00524           }
00525         }
00526 
00527         else {
00528           
00529           // one time divided face
00530           vh[0] = Base::mesh_.TVH(Base::mesh_.NHEH(Base::mesh_.OHEH(red_hh)));
00531           vh[1] = Base::mesh_.TVH(red_hh);
00532           vh[2] = Base::mesh_.TVH(Base::mesh_.NHEH(red_hh));
00533       
00534           new_vh[0] = Base::mesh_.FVH(red_hh);
00535           new_vh[1] = Base::mesh_.add_vertex(zero_point);
00536           new_vh[2] = Base::mesh_.add_vertex(zero_point);
00537 
00538           hh[0] = Base::mesh_.PHEH(red_hh);
00539           hh[1] = Base::mesh_.PHEH(Base::mesh_.OHEH(red_hh));
00540           hh[2] = Base::mesh_.NHEH(red_hh);
00541 
00542           // split two edges
00543           eh = Base::mesh_.EH(red_hh);
00544 
00545           split_edge(hh[1], new_vh[1], _target_state);
00546           split_edge(hh[2], new_vh[2], _target_state);
00547 
00548           assert(Base::mesh_.FVH(hh[2]) == vh[1]);
00549           assert(Base::mesh_.FVH(hh[1]) == vh[0]);
00550           assert(Base::mesh_.FVH(hh[0]) == vh[2]);
00551 
00552           if (Base::mesh_.FH(Base::mesh_.OHEH(hh[1])).is_valid()) {
00553 
00554             temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[1])));
00555             if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
00556               halfedge_vector.push_back(temp_hh);
00557           }
00558 
00559           if (Base::mesh_.FH(Base::mesh_.OHEH(hh[2])).is_valid()) {
00560 
00561             temp_hh = Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(hh[2])));
00562             if (MOBJ(Base::mesh_.FH(temp_hh)).red_halfedge() != temp_hh)
00563               halfedge_vector.push_back(temp_hh);
00564           }
00565         }
00566       }
00567     }
00568 
00569     // continue here for all cases
00570         
00571     // flip edge
00572     if (Base::mesh_.is_flip_ok(eh)) {
00573       
00574       Base::mesh_.flip(eh);
00575     }
00576     
00577     // search new faces
00578     fh[0] = Base::mesh_.FH(hh[0]);
00579     fh[1] = Base::mesh_.FH(hh[1]);
00580     fh[2] = Base::mesh_.FH(hh[2]);
00581     fh[3] = Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(hh[0])));
00582 
00583     assert(_fh == fh[0] || _fh == fh[1] || _fh == fh[2] || _fh == fh[3]); 
00584 
00585     // init new faces
00586     for (int i = 0; i <= 3; ++i) {
00587 
00588       MOBJ(fh[i]).set_state(_target_state);
00589       MOBJ(fh[i]).set_final();
00590       MOBJ(fh[i]).set_position(_target_state, face_position);
00591       MOBJ(fh[i]).set_red_halfedge(Base::mesh_.InvalidHalfedgeHandle);
00592     }
00593 
00594     // init new vertices and edges
00595     for (int i = 0; i <= 2; ++i) {
00596         
00597       MOBJ(new_vh[i]).set_position(_target_state, zero_point);
00598       MOBJ(new_vh[i]).set_state(_target_state);
00599       MOBJ(new_vh[i]).set_not_final();
00600 
00601       Base::mesh_.set_point(new_vh[i], (Base::mesh_.point(vh[i]) + Base::mesh_.point(vh[(i + 2) % 3])) * 0.5);
00602 
00603       MOBJ(Base::mesh_.EH(hh[i])).set_state(_target_state);
00604       MOBJ(Base::mesh_.EH(hh[i])).set_position(_target_state, zero_point);
00605       MOBJ(Base::mesh_.EH(hh[i])).set_final();
00606 
00607       MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh[i]))).set_state(_target_state);
00608       MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh[i]))).set_position(_target_state, zero_point);
00609       MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh[i]))).set_final();
00610 
00611       MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh[i]))).set_state(_target_state);
00612       MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh[i]))).set_position(_target_state, zero_point);
00613       MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh[i]))).set_final();
00614     }
00615 
00616     // check, if opposite triangle needs splitting
00617     while (!halfedge_vector.empty()) {
00618 
00619       temp_hh = halfedge_vector.back();
00620       halfedge_vector.pop_back();
00621 
00622       check_edge(temp_hh, _target_state);
00623     }
00624 
00625     assert(MOBJ(fh[0]).state() == _target_state);
00626     assert(MOBJ(fh[1]).state() == _target_state);
00627     assert(MOBJ(fh[2]).state() == _target_state);
00628     assert(MOBJ(fh[3]).state() == _target_state);
00629   }
00630 }
00631 
00632 
00633 template<class M>
00634 void
00635 Tvv4<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
00636 {
00637  
00638   if (MOBJ(_vh).state() < _target_state)
00639   {
00640 
00641     update(_vh, _target_state);
00642 
00643     // multiply old position by 4
00644     MOBJ(_vh).set_position(_target_state, MOBJ(_vh).position(_target_state - 1) * 4.0);
00645 
00646     MOBJ(_vh).inc_state();
00647   }
00648 }
00649 
00650 
00651 template<class M>
00652 void
00653 Tvv4<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) 
00654 { 
00655   if (MOBJ(_eh).state() < _target_state) 
00656   {
00657     update(_eh, _target_state);
00658 
00659     typename M::FaceHandle fh(Base::mesh_.FH(Base::mesh_.HEH(_eh, 0)));
00660 
00661     if (!fh.is_valid())
00662       fh=Base::mesh_.FH(Base::mesh_.HEH(_eh, 1));
00663 
00664     raise(fh, _target_state);
00665 
00666     assert(MOBJ(_eh).state() == _target_state);
00667   }
00668 }
00669 
00670 #ifndef DOXY_IGNORE_THIS
00671 
00672 template<class M>
00673 void
00674 Tvv4<M>::split_edge(typename M::HalfedgeHandle &_hh, 
00675                        typename M::VertexHandle   &_vh, 
00676                        state_t _target_state) 
00677 {
00678   typename M::HalfedgeHandle temp_hh;
00679 
00680   if (Base::mesh_.FH(Base::mesh_.OHEH(_hh)).is_valid()) 
00681   {
00682     if (!MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).final()) 
00683     {    
00684       if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).red_halfedge().is_valid()) 
00685       {        
00686         temp_hh = MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).red_halfedge();
00687       }
00688       else 
00689       {
00690         // two cases for divided, but not visited face
00691         if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh))))).state() 
00692             == MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).state()) 
00693         {
00694           temp_hh = Base::mesh_.PHEH(Base::mesh_.OHEH(_hh));
00695         }
00696 
00697         else if (MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(_hh))))).state() 
00698                  == MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).state())
00699         {
00700           temp_hh = Base::mesh_.NHEH(Base::mesh_.OHEH(_hh)); 
00701         }
00702       }
00703     }
00704     else
00705       temp_hh = Base::mesh_.InvalidHalfedgeHandle;
00706   }
00707   else
00708     temp_hh = Base::mesh_.InvalidHalfedgeHandle;
00709          
00710   // split edge
00711   Base::mesh_.split(Base::mesh_.EH(_hh), _vh);
00712           
00713   if (Base::mesh_.FVH(_hh) == _vh) 
00714   {         
00715     MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(_hh))))).set_state(MOBJ(Base::mesh_.EH(_hh)).state());
00716     _hh = Base::mesh_.PHEH(Base::mesh_.OHEH(Base::mesh_.PHEH(_hh)));
00717   }
00718 
00719   if (Base::mesh_.FH(Base::mesh_.OHEH(_hh)).is_valid()) {
00720           
00721     MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh)))).set_not_final();
00722     MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).set_state(_target_state-1);
00723     MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh)))))).set_state(_target_state-1);
00724 
00725     MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).set_not_final();
00726     MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh)))))).set_not_final();
00727 
00728     MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh)))).set_state(_target_state);
00729 
00730     if (temp_hh.is_valid()) {
00731         
00732       MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(_hh))).set_red_halfedge(temp_hh);
00733       MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh)))))).set_red_halfedge(temp_hh);
00734     } 
00735     else {
00736 
00737       typename M::FaceHandle 
00738         fh1(Base::mesh_.FH(Base::mesh_.OHEH(_hh))),
00739         fh2(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh))))));
00740 
00741       MOBJ(fh1).set_red_halfedge(Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh))));
00742       MOBJ(fh2).set_red_halfedge(Base::mesh_.OHEH(Base::mesh_.PHEH(Base::mesh_.OHEH(_hh))));
00743 
00744       const typename M::Point zero_point(0.0, 0.0, 0.0);
00745 
00746       MOBJ(fh1).set_position(_target_state - 1, zero_point);
00747       MOBJ(fh2).set_position(_target_state - 1, zero_point);
00748     }
00749   }
00750 
00751   // init edges
00752   MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh))))).set_state(_target_state - 1);
00753   MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(Base::mesh_.OHEH(Base::mesh_.NHEH(_hh))))).set_final();
00754 
00755   MOBJ(Base::mesh_.EH(_hh)).set_state(_target_state - 1);
00756   MOBJ(Base::mesh_.EH(_hh)).set_final();
00757 }  
00758 
00759 
00760 template<class M>
00761 void Tvv4<M>::check_edge(const typename M::HalfedgeHandle& _hh, 
00762                          state_t _target_state) 
00763 {
00764   typename M::FaceHandle fh1(Base::mesh_.FH(_hh)),
00765     fh2(Base::mesh_.FH(Base::mesh_.OHEH(_hh)));
00766 
00767   assert(fh1.is_valid());
00768   assert(fh2.is_valid());
00769 
00770   typename M::HalfedgeHandle red_hh(MOBJ(fh1).red_halfedge());
00771   
00772   if (!MOBJ(fh1).final()) {
00773 
00774     assert (MOBJ(fh1).final() == MOBJ(fh2).final());
00775     assert (!MOBJ(fh1).final());
00776     assert (MOBJ(fh1).red_halfedge() == MOBJ(fh2).red_halfedge());
00777 
00778     const typename M::Point zero_point(0.0, 0.0, 0.0);
00779 
00780     MOBJ(fh1).set_position(_target_state - 1, zero_point);
00781     MOBJ(fh2).set_position(_target_state - 1, zero_point);
00782 
00783     assert(red_hh.is_valid());
00784 
00785     if (!red_hh.is_valid()) {
00786 
00787       MOBJ(fh1).set_state(_target_state - 1);
00788       MOBJ(fh2).set_state(_target_state - 1);
00789 
00790       MOBJ(fh1).set_red_halfedge(_hh);
00791       MOBJ(fh2).set_red_halfedge(_hh);
00792 
00793       MOBJ(Base::mesh_.EH(_hh)).set_not_final();
00794       MOBJ(Base::mesh_.EH(_hh)).set_state(_target_state - 1);
00795     }
00796 
00797     else {
00798 
00799       MOBJ(Base::mesh_.EH(_hh)).set_not_final();
00800       MOBJ(Base::mesh_.EH(_hh)).set_state(_target_state - 1);
00801 
00802       raise(fh1, _target_state);
00803 
00804       assert(MOBJ(fh1).state() == _target_state);
00805     }
00806   }
00807 }
00808 
00809 
00810 // -------------------------------------------------------------------- VF ----
00811 
00812 
00813 template<class M>
00814 void VF<M>::raise(typename M::FaceHandle& _fh, state_t _target_state) 
00815 {
00816   if (MOBJ(_fh).state() < _target_state) {
00817 
00818     update(_fh, _target_state);
00819 
00820     // raise all neighbour vertices to level x-1
00821     typename M::FaceVertexIter            fv_it;
00822     typename M::VertexHandle              vh;
00823     std::vector<typename M::VertexHandle> vertex_vector;
00824 
00825     if (_target_state > 1) {
00826 
00827       for (fv_it = Base::mesh_.fv_iter(_fh); fv_it; ++fv_it) {
00828 
00829         vertex_vector.push_back(fv_it.handle());
00830       }
00831 
00832       while (!vertex_vector.empty()) {
00833 
00834         vh = vertex_vector.back();
00835         vertex_vector.pop_back();
00836 
00837         Base::prev_rule()->raise(vh, _target_state - 1);
00838       }
00839     }
00840 
00841     // calculate new position
00842     typename M::Point position(0.0, 0.0, 0.0);
00843     int                  valence(0);
00844 
00845     for (fv_it = Base::mesh_.fv_iter(_fh); fv_it; ++fv_it) {
00846 
00847       ++valence;
00848       position += Base::mesh_.data(fv_it).position(_target_state - 1);
00849     }
00850 
00851     position /= valence;
00852 
00853     // boundary rule
00854     if (Base::number() == Base::subdiv_rule()->number() + 1 && 
00855         Base::mesh_.is_boundary(_fh)                  && 
00856         !MOBJ(_fh).final())
00857       position *= 0.5;
00858 
00859     MOBJ(_fh).set_position(_target_state, position);
00860     MOBJ(_fh).inc_state();
00861 
00862     assert(_target_state == MOBJ(_fh).state());
00863   }
00864 }
00865 
00866 
00867 // -------------------------------------------------------------------- FF ----
00868 
00869 
00870 template<class M>
00871 void FF<M>::raise(typename M::FaceHandle& _fh, state_t _target_state) {
00872 
00873   if (MOBJ(_fh).state() < _target_state) {
00874 
00875     update(_fh, _target_state);
00876 
00877     // raise all neighbour faces to level x-1
00878     typename M::FaceFaceIter              ff_it;
00879     typename M::FaceHandle                fh;
00880     std::vector<typename M::FaceHandle>   face_vector;
00881 
00882     if (_target_state > 1) {
00883 
00884       for (ff_it = Base::mesh_.ff_iter(_fh); ff_it; ++ff_it) {
00885 
00886         face_vector.push_back(ff_it.handle());
00887       }
00888 
00889       while (!face_vector.empty()) {
00890 
00891         fh = face_vector.back();
00892         face_vector.pop_back();
00893 
00894         Base::prev_rule()->raise(fh, _target_state - 1);
00895       }
00896 
00897       for (ff_it = Base::mesh_.ff_iter(_fh); ff_it; ++ff_it) {
00898 
00899         face_vector.push_back(ff_it.handle());
00900       }
00901 
00902       while (!face_vector.empty()) {
00903 
00904         fh = face_vector.back();
00905         face_vector.pop_back();
00906 
00907         while (MOBJ(fh).state() < _target_state - 1)
00908           Base::prev_rule()->raise(fh, _target_state - 1);
00909       }
00910     }
00911 
00912     // calculate new position
00913     typename M::Point position(0.0, 0.0, 0.0);
00914     int                  valence(0);
00915 
00916     for (ff_it = Base::mesh_.ff_iter(_fh); ff_it; ++ff_it) {
00917 
00918       ++valence;
00919 
00920       position += Base::mesh_.data(ff_it).position(_target_state - 1);
00921     }
00922 
00923     position /= valence;
00924 
00925     MOBJ(_fh).set_position(_target_state, position);
00926     MOBJ(_fh).inc_state();
00927   }
00928 }
00929 
00930 
00931 // ------------------------------------------------------------------- FFc ----
00932 
00933 
00934 template<class M>
00935 void FFc<M>::raise(typename M::FaceHandle& _fh, state_t _target_state) 
00936 {
00937   if (MOBJ(_fh).state() < _target_state) {
00938 
00939     update(_fh, _target_state);
00940 
00941     // raise all neighbour faces to level x-1
00942     typename M::FaceFaceIter              ff_it(Base::mesh_.ff_iter(_fh));
00943     typename M::FaceHandle                fh;
00944     std::vector<typename M::FaceHandle>   face_vector;
00945 
00946     if (_target_state > 1) 
00947     {
00948       for (; ff_it; ++ff_it) 
00949         face_vector.push_back(ff_it.handle());
00950 
00951       while (!face_vector.empty()) 
00952       {
00953         fh = face_vector.back();
00954         face_vector.pop_back();
00955         Base::prev_rule()->raise(fh, _target_state - 1);
00956       }
00957 
00958       for (ff_it = Base::mesh_.ff_iter(_fh); ff_it; ++ff_it) 
00959         face_vector.push_back(ff_it.handle());
00960 
00961       while (!face_vector.empty()) {
00962 
00963         fh = face_vector.back();
00964         face_vector.pop_back();
00965 
00966         while (MOBJ(fh).state() < _target_state - 1)
00967           Base::prev_rule()->raise(fh, _target_state - 1);
00968       }
00969     }
00970 
00971     // calculate new position
00972     typename M::Point position(0.0, 0.0, 0.0);
00973     int                  valence(0);
00974 
00975     for (ff_it = Base::mesh_.ff_iter(_fh); ff_it; ++ff_it) 
00976     {
00977       ++valence;
00978       position += Base::mesh_.data(ff_it).position(_target_state - 1);
00979     }
00980 
00981     position /= valence;
00982 
00983     // choose coefficient c
00984     typename M::Scalar c = Base::coeff();
00985 
00986     position *= (1.0 - c);
00987     position += MOBJ(_fh).position(_target_state - 1) * c;
00988 
00989     MOBJ(_fh).set_position(_target_state, position);
00990     MOBJ(_fh).inc_state();
00991   }
00992 }
00993 
00994 
00995 // -------------------------------------------------------------------- FV ----
00996 
00997 
00998 template<class M>
00999 void FV<M>::raise(typename M::VertexHandle& _vh, state_t _target_state) 
01000 {
01001 
01002   if (MOBJ(_vh).state() < _target_state) {
01003 
01004     update(_vh, _target_state);
01005 
01006     // raise all neighbour vertices to level x-1
01007     typename M::VertexFaceIter            vf_it(Base::mesh_.vf_iter(_vh));
01008     typename M::FaceHandle                fh;
01009     std::vector<typename M::FaceHandle>   face_vector;
01010 
01011     if (_target_state > 1) {
01012 
01013       for (; vf_it; ++vf_it) {
01014 
01015         face_vector.push_back(vf_it.handle());
01016       }
01017 
01018       while (!face_vector.empty()) {
01019 
01020         fh = face_vector.back();
01021         face_vector.pop_back();
01022 
01023         Base::prev_rule()->raise(fh, _target_state - 1);
01024       }
01025 
01026       for (vf_it = Base::mesh_.vf_iter(_vh); vf_it; ++vf_it) {
01027 
01028         face_vector.push_back(vf_it.handle());
01029       }
01030 
01031       while (!face_vector.empty()) {
01032 
01033         fh = face_vector.back();
01034         face_vector.pop_back();
01035 
01036         while (MOBJ(fh).state() < _target_state - 1)
01037           Base::prev_rule()->raise(fh, _target_state - 1);
01038       }
01039     }
01040 
01041     // calculate new position
01042     typename M::Point position(0.0, 0.0, 0.0);
01043     int                  valence(0);
01044 
01045     for (vf_it = Base::mesh_.vf_iter(_vh); vf_it; ++vf_it) {
01046 
01047       ++valence;
01048       position += Base::mesh_.data(vf_it).position(_target_state - 1);
01049     }
01050 
01051     position /= valence;
01052 
01053     MOBJ(_vh).set_position(_target_state, position);
01054     MOBJ(_vh).inc_state();
01055 
01056     if (Base::number() == Base::n_rules() - 1) {
01057 
01058       Base::mesh_.set_point(_vh, position);
01059       MOBJ(_vh).set_final();
01060     }
01061   }
01062 }
01063 
01064 
01065 // ------------------------------------------------------------------- FVc ----
01066 
01067 
01068 template<class M>
01069 void FVc<M>::raise(typename M::VertexHandle& _vh, state_t _target_state) 
01070 {
01071   if (MOBJ(_vh).state() < _target_state) {
01072 
01073     update(_vh, _target_state);
01074 
01075     typename M::VertexOHalfedgeIter       voh_it;
01076     typename M::FaceHandle                fh;
01077     std::vector<typename M::FaceHandle>   face_vector;
01078     int                                      valence(0);
01079 
01080     face_vector.clear();
01081 
01082     // raise all neighbour faces to level x-1
01083     if (_target_state > 1) {
01084 
01085       for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it) {
01086 
01087         if (Base::mesh_.FH(voh_it.handle()).is_valid()) {
01088 
01089           face_vector.push_back(Base::mesh_.FH(voh_it.handle()));
01090 
01091           if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle()))).is_valid()) {
01092 
01093             face_vector.push_back(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle()))));
01094           }
01095         }
01096       }
01097 
01098       while (!face_vector.empty()) {
01099 
01100         fh = face_vector.back();
01101         face_vector.pop_back();
01102 
01103         Base::prev_rule()->raise(fh, _target_state - 1);
01104       }
01105 
01106       for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it) {
01107 
01108         if (Base::mesh_.FH(voh_it.handle()).is_valid()) {
01109 
01110           face_vector.push_back(Base::mesh_.FH(voh_it.handle()));
01111 
01112           if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle()))).is_valid()) {
01113 
01114             face_vector.push_back(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle()))));
01115           }
01116         }
01117       }
01118 
01119       while (!face_vector.empty()) {
01120 
01121         fh = face_vector.back();
01122         face_vector.pop_back();
01123 
01124         while (MOBJ(fh).state() < _target_state - 1)
01125           Base::prev_rule()->raise(fh, _target_state - 1);
01126       }
01127 
01128       for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it) {
01129 
01130         if (Base::mesh_.FH(voh_it.handle()).is_valid()) {
01131 
01132           face_vector.push_back(Base::mesh_.FH(voh_it.handle()));
01133 
01134           if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle()))).is_valid()) {
01135 
01136             face_vector.push_back(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle()))));
01137           }
01138         }
01139       }
01140 
01141       while (!face_vector.empty()) {
01142 
01143         fh = face_vector.back();
01144         face_vector.pop_back();
01145 
01146         while (MOBJ(fh).state() < _target_state - 1)
01147           Base::prev_rule()->raise(fh, _target_state - 1);
01148       }
01149     }
01150 
01151     // calculate new position
01152     typename M::Point               position(0.0, 0.0, 0.0);
01153     typename M::Scalar              c;
01154 #if 0
01155     const typename M::Scalar        _2pi(2.0*M_PI);
01156     const typename M::Scalar        _2over3(2.0/3.0);
01157 
01158     for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it) 
01159     {     
01160       ++valence;
01161     }
01162 
01163     // choose coefficient c
01164     c = _2over3 * ( cos( _2pi / valence) + 1.0);
01165 #else
01166     valence = Base::mesh_.valence(_vh);
01167     c       = coeff(valence);
01168 #endif
01169 
01170 
01171     for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it) {
01172 
01173       fh = Base::mesh_.FH(voh_it.handle());
01174       if (fh.is_valid())
01175         Base::prev_rule()->raise(fh, _target_state - 1);
01176 
01177       fh = Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle())));
01178       if (fh.is_valid())
01179         Base::prev_rule()->raise(fh, _target_state - 1);
01180 
01181       if (Base::mesh_.FH(voh_it.handle()).is_valid()) {
01182 
01183         if (Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle()))).is_valid()) {
01184 
01185           position += MOBJ(Base::mesh_.FH(voh_it.handle())).position(_target_state - 1) * c;
01186 
01187           position += MOBJ(Base::mesh_.FH(Base::mesh_.OHEH(Base::mesh_.NHEH(voh_it.handle())))).position(_target_state - 1) * (1.0 - c);
01188         }
01189         else {
01190 
01191           position += MOBJ(Base::mesh_.FH(voh_it.handle())).position(_target_state - 1);
01192         }
01193       }
01194 
01195       else {
01196 
01197         --valence;
01198       }
01199     } 
01200     
01201     position /= valence;
01202 
01203     MOBJ(_vh).set_position(_target_state, position);
01204     MOBJ(_vh).inc_state();
01205 
01206     assert(MOBJ(_vh).state() == _target_state);
01207     
01208     // check if last rule
01209     if (Base::number() == Base::n_rules() - 1) {
01210 
01211       Base::mesh_.set_point(_vh, position);
01212       MOBJ(_vh).set_final();
01213     }
01214   }
01215 }
01216 
01217 template<class M>
01218 std::vector<double> FVc<M>::coeffs_;
01219 
01220 template <class M>
01221 void FVc<M>::init_coeffs(size_t _max_valence)
01222 {
01223   if ( coeffs_.size() == _max_valence+1 )
01224     return;
01225 
01226   if ( coeffs_.size() < _max_valence+1 )
01227   {
01228     const double _2pi(2.0*M_PI);
01229     const double _2over3(2.0/3.0);
01230 
01231     if (coeffs_.empty())
01232       coeffs_.push_back(0.0); // dummy for valence 0
01233       
01234     for(size_t v=coeffs_.size(); v <= _max_valence; ++v)
01235       coeffs_.push_back(_2over3 * ( cos( _2pi / v) + 1.0));    
01236   }
01237 }
01238 
01239 
01240 // -------------------------------------------------------------------- VV ----
01241 
01242 
01243 template<class M>
01244 void VV<M>::raise(typename M::VertexHandle& _vh, state_t _target_state) 
01245 {
01246   if (MOBJ(_vh).state() < _target_state) 
01247   {
01248     update(_vh, _target_state);
01249 
01250     // raise all neighbour vertices to level x-1
01251     typename M::VertexVertexIter              vv_it(Base::mesh_.vv_iter(_vh));
01252     typename M::VertexHandle                  vh;
01253     std::vector<typename M::VertexHandle>     vertex_vector;
01254 
01255     if (_target_state > 1) {
01256 
01257       for (; vv_it; ++vv_it) {
01258 
01259         vertex_vector.push_back(vv_it.handle());
01260       }
01261 
01262       while (!vertex_vector.empty()) {
01263 
01264         vh = vertex_vector.back();
01265         vertex_vector.pop_back();
01266 
01267         Base::prev_rule()->raise(vh, _target_state - 1);
01268       }
01269 
01270       for (; vv_it; ++vv_it) {
01271 
01272         vertex_vector.push_back(vv_it.handle());
01273       }
01274 
01275       while (!vertex_vector.empty()) {
01276 
01277         vh = vertex_vector.back();
01278         vertex_vector.pop_back();
01279 
01280         Base::prev_rule()->raise(vh, _target_state - 1);
01281       }
01282     }
01283 
01284     // calculate new position
01285     typename M::Point position(0.0, 0.0, 0.0);
01286     int                  valence(0);
01287 
01288     for (vv_it = Base::mesh_.vv_iter(_vh); vv_it; ++vv_it) {
01289 
01290       ++valence;
01291 
01292       position += Base::mesh_.data(vv_it).position(_target_state - 1);
01293     }
01294 
01295     position /= valence;
01296 
01297     MOBJ(_vh).set_position(_target_state, position);
01298     MOBJ(_vh).inc_state();
01299 
01300     // check if last rule
01301     if (Base::number() == Base::n_rules() - 1) {
01302 
01303       Base::mesh_.set_point(_vh, position);
01304       MOBJ(_vh).set_final();
01305     }
01306   }
01307 }
01308 
01309 
01310 // ------------------------------------------------------------------- VVc ----
01311 
01312 
01313 template<class M>
01314 void VVc<M>::raise(typename M::VertexHandle& _vh, state_t _target_state)
01315 {
01316   if (MOBJ(_vh).state() < _target_state) {
01317 
01318     update(_vh, _target_state);
01319 
01320     // raise all neighbour vertices to level x-1
01321     typename M::VertexVertexIter              vv_it(Base::mesh_.vv_iter(_vh));
01322     typename M::VertexHandle                  vh;
01323     std::vector<typename M::VertexHandle>     vertex_vector;
01324 
01325     if (_target_state > 1) {
01326 
01327       for (; vv_it; ++vv_it) {
01328 
01329         vertex_vector.push_back(vv_it.handle());
01330       }
01331 
01332       while (!vertex_vector.empty()) {
01333 
01334         vh = vertex_vector.back();
01335         vertex_vector.pop_back();
01336 
01337         Base::prev_rule()->raise(vh, _target_state - 1);
01338       }
01339 
01340       for (; vv_it; ++vv_it) {
01341 
01342         vertex_vector.push_back(vv_it.handle());
01343       }
01344 
01345       while (!vertex_vector.empty()) {
01346 
01347         vh = vertex_vector.back();
01348         vertex_vector.pop_back();
01349 
01350         Base::prev_rule()->raise(vh, _target_state - 1);
01351       }
01352     }
01353 
01354     // calculate new position
01355     typename M::Point position(0.0, 0.0, 0.0);
01356     int                  valence(0);
01357     typename M::Scalar c;
01358 
01359     for (vv_it = Base::mesh_.vv_iter(_vh); vv_it; ++vv_it) 
01360     {
01361       ++valence;
01362       position += Base::mesh_.data(vv_it).position(_target_state - 1);
01363     }
01364 
01365     position /= valence;
01366 
01367     // choose coefficcient c
01368     c = Base::coeff();
01369 
01370     position *= (1.0 - c);
01371     position += MOBJ(_vh).position(_target_state - 1) * c;
01372 
01373     MOBJ(_vh).set_position(_target_state, position);
01374     MOBJ(_vh).inc_state();
01375 
01376     if (Base::number() == Base::n_rules() - 1) 
01377     {
01378       Base::mesh_.set_point(_vh, position);
01379       MOBJ(_vh).set_final();
01380     }
01381   }
01382 }
01383 
01384 
01385 // -------------------------------------------------------------------- VE ----
01386 
01387 
01388 template<class M>
01389 void VE<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) 
01390 {
01391   if (MOBJ(_eh).state() < _target_state) {
01392 
01393     update(_eh, _target_state);
01394 
01395     // raise all neighbour vertices to level x-1
01396     typename M::VertexHandle          vh;
01397     typename M::HalfedgeHandle        hh1(Base::mesh_.HEH(_eh, 0)),
01398                                          hh2(Base::mesh_.HEH(_eh, 1));
01399 
01400     if (_target_state > 1) {
01401 
01402       vh = Base::mesh_.TVH(hh1);
01403       
01404       Base::prev_rule()->raise(vh, _target_state - 1);
01405 
01406       vh = Base::mesh_.TVH(hh2);
01407       
01408       Base::prev_rule()->raise(vh, _target_state - 1);
01409     }
01410 
01411     // calculate new position
01412     typename M::Point position(0.0, 0.0, 0.0);
01413     int                  valence(0);
01414 
01415     valence = 2;
01416     position += MOBJ(Base::mesh_.TVH(hh1)).position(_target_state - 1);
01417     position += MOBJ(Base::mesh_.TVH(hh2)).position(_target_state - 1);
01418 
01419     position /= valence;
01420 
01421     MOBJ(_eh).set_position(_target_state, position);
01422     MOBJ(_eh).inc_state();
01423   }
01424 }
01425 
01426 
01427 // ------------------------------------------------------------------- VdE ----
01428 
01429 
01430 template<class M>
01431 void VdE<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) 
01432 {
01433   if (MOBJ(_eh).state() < _target_state) 
01434   {
01435     update(_eh, _target_state);
01436 
01437     // raise all neighbour vertices to level x-1
01438     typename M::VertexHandle             vh;
01439     typename M::HalfedgeHandle           hh1(Base::mesh_.HEH(_eh, 0)),
01440                                             hh2(Base::mesh_.HEH(_eh, 1));
01441     std::vector<typename M::VertexHandle> vertex_vector;
01442     typename M::FaceHandle                fh1, fh2;
01443 
01444     if (_target_state > 1) {
01445 
01446       fh1 = Base::mesh_.FH(hh1);
01447       fh2 = Base::mesh_.FH(hh2);
01448 
01449       if (fh1.is_valid()) {
01450 
01451         Base::prev_rule()->raise(fh1, _target_state - 1);
01452 
01453         vh = Base::mesh_.TVH(Base::mesh_.NHEH(hh1));
01454         Base::prev_rule()->raise(vh, _target_state - 1);
01455       }
01456                            
01457       if (fh2.is_valid()) {
01458 
01459         Base::prev_rule()->raise(fh2, _target_state - 1);
01460 
01461         vh = Base::mesh_.TVH(Base::mesh_.NHEH(hh2));
01462         Base::prev_rule()->raise(vh, _target_state - 1);
01463       }
01464                            
01465       vh = Base::mesh_.TVH(hh1);
01466       Base::prev_rule()->raise(vh, _target_state - 1);
01467 
01468       vh = Base::mesh_.TVH(hh2);
01469       Base::prev_rule()->raise(vh, _target_state - 1);
01470     }
01471 
01472     // calculate new position
01473     typename M::Point position(0.0, 0.0, 0.0);
01474     int                  valence(0);
01475 
01476     valence = 2;
01477     position += MOBJ(Base::mesh_.TVH(hh1)).position(_target_state - 1);
01478     position += MOBJ(Base::mesh_.TVH(hh2)).position(_target_state - 1);
01479 
01480     if (fh1.is_valid()) {
01481 
01482       position += MOBJ(Base::mesh_.TVH(Base::mesh_.NHEH(hh1))).position(_target_state - 1);
01483       ++valence;
01484     }
01485 
01486     if (fh2.is_valid()) {
01487       
01488       position += MOBJ(Base::mesh_.TVH(Base::mesh_.NHEH(hh2))).position(_target_state - 1);
01489       ++valence;
01490     }
01491 
01492     if (Base::number() == Base::subdiv_rule()->Base::number() + 1) 
01493       valence = 4;
01494 
01495     position /= valence;
01496 
01497     MOBJ(_eh).set_position(_target_state, position);
01498     MOBJ(_eh).inc_state();
01499   }
01500 }
01501 
01502 
01503 // ------------------------------------------------------------------ VdEc ----
01504 
01505 
01506 template<class M>
01507 void
01508 VdEc<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) 
01509 {
01510   if (MOBJ(_eh).state() < _target_state) 
01511   {
01512     update(_eh, _target_state);
01513 
01514     // raise all neighbour vertices to level x-1
01515     typename M::VertexHandle             vh;
01516     typename M::HalfedgeHandle           hh1(Base::mesh_.HEH(_eh, 0)),
01517                                             hh2(Base::mesh_.HEH(_eh, 1));
01518     std::vector<typename M::VertexHandle> vertex_vector;
01519     typename M::FaceHandle                fh1, fh2;
01520 
01521     if (_target_state > 1) {
01522 
01523       fh1 = Base::mesh_.FH(Base::mesh_.HEH(_eh, 0));
01524       fh2 = Base::mesh_.FH(Base::mesh_.HEH(_eh, 1));
01525 
01526       Base::prev_rule()->raise(fh1, _target_state - 1);
01527       Base::prev_rule()->raise(fh2, _target_state - 1);
01528 
01529       vertex_vector.push_back(Base::mesh_.TVH(hh1));
01530       vertex_vector.push_back(Base::mesh_.TVH(hh2));
01531 
01532       vertex_vector.push_back(Base::mesh_.TVH(Base::mesh_.NHEH(hh1)));
01533       vertex_vector.push_back(Base::mesh_.TVH(Base::mesh_.NHEH(hh2)));
01534 
01535       while (!vertex_vector.empty()) {
01536 
01537         vh = vertex_vector.back();
01538         vertex_vector.pop_back();
01539 
01540         Base::prev_rule()->raise(vh, _target_state - 1);
01541       }
01542 
01543       vertex_vector.push_back(Base::mesh_.TVH(hh1));
01544       vertex_vector.push_back(Base::mesh_.TVH(hh2));
01545 
01546       vertex_vector.push_back(Base::mesh_.TVH(Base::mesh_.NHEH(hh1)));
01547       vertex_vector.push_back(Base::mesh_.TVH(Base::mesh_.NHEH(hh2)));
01548 
01549       while (!vertex_vector.empty()) {
01550 
01551         vh = vertex_vector.back();
01552         vertex_vector.pop_back();
01553 
01554         Base::prev_rule()->raise(vh, _target_state - 1);
01555       }
01556     }
01557 
01558     // calculate new position
01559     typename M::Point position(0.0, 0.0, 0.0);
01560     int                  valence(0);
01561     typename M::Scalar c;
01562 
01563     // choose coefficient c
01564     c = Base::coeff();
01565 
01566     valence = 4;
01567     position += MOBJ(Base::mesh_.TVH(hh1)).position(_target_state - 1) * c;
01568     position += MOBJ(Base::mesh_.TVH(hh2)).position(_target_state - 1) * c;
01569     position += MOBJ(Base::mesh_.TVH(Base::mesh_.NHEH(hh1))).position(_target_state - 1) * (0.5 - c);
01570     position += MOBJ(Base::mesh_.TVH(Base::mesh_.NHEH(hh2))).position(_target_state - 1) * (0.5 - c);
01571 
01572     position /= valence;
01573 
01574     MOBJ(_eh).set_position(_target_state, position);
01575     MOBJ(_eh).inc_state();
01576   }
01577 }
01578 
01579 
01580 // -------------------------------------------------------------------- EV ----
01581 
01582 
01583 template<class M>
01584 void EV<M>::raise(typename M::VertexHandle& _vh, state_t _target_state) 
01585 {
01586   if (MOBJ(_vh).state() < _target_state) {
01587 
01588     update(_vh, _target_state);
01589 
01590     // raise all neighbour vertices to level x-1
01591     typename M::VertexEdgeIter            ve_it(Base::mesh_.ve_iter(_vh));
01592     typename M::EdgeHandle                eh;
01593     std::vector<typename M::EdgeHandle>   edge_vector;
01594 
01595     if (_target_state > 1) {
01596 
01597       for (; ve_it; ++ve_it) {
01598 
01599         edge_vector.push_back(ve_it.handle());
01600       }
01601 
01602       while (!edge_vector.empty()) {
01603 
01604         eh = edge_vector.back();
01605         edge_vector.pop_back();
01606 
01607         Base::prev_rule()->raise(eh, _target_state - 1);
01608       }
01609 
01610       for (ve_it = Base::mesh_.ve_iter(_vh); ve_it; ++ve_it) {
01611 
01612         edge_vector.push_back(ve_it.handle());
01613       }
01614 
01615       while (!edge_vector.empty()) {
01616 
01617         eh = edge_vector.back();
01618         edge_vector.pop_back();
01619 
01620         while (MOBJ(eh).state() < _target_state - 1)
01621           Base::prev_rule()->raise(eh, _target_state - 1);
01622       }
01623     }
01624 
01625     // calculate new position
01626     typename M::Point position(0.0, 0.0, 0.0);
01627     int                  valence(0);
01628 
01629     for (ve_it = Base::mesh_.ve_iter(_vh); ve_it; ++ve_it) {
01630 
01631       if (Base::mesh_.data(ve_it).final()) {
01632 
01633         ++valence;
01634 
01635         position += Base::mesh_.data(ve_it).position(_target_state - 1);
01636       }
01637     }
01638 
01639     position /= valence;
01640 
01641     MOBJ(_vh).set_position(_target_state, position);
01642     MOBJ(_vh).inc_state();
01643 
01644     // check if last rule
01645     if (Base::number() == Base::n_rules() - 1) {
01646 
01647       Base::mesh_.set_point(_vh, position);
01648       MOBJ(_vh).set_final();
01649     }
01650   }
01651 }
01652 
01653 
01654 // ------------------------------------------------------------------- EVc ----
01655 
01656 template<class M>
01657 std::vector<double> EVc<M>::coeffs_;
01658 
01659 template<class M>
01660 void EVc<M>::raise(typename M::VertexHandle& _vh, state_t _target_state) 
01661 {
01662   if (MOBJ(_vh).state() < _target_state) 
01663   {
01664     update(_vh, _target_state);
01665 
01666     // raise all neighbour vertices to level x-1
01667     typename M::VertexOHalfedgeIter       voh_it;
01668     typename M::EdgeHandle                eh;
01669     typename M::FaceHandle                fh;
01670     std::vector<typename M::EdgeHandle>   edge_vector;
01671     std::vector<typename M::FaceHandle>   face_vector;
01672 
01673     if (_target_state > 1) {
01674 
01675       for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it) {
01676 
01677         face_vector.push_back(Base::mesh_.FH(voh_it.handle()));
01678       }
01679 
01680       while (!face_vector.empty()) {
01681 
01682         fh = face_vector.back();
01683         face_vector.pop_back();
01684 
01685         if (fh.is_valid())
01686           Base::prev_rule()->raise(fh, _target_state - 1);
01687       }
01688 
01689       for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it) {
01690 
01691         edge_vector.push_back(Base::mesh_.EH(voh_it.handle()));
01692 
01693         edge_vector.push_back(Base::mesh_.EH(Base::mesh_.NHEH(voh_it.handle())));
01694       }
01695 
01696       while (!edge_vector.empty()) {
01697 
01698         eh = edge_vector.back();
01699         edge_vector.pop_back();
01700 
01701         while (MOBJ(eh).state() < _target_state - 1)
01702           Base::prev_rule()->raise(eh, _target_state - 1);
01703       }
01704     }
01705 
01706 
01707     // calculate new position
01708     typename M::Point        position(0.0, 0.0, 0.0);
01709     typename M::Scalar       c;
01710     typename M::Point        zero_point(0.0, 0.0, 0.0);
01711     int                         valence(0);
01712 
01713     valence = Base::mesh_.valence(_vh);
01714     c       = coeff( valence );
01715 
01716     for (voh_it = Base::mesh_.voh_iter(_vh); voh_it; ++voh_it) 
01717     {
01718       if (MOBJ(Base::mesh_.EH(voh_it.handle())).final()) 
01719       {
01720         position += MOBJ(Base::mesh_.EH(voh_it.handle())).position(_target_state-1)*c;
01721 
01722         if ( Base::mesh_.FH(voh_it.handle()).is_valid() && 
01723              MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(voh_it.handle()))).final() && 
01724              MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(voh_it.handle()))).position(_target_state - 1) != zero_point) 
01725         {
01726           position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(voh_it.handle()))).position(_target_state-1) * (1.0-c);
01727         }
01728         else {
01729           position += MOBJ(Base::mesh_.EH(voh_it.handle())).position(_target_state - 1) * (1.0 - c);
01730         }
01731       }
01732       else {
01733         --valence;
01734       }
01735     }
01736 
01737     position /= valence;
01738 
01739     MOBJ(_vh).set_position(_target_state, position);
01740     MOBJ(_vh).inc_state();
01741 
01742     // check if last rule
01743     if (Base::number() == Base::n_rules() - 1) {
01744 
01745       Base::mesh_.set_point(_vh, position);
01746       MOBJ(_vh).set_final();
01747     }
01748   }
01749 }
01750 
01751 
01752 template <class M>
01753 void 
01754 EVc<M>::init_coeffs(size_t _max_valence)
01755 {
01756   if ( coeffs_.size() == _max_valence+1 ) // equal? do nothing
01757     return;
01758 
01759   if (coeffs_.size() < _max_valence+1) // less than? add additional valences
01760   {
01761     const double _2pi = 2.0*M_PI;
01762     double c;
01763   
01764     if (coeffs_.empty())
01765       coeffs_.push_back(0.0); // dummy for invalid valences 0,1,2
01766 
01767     for(size_t v=coeffs_.size(); v <= _max_valence; ++v)
01768     {
01769       // ( 3/2 + cos ( 2 PI / valence ) )� / 2 - 1
01770       c = 1.5 + cos( _2pi / v );
01771       c = c * c * 0.5 - 1.0;
01772       coeffs_.push_back(c);
01773     }
01774   }
01775 }
01776 
01777 
01778 // -------------------------------------------------------------------- EF ----
01779 
01780 template<class M>
01781 void
01782 EF<M>::raise(typename M::FaceHandle& _fh, state_t _target_state) {
01783 
01784   if (MOBJ(_fh).state() < _target_state) {
01785 
01786     update(_fh, _target_state);
01787 
01788     // raise all neighbour edges to level x-1
01789     typename M::FaceEdgeIter              fe_it(Base::mesh_.fe_iter(_fh));
01790     typename M::EdgeHandle                eh;
01791     std::vector<typename M::EdgeHandle>   edge_vector;
01792 
01793     if (_target_state > 1) {
01794 
01795       for (; fe_it; ++fe_it) {
01796 
01797         edge_vector.push_back(fe_it.handle());
01798       }
01799 
01800       while (!edge_vector.empty()) {
01801 
01802         eh = edge_vector.back();
01803         edge_vector.pop_back();
01804 
01805         Base::prev_rule()->raise(eh, _target_state - 1);
01806       }
01807 
01808       for (fe_it = Base::mesh_.fe_iter(_fh); fe_it; ++fe_it) {
01809 
01810         edge_vector.push_back(fe_it.handle());
01811       }
01812 
01813       while (!edge_vector.empty()) {
01814 
01815         eh = edge_vector.back();
01816         edge_vector.pop_back();
01817 
01818         while (MOBJ(eh).state() < _target_state - 1)
01819           Base::prev_rule()->raise(eh, _target_state - 1);
01820       }
01821     }
01822 
01823     // calculate new position
01824     typename M::Point position(0.0, 0.0, 0.0);
01825     int                  valence(0);
01826 
01827     for (fe_it = Base::mesh_.fe_iter(_fh); fe_it; ++fe_it) {
01828 
01829       if (Base::mesh_.data(fe_it).final()) {
01830 
01831         ++valence;
01832 
01833         position += Base::mesh_.data(fe_it).position(_target_state - 1);
01834       }
01835     }
01836 
01837     assert (valence == 3);
01838 
01839     position /= valence;
01840 
01841     MOBJ(_fh).set_position(_target_state, position);
01842     MOBJ(_fh).inc_state();
01843   }
01844 }
01845 
01846 
01847 // -------------------------------------------------------------------- FE ----
01848 
01849 
01850 template<class M>
01851 void
01852 FE<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) {
01853 
01854   if (MOBJ(_eh).state() < _target_state) {
01855 
01856     update(_eh, _target_state);
01857 
01858     // raise all neighbour faces to level x-1
01859     typename M::FaceHandle                fh;
01860 
01861     if (_target_state > 1) {
01862 
01863       fh = Base::mesh_.FH(Base::mesh_.HEH(_eh, 0));
01864       Base::prev_rule()->raise(fh, _target_state - 1);
01865 
01866       fh = Base::mesh_.FH(Base::mesh_.HEH(_eh, 1));
01867       Base::prev_rule()->raise(fh, _target_state - 1);
01868 
01869       fh = Base::mesh_.FH(Base::mesh_.HEH(_eh, 0));
01870       Base::prev_rule()->raise(fh, _target_state - 1);
01871 
01872       fh = Base::mesh_.FH(Base::mesh_.HEH(_eh, 1));
01873       Base::prev_rule()->raise(fh, _target_state - 1);
01874     }
01875 
01876     // calculate new position
01877     typename M::Point position(0.0, 0.0, 0.0);
01878     int                  valence(2);
01879 
01880     position += MOBJ(Base::mesh_.FH(Base::mesh_.HEH(_eh, 0))).position(_target_state - 1);
01881     
01882     position += MOBJ(Base::mesh_.FH(Base::mesh_.HEH(_eh, 1))).position(_target_state - 1);
01883     
01884     position /= valence;
01885 
01886     MOBJ(_eh).set_position(_target_state, position);
01887     MOBJ(_eh).inc_state();
01888   }
01889 }
01890 
01891 
01892 // ------------------------------------------------------------------- EdE ----
01893 
01894 
01895 template<class M>
01896 void
01897 EdE<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) {
01898 
01899   if (MOBJ(_eh).state() < _target_state) {
01900 
01901     update(_eh, _target_state);
01902 
01903     // raise all neighbour faces and edges to level x-1
01904     typename M::HalfedgeHandle            hh1, hh2;
01905     typename M::FaceHandle                fh;
01906     typename M::EdgeHandle                eh;
01907 
01908     hh1 = Base::mesh_.HEH(_eh, 0);
01909     hh2 = Base::mesh_.HEH(_eh, 1);
01910 
01911     if (_target_state > 1) {
01912 
01913       fh = Base::mesh_.FH(hh1);
01914       Base::prev_rule()->raise(fh, _target_state - 1);
01915 
01916       fh = Base::mesh_.FH(hh2);
01917       Base::prev_rule()->raise(fh, _target_state - 1);
01918 
01919       eh = Base::mesh_.EH(Base::mesh_.NHEH(hh1));
01920       Base::prev_rule()->raise(eh, _target_state - 1);
01921 
01922       eh = Base::mesh_.EH(Base::mesh_.PHEH(hh1));
01923       Base::prev_rule()->raise(eh, _target_state - 1);
01924 
01925       eh = Base::mesh_.EH(Base::mesh_.NHEH(hh2));
01926       Base::prev_rule()->raise(eh, _target_state - 1);
01927 
01928       eh = Base::mesh_.EH(Base::mesh_.PHEH(hh2));
01929       Base::prev_rule()->raise(eh, _target_state - 1);
01930     }
01931 
01932     // calculate new position
01933     typename M::Point position(0.0, 0.0, 0.0);
01934     int                  valence(4);
01935 
01936     position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh1))).position(_target_state - 1);
01937     position += MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh1))).position(_target_state - 1);
01938     position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh2))).position(_target_state - 1);
01939     position += MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh2))).position(_target_state - 1);
01940 
01941     position /= valence;
01942 
01943     MOBJ(_eh).set_position(_target_state, position);
01944     MOBJ(_eh).inc_state();
01945   }
01946 }
01947 
01948 
01949 // ------------------------------------------------------------------ EdEc ----
01950 
01951 
01952 template<class M>
01953 void
01954 EdEc<M>::raise(typename M::EdgeHandle& _eh, state_t _target_state) 
01955 {
01956   if (MOBJ(_eh).state() < _target_state) {
01957 
01958     update(_eh, _target_state);
01959 
01960     // raise all neighbour faces and edges to level x-1
01961     typename M::HalfedgeHandle            hh1, hh2;
01962     typename M::FaceHandle                fh;
01963     typename M::EdgeHandle                eh;
01964 
01965     hh1 = Base::mesh_.HEH(_eh, 0);
01966     hh2 = Base::mesh_.HEH(_eh, 1);
01967 
01968     if (_target_state > 1) {
01969 
01970       fh = Base::mesh_.FH(hh1);
01971       Base::prev_rule()->raise(fh, _target_state - 1);
01972 
01973       fh = Base::mesh_.FH(hh2);
01974       Base::prev_rule()->raise(fh, _target_state - 1);
01975 
01976       eh = Base::mesh_.EH(Base::mesh_.NHEH(hh1));
01977       Base::prev_rule()->raise(eh, _target_state - 1);
01978 
01979       eh = Base::mesh_.EH(Base::mesh_.PHEH(hh1));
01980       Base::prev_rule()->raise(eh, _target_state - 1);
01981 
01982       eh = Base::mesh_.EH(Base::mesh_.NHEH(hh2));
01983       Base::prev_rule()->raise(eh, _target_state - 1);
01984 
01985       eh = Base::mesh_.EH(Base::mesh_.PHEH(hh2));
01986       Base::prev_rule()->raise(eh, _target_state - 1);
01987     }
01988 
01989     // calculate new position
01990     typename M::Point  position(0.0, 0.0, 0.0);
01991     int                valence(4);
01992     typename M::Scalar c;
01993 
01994     position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh1))).position(_target_state - 1);
01995     position += MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh1))).position(_target_state - 1);
01996     position += MOBJ(Base::mesh_.EH(Base::mesh_.NHEH(hh2))).position(_target_state - 1);
01997     position += MOBJ(Base::mesh_.EH(Base::mesh_.PHEH(hh2))).position(_target_state - 1);
01998 
01999     position /= valence;
02000 
02001     // choose coefficient c
02002     c = Base::coeff();
02003 
02004     position *= (1.0 - c);
02005 
02006     position += MOBJ(_eh).position(_target_state - 1) * c;
02007 
02008     MOBJ(_eh).set_position(_target_state, position);
02009     MOBJ(_eh).inc_state();
02010   }
02011 }
02012 
02013 #endif // DOXY_IGNORE_THIS
02014 
02015 #undef FH 
02016 #undef VH 
02017 #undef EH 
02018 #undef HEH
02019 #undef M
02020 #undef MOBJ
02021 
02022 //=============================================================================
02023 } // END_NS_ADAPTIVE
02024 } // END_NS_SUBDIVIDER
02025 } // END_NS_OPENMESH
02026 //=============================================================================