00001 #ifndef OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_COMPUTE_VERTEX_GAUSSIAN_CURVATURE_H 00002 #define OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_COMPUTE_VERTEX_GAUSSIAN_CURVATURE_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_is_angle_obtuse.h> 00013 #include <OpenTissue/core/geometry/geometry_is_triangle_obtuse.h> 00014 #include <OpenTissue/core/geometry/geometry_angle_from_cot.h> 00015 #include <OpenTissue/core/geometry/geometry_compute_area_mixed.h> 00016 00017 #include <OpenTissue/core/containers/mesh/polymesh/util/polymesh_is_boundary.h> 00018 #include <OpenTissue/core/containers/mesh/polymesh/util/polymesh_is_vertex_neighbors_triangular.h> 00019 00020 #include <OpenTissue/core/math/math_precision.h> // Needed for math::machine_precision 00021 #include <OpenTissue/core/math/math_constants.h> 00022 00023 #include <cassert> // Needed for assert 00024 00025 namespace OpenTissue 00026 { 00027 namespace polymesh 00028 { 00029 00042 template<typename mesh_type, typename real_type2> 00043 bool compute_vertex_gaussian_curvature( PolyMeshVertex<mesh_type> const & v, real_type2 & Kg) 00044 { 00045 typedef PolyMeshVertex<mesh_type> vertex_type; 00046 typedef typename mesh_type::vertex_iterator vertex_iterator; 00047 typedef typename mesh_type::face_iterator face_iterator; 00048 typedef typename mesh_type::vertex_face_circulator vertex_face_circulator; 00049 typedef typename mesh_type::vertex_halfedge_circulator vertex_halfedge_circulator; 00050 00051 typedef typename mesh_type::math_types math_types; 00052 typedef typename math_types::vector3_type vector3_type; 00053 typedef typename math_types::real_type real_type; 00054 00055 static real_type const epsilon = math::machine_precision<real_type>(); 00056 static real_type const zero = real_type(); //--- by standard default constructed integral types is zero 00057 static real_type const two = boost::numeric_cast<real_type>(2.0); 00058 00059 if(is_boundary(v)) 00060 { 00061 std::cerr << "compute_vertex_gaussian_curvature(): boundary vertex encountered" << std::endl; 00062 return false; 00063 } 00064 if(!is_vertex_neighbors_triangular(v)) 00065 { 00066 std::cout << "compute_vertex_gaussian_curvature(): non-triangular face encountered" << std::endl; 00067 return false; 00068 } 00069 real_type area = zero; 00070 real_type angle_sum = zero; 00071 vector3_type p = v.m_coord; 00072 //--- 00073 //--- Orignally we have 00074 //--- 00075 //--- p // 00076 //--- /|\ // 00077 //--- / | \ // 00078 //--- alpha_k / | \ beta_k K_p = 1/A_p * sum_k (cot(alpha_k)+cot(beta_k))*(p-p_k) // 00079 //--- \ | / // 00080 //--- \ | / // 00081 //--- \|/ // 00082 //--- p_k // 00083 //--- 00084 //--- But we could just as easily re-write this as 00085 //--- 00086 //--- p // 00087 //--- /|\ // 00088 //--- / | \ // 00089 //--- alpha_j / | \ beta_i K_p = 1/A_p * sum_i ( cot(alpha_j)(p-p_j) + cot(beta_i))*(p-p_i) ) // 00090 //--- /-------\ // 00091 //--- p_i p_j // 00092 //--- 00093 //--- which is the forumula we have implemented below 00094 vertex_halfedge_circulator h(v),hend; 00095 for(;h!=hend;++h) 00096 { 00097 assert( h->get_next_iterator()->get_origin_iterator()->get_handle().get_idx() != v.get_handle().get_idx() || !"compute_vertex_mean_curvature_normal(): illegal topology"); 00098 assert( h->get_next_iterator()->get_destination_iterator()->get_handle().get_idx() != v.get_handle().get_idx() || !"compute_vertex_mean_curvature_normal(): illegal topology"); 00099 vector3_type pi = h->get_next_iterator()->get_origin_iterator()->m_coord; 00100 vector3_type pj = h->get_next_iterator()->get_destination_iterator()->m_coord; 00101 area += OpenTissue::geometry::compute_area_mixed(p,pi,pj); 00102 //--- 00103 //--- we could have used the cosine formula 00104 //--- 00105 //--- (pi-p) \cdot (pj-p) 00106 //--- \cos(\theta) = ---------------------- 00107 //--- ||pi-p|| || pj-p || 00108 //--- 00109 //--- However, this involves computing the length of the vectors (pi-p) and (pj-p) 00110 //--- 00111 //--- The angle_from_coton avoids this (no square root computations), at the cost of a atan2 call. 00112 //--- 00113 angle_sum += OpenTissue::geometry::angle_from_cot(p,pi,pj); 00114 } 00115 if (area > epsilon) 00116 Kg = boost::numeric_cast<real_type2> ( (two*math::detail::pi<real_type>() - angle_sum)/area ); 00117 else 00118 Kg = boost::numeric_cast<real_type2>( zero ); 00119 return true; 00120 } 00121 00122 } // namespace polymesh 00123 } // namespace OpenTissue 00124 00125 //OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_COMPUTE_VERTEX_GAUSSIAN_CURVATURE_H 00126 #endif
1.7.1