Go to the documentation of this file.00001 #ifndef OPENTISSUE_COLLISION_SDF_SDF_COMPUTE_POINT_SAMPLING_H
00002 #define OPENTISSUE_COLLISION_SDF_SDF_COMPUTE_POINT_SAMPLING_H
00003
00004
00005
00006
00007
00008
00009
00010 #include <OpenTissue/configuration.h>
00011
00012 #include <OpenTissue/core/containers/mesh/mesh.h>
00013
00014 #include <boost/cast.hpp>
00015
00016 #include <list>
00017
00018 namespace OpenTissue
00019 {
00020 namespace sdf
00021 {
00022
00040 template<typename mesh_type,typename grid_type, typename point_container>
00041 void compute_point_sampling(
00042 mesh_type & mesh
00043 , grid_type const & phi
00044 , double edge_resolution
00045 , bool face_sampling
00046 , point_container & points
00047 )
00048 {
00049 using std::min;
00050 using std::max;
00051 using std::sqrt;
00052
00053 typedef typename mesh_type::vertex_iterator vertex_iterator;
00054 typedef typename mesh_type::halfedge_iterator halfedge_iterator;
00055 typedef typename mesh_type::face_iterator face_iterator;
00056 typedef typename mesh_type::face_type face_type;
00057 typedef typename mesh_type::halfedge_type halfedge_type;
00058 typedef typename mesh_type::face_halfedge_circulator face_halfedge_circulator;
00059 typedef typename std::list<face_type*> face_queue;
00060
00061 typedef typename mesh_type::math_types math_types;
00062 typedef typename math_types::vector3_type vector3_type;
00063 typedef typename math_types::real_type real_type;
00064
00065 assert(edge_resolution>=0 || !"compute_point_sampling(): edge resolution was negative");
00066
00067 mesh::clear_vertex_tags( mesh);
00068 mesh::clear_halfedge_tags( mesh);
00069 mesh::clear_face_tags( mesh);
00070
00071 points.clear();
00072
00073
00074 for(vertex_iterator v = mesh.vertex_begin();v!=mesh.vertex_end();++v)
00075 {
00076 v->m_tag = 1;
00077 if(!is_convex( *v ) )
00078 continue;
00079 points.push_back( v->m_coord );
00080 }
00081
00082
00083
00084 real_type tmp = boost::numeric_cast<real_type>( edge_resolution );
00085 real_type threshold = max(tmp, sqrt( phi.dx()*phi.dx() + phi.dy()*phi.dy() + phi.dz()*phi.dz() ));
00086
00087 for(halfedge_iterator h = mesh.halfedge_begin();h!=mesh.halfedge_end();++h)
00088 {
00089 if(h->m_tag)
00090 continue;
00091 h->m_tag = 1;
00092 h->get_twin_iterator()->m_tag = 1;
00093 if(!is_convex( *h ) )
00094 continue;
00095 vector3_type u = h->get_destination_iterator()->m_coord - h->get_origin_iterator()->m_coord;
00096 real_type lgth = sqrt(u*u);
00097 if(lgth>threshold)
00098 {
00099 u /= lgth;
00100 vector3_type p = h->get_origin_iterator()->m_coord;
00101 real_type t = threshold;
00102 while(t<lgth)
00103 {
00104 p += u*threshold;
00105 t += threshold;
00106 points.push_back( p );
00107 }
00108 }
00109 }
00110
00111
00112
00113 if(face_sampling)
00114 {
00115 vector3_type Ai,ai,ei;
00116 real_type area_test = max( phi.dx()*phi.dy(), max(phi.dx()*phi.dz(),phi.dy()*phi.dz()));
00117
00118 for(face_iterator face = mesh.face_begin();face!=mesh.face_end();++face)
00119 {
00120 if(face->m_tag)
00121 continue;
00122 real_type area = 0;
00123 vector3_type centroid = vector3_type(0,0,0);
00124 unsigned int size = 0;
00125 face_queue Q;
00126 Q.push_back( &(*face) );
00127 face->m_tag = 1;
00128 while(!Q.empty())
00129 {
00130 face_type * cur = Q.front();Q.pop_front();
00131 face_halfedge_circulator h(*cur),hend;
00132 for(;h!=hend;++h)
00133 {
00134 ai = h->get_origin_iterator()->m_coord;
00135 ei = h->get_destination_iterator()->m_coord - ai;
00136 Ai = ai % ei;
00137 area += 0.5*sqrt(Ai*Ai);
00138 ++size;
00139 centroid += h->get_origin_iterator()->m_coord;
00140
00141 if(h->get_twin_iterator()->get_face_handle().is_null())
00142 continue;
00143
00144 face_type * neighbor = &(*h->get_twin_iterator()->get_face_iterator());
00145 bool unseen = !neighbor->m_tag;
00146
00147 bool coplanar = is_planar(*h);
00148 if(unseen && coplanar)
00149 {
00150 neighbor->m_tag = 1;
00151 Q.push_back(neighbor);
00152 }
00153 }
00154 }
00155 if(size && area > area_test)
00156 {
00157 centroid /= size;
00158 points.push_back( centroid );
00159 }
00160 }
00161 }
00162 }
00163
00164 }
00165
00166 }
00167
00168
00169 #endif