00001 #ifndef OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_COMPUTE_VERTEX_MEAN_CURVATURE_NORMAL_H 00002 #define OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_COMPUTE_VERTEX_MEAN_CURVATURE_NORMAL_H 00003 // 00004 // OpenTissue Template Library 00005 // - A generic toolbox for physics-based modeling and simulation. 00006 // Copyright (C) 2008 Department of Computer Science, University of Copenhagen. 00007 // 00008 // OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php 00009 // 00010 #include <OpenTissue/configuration.h> 00011 00012 #include <OpenTissue/core/geometry/geometry_cot.h> 00013 #include <OpenTissue/core/geometry/geometry_compute_area_mixed.h> 00014 #include <OpenTissue/core/containers/mesh/polymesh/util/polymesh_is_boundary.h> 00015 #include <OpenTissue/core/containers/mesh/polymesh/util/polymesh_is_vertex_neighbors_triangular.h> 00016 #include <OpenTissue/core/math/math_precision.h> // Needed for math::machine_precision 00017 00018 #include <cassert> // Needed for assert 00019 00020 namespace OpenTissue 00021 { 00022 namespace polymesh 00023 { 00024 00037 template<typename mesh_type> 00038 bool compute_vertex_mean_curvature_normal(PolyMeshVertex<mesh_type> const & v, typename mesh_type::math_types::vector3_type & Kh) 00039 { 00040 typedef PolyMeshVertex<mesh_type> vertex_type; 00041 typedef typename mesh_type::vertex_iterator vertex_iterator; 00042 typedef typename mesh_type::face_iterator face_iterator; 00043 typedef typename mesh_type::vertex_face_circulator vertex_face_circulator; 00044 typedef typename mesh_type::vertex_halfedge_circulator vertex_halfedge_circulator; 00045 00046 typedef typename mesh_type::math_types math_types; 00047 typedef typename math_types::vector3_type vector3_type; 00048 typedef typename math_types::real_type real_type; 00049 00050 static real_type const epsilon = math::machine_precision<real_type>(); 00051 static real_type const zero = real_type(); //--- by standard default constructed integral types is zero 00052 static real_type const two = boost::numeric_cast<real_type>(2.0); 00053 00054 if(is_boundary(v)) 00055 { 00056 std::cerr << "compute_vertex_mean_curvature_normal(): boundary vertex encountered" << std::endl; 00057 return false; 00058 } 00059 if(!is_vertex_neighbors_triangular(v)) 00060 { 00061 std::cout << "compute_vertex_mean_curvature_normal(): non-triangular face encountered" << std::endl; 00062 return false; 00063 } 00064 00065 real_type area = zero; 00066 Kh.clear(); // By OT convention should set all elements to zero. 00067 vector3_type p = v.m_coord; 00068 00069 //using std::max; 00070 //using std::min; 00071 //{ 00072 // real_type one = 1.0; 00073 // vertex_halfedge_circulator h(v),hend; 00074 // for(;h!=hend;++h) 00075 // { 00076 // vector3_type nbr = h->get_destination_iterator()->m_coord; 00077 // vector3_type left = h->get_twin_iterator()->get_prev_iterator()->get_origin_iterator()->m_coord; 00078 // vector3_type right = h->get_next_iterator()->get_destination_iterator()->m_coord; 00079 // real_type d_left = unit(nbr-left)*unit(p-left); 00080 // real_type d_right = unit(nbr-right)*unit(p-right); 00081 // real_type a_left = acos( min(one, max(-one, d_left )) ); 00082 // real_type a_right = acos( min(one, max(-one, d_right)) ); 00083 // real_type w = 1.0/tan(a_left) + 1.0/tan(a_right); 00084 // Kh += w * (p-nbr); 00085 // } 00086 //} 00087 //{ 00088 // vertex_halfedge_circulator h(v),hend; 00089 // for(;h!=hend;++h) 00090 // { 00091 // vector3_type pi = h->get_next_iterator()->get_origin_iterator()->m_coord; 00092 // vector3_type pj = h->get_next_iterator()->get_destination_iterator()->m_coord; 00093 // area += compute_area_mixed(p,pi,pj); 00094 // } 00095 //} 00096 //--- 00097 //--- Orignally we have 00098 //--- 00099 //--- p // 00100 //--- /|\ // 00101 //--- / | \ // 00102 //--- alpha_k / | \ beta_k K_p = 1/A_p * sum_k (cot(alpha_k)+cot(beta_k))*(p-p_k) // 00103 //--- \ | / // 00104 //--- \ | / // 00105 //--- \|/ // 00106 //--- p_k // 00107 //--- 00108 //--- But we could just as easily re-write this as 00109 //--- 00110 //--- p // 00111 //--- /|\ // 00112 //--- / | \ // 00113 //--- alpha_j / | \ beta_i K_p = 1/A_p * sum_i ( cot(alpha_j)(p-p_j) + cot(beta_i))*(p-p_i) ) // 00114 //--- /---+---\ // 00115 //--- p_i p_j // 00116 //--- 00117 //--- which is the forumula we have implemented below 00118 vertex_halfedge_circulator h(v),hend; 00119 for(;h!=hend;++h) 00120 { 00121 assert( h->get_next_iterator()->get_origin_iterator()->get_handle().get_idx() != v.get_handle().get_idx() || !"compute_vertex_mean_curvature_normal(): illegal topology"); 00122 assert( h->get_next_iterator()->get_destination_iterator()->get_handle().get_idx() != v.get_handle().get_idx() || !"compute_vertex_mean_curvature_normal(): illegal topology"); 00123 vector3_type pi = h->get_next_iterator()->get_origin_iterator()->m_coord; 00124 vector3_type pj = h->get_next_iterator()->get_destination_iterator()->m_coord; 00125 area += OpenTissue::geometry::compute_area_mixed(p,pi,pj); 00126 real_type cotan_alpha = OpenTissue::geometry::cot( pi, pj, p ); 00127 real_type cotan_beta = OpenTissue::geometry::cot( pj, p, pi ); 00128 Kh += (cotan_alpha*(p-pj)) + (cotan_beta*(p-pi)); 00129 } 00130 if (area > epsilon) 00131 Kh /= two*area; //--- hmmm, I thinnk there is a bug here, it should read: Kh /= four*area;??? 00132 else 00133 { 00134 //std::cout << "compute_vertex_mean_curvature_normal(): zero area clamping" << std::endl; 00135 Kh.clear(); 00136 return false; 00137 } 00138 if((Kh*Kh)<epsilon) 00139 { 00140 //std::cout << "compute_vertex_mean_curvature_normal(): too small curvature normal clamping" << std::endl; 00141 Kh.clear(); 00142 } 00143 return true; 00144 } 00145 00146 } // namespace polymesh 00147 } // namespace OpenTissue 00148 00149 //OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_COMPUTE_VERTEX_MEAN_CURVATURE_NORMAL_H 00150 #endif
1.7.1