00001 #ifndef OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_DIHEDRAL_ANGLE_H 00002 #define OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_DIHEDRAL_ANGLE_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/containers/mesh/polymesh/polymesh_halfedge.h> 00013 00014 namespace OpenTissue 00015 { 00016 namespace polymesh 00017 { 00018 00021 template<typename mesh_type,typename real_type> 00022 void compute_dihedral_angle(PolyMeshHalfEdge<mesh_type> const & e, real_type & radian) 00023 { 00024 typedef typename mesh_type::math_types math_types; 00025 typedef typename math_types::vector3_type vector3_type; 00026 typedef typename mesh_type::vertex_iterator vertex_iterator; 00027 typedef typename mesh_type::face_iterator face_iterator; 00028 static real_type threshold = math::working_precision<real_type>(); 00029 00030 if( e.get_face_handle().is_null() || e.get_twin_iterator()->get_face_handle().is_null() ) 00031 { 00032 std::cout << "compute_dihedral_angle(): boundary edge encountered!" << std::endl; 00033 radian = 0; 00034 return; 00035 } 00036 00037 face_iterator face = e.get_face_iterator(); 00038 face_iterator twin_face = e.get_twin_iterator()->get_face_iterator(); 00039 00040 vector3_type n1,n2; 00041 compute_face_normal(*face, n1); 00042 compute_face_normal(*twin_face, n2); 00043 00044 real_type denom = length(n1)*length(n2); 00045 if ( denom < threshold) 00046 { 00047 radian = 0; 00048 return; 00049 } 00050 real_type cos_angle = n1*n2 / denom; 00051 cos_angle = (cos_angle>1.0)?1.0: ( (cos_angle<-1.0)?-1.0:cos_angle ); 00052 00053 vertex_iterator d = e.get_destination_iterator(); 00054 vertex_iterator o = e.get_origin_iterator(); 00055 vector3_type u = d->m_coord - o->m_coord; 00056 real_type sin_angle = (n1 % n2) * u; 00057 sin_angle = (sin_angle>1.0)?1.0: ( (sin_angle<-1.0)?-1.0:sin_angle ); 00058 radian = sin_angle >= 0 ? acos(cos_angle) : -acos(cos_angle); 00059 } 00060 00061 template<typename mesh_type,typename real_type> 00062 void compute_dihedral_angle(PolyMeshEdge<mesh_type> const & e, real_type & radian) 00063 { 00064 compute_dihedral_angle( *(e.get_halfedge0_iterator()), radian ); 00065 } 00066 00067 00068 } // namespace polymesh 00069 } // namespace OpenTissue 00070 00071 //OPENTISSUE_CORE_CONTAINERS_MESH_POLYMESH_UTIL_POLYMESH_DIHEDRAL_ANGLE_H 00072 #endif 00073
1.7.1