|
shape_doc 0.1
|
|
shape_doc 0.1
|
00001 #ifndef FIND_MINIMAL_EDGES_HPP 00002 #define FIND_MINIMAL_EDGES_HPP 00003 00004 /***************************************************************************** 00005 * M a t h e m a g i x 00006 ***************************************************************************** 00007 * Find Minimal Edges 00008 * 2009-04-29 00009 * Johan Seland & Bernard Mourrain 00010 ***************************************************************************** 00011 * Copyright (C) 2009 INRIA Sophia-Antipolis 00012 ***************************************************************************** 00013 * Comments : 00014 * find_minimal_edges implement a recursive algorithm to find the 00015 * minimal edges of an kd-tree. Minimal edges are edges of leaf cubes that do not 00016 * properly contain an edge of a neighboring leaf. 00017 * 00018 * 00019 ****************************************************************************/ 00020 00021 #include <tr1/array> 00022 #include <vector> 00023 #include <algorithm> 00024 #include <cassert> 00025 #include <iostream> 00026 00027 #include "shape/node.hpp" 00028 #include "shape/bounding_box.hpp" 00029 00030 //#define ELB_DEBUG 00031 //#define ELB_RUNTIME_CHECK 00032 00033 namespace mmx { namespace shape { 00034 template<class node_t> class EdgeListBuilder { 00035 public: 00036 typedef std::tr1::array<const node_t*, 2> faceTuple_t; 00037 typedef std::vector<faceTuple_t> faceList_t; 00038 00039 typedef std::tr1::array<const node_t*, 3> edgeTuple_t; 00040 typedef std::vector<edgeTuple_t> edgeList_t; 00041 00042 typedef mmx::shape::bounding_box<double> boundingBox_t; 00043 00044 EdgeListBuilder() : m_os( std::cerr ) {} 00045 00046 void cellProc( const node_t* n ); 00047 bool faceProc( const node_t* n0, const node_t* n1, int split_dir ); 00048 void edgeProc( const node_t* n0, const node_t* n1, const node_t* n2, int split_dir ); 00049 00051 int dist( const node_t* n0, const node_t* n1 ) const { 00052 return std::abs( n0->depth() - n1->depth() ); 00053 } 00054 00056 const node_t* getChild( const node_t* n, typename node_t::NODE_TYPE nt ) const { 00057 00058 assert( !n->is_leaf() && "Trying to take children of leaf node." ); 00059 00060 switch ( nt ) { 00061 case node_t::LEFT : return n->left(); break; 00062 case node_t::RIGHT: return n->right(); break; 00063 default: return 0; 00064 } 00065 } 00066 00068 bool commonFace( const faceTuple_t& f ) const; 00069 00070 boundingBox_t computeCommonFace( const faceTuple_t& t ) const; 00071 00074 bool overlapAlongAxis( const node_t* n0, const node_t* n1, int axis, bool strict ) const; 00075 00077 void verifyFaceList() const; 00078 00080 typename node_t::NODE_TYPE getSplitDirParentType( const node_t* n, const int split_dir ) const { 00081 while ( n->parent() != 0 && n->parent()->split_dir() != split_dir ) { 00082 n = n->parent(); 00083 } 00084 if ( n->parent() == 0 ) 00085 std::cerr << "Could not find a parent that was split in " 00086 << split_dir << " direction.\n"; 00087 00088 return n->type(); 00089 } 00090 00092 const edgeList_t& edgeList() const { return m_edgelist; } 00093 00095 edgeList_t& edgeList() { return m_edgelist; } 00096 00098 const faceList_t& faceList() const { return m_facelist; } 00099 00101 faceList_t& faceList() { return m_facelist; } 00102 00103 protected: 00104 void edgeProc_helper( const node_t* n0, const node_t* n1, const node_t* n2, int split_dir ); 00105 edgeList_t m_edgelist; 00106 faceList_t m_facelist; 00107 std::vector<int> m_facesplitdir; 00108 00109 std::ostream& m_os; 00110 }; 00111 00112 template<class node_t> void EdgeListBuilder<node_t>::verifyFaceList() const { 00113 size_t errs = 0; 00114 for ( size_t i = 0; i < m_facelist.size(); ++i ) { 00115 00116 if ( !commonFace( m_facelist[i] ) ) { 00117 errs++; 00118 m_os << "Invalid facetuple found: \n"; 00119 m_os << m_facelist[i][0]->get_cell()->boundingBox() << std::endl; 00120 m_os << m_facelist[i][1]->get_cell()->boundingBox() << std::endl; 00121 m_os << m_facelist[i][0]->split_dir() << ", " 00122 << m_facelist[i][1]->split_dir() << "\n"; 00123 m_os << "Split dir: " << m_facesplitdir[i] << std::endl; 00124 const node_t* n0 = m_facelist[i][0]; 00125 const node_t* n1 = m_facelist[i][1]; 00126 00127 while ( n0->parent() != 0 ) { 00128 m_os << n0->get_cell()->boundingBox() << std::endl; 00129 n0 = n0->parent(); 00130 } 00131 m_os << std::endl; 00132 while ( n1->parent() != 0 ) { 00133 m_os << n1->get_cell()->boundingBox() << std::endl; 00134 n1 = n1->parent(); 00135 } 00136 exit ( 0 ); 00137 } else { 00138 m_os << m_facelist[i][0]->get_cell()->boundingBox() << std::endl; 00139 m_os << m_facelist[i][1]->get_cell()->boundingBox() << std::endl; 00140 boundingBox_t c = computeCommonFace( m_facelist[i] ); 00141 m_os << "[" << c.xmin() << ", " << c.xmax() << "] x [" << c.ymin() << ", " << c.ymax() << "] x [" << c.zmin() << ", " << c.zmax() << "]\n\n" ; 00142 } 00143 00144 00145 } 00146 00147 m_os << "Errors found: " << errs << " (total " << m_facelist.size() << ")\n"; 00148 00149 } 00150 00151 template<class node_t> bool EdgeListBuilder<node_t>::commonFace 00152 ( const faceTuple_t& f ) const { 00153 00154 00155 int strict = 0; 00156 int nonstrict = 0; 00157 00158 for ( int i = 0; i < 3; ++i ) { 00159 if ( overlapAlongAxis( f[0], f[1], i, true ) ) 00160 strict++; 00161 if ( overlapAlongAxis( f[0], f[1], i, false ) ) 00162 nonstrict++; 00163 } 00164 00165 return ( strict >= 2 && nonstrict >= 3 ) ? true : false; 00166 } 00167 00168 template<class node_t> typename EdgeListBuilder<node_t>::boundingBox_t EdgeListBuilder<node_t>::computeCommonFace ( const faceTuple_t& f ) const { 00169 00170 using std::min; using std::max; 00171 00172 boundingBox_t b0 = f[0]->get_cell()->boundingBox(); 00173 boundingBox_t b1 = f[1]->get_cell()->boundingBox(); 00174 00175 boundingBox_t b( max( b0.xmin(), b1.xmin() ), 00176 min( b0.xmax(), b1.xmax() ), 00177 max( b0.ymin(), b1.ymin() ), 00178 min( b0.ymax(), b1.ymax() ), 00179 max( b0.zmin(), b1.zmin() ), 00180 min( b0.zmax(), b1.zmax() ) ); 00181 00182 return b; 00183 } 00184 00185 00186 00187 template<class node_t> bool EdgeListBuilder<node_t>::overlapAlongAxis 00188 ( const node_t* n0, const node_t* n1, const int axis, bool strict ) const { 00189 00190 // Intervals for the boundingbox 00191 00192 std::tr1::array<double, 2> a, b; 00193 00194 switch( axis ) { 00195 case 0: 00196 a[0] = n0->get_cell()->boundingBox().xmin(); 00197 b[0] = n1->get_cell()->boundingBox().xmin(); 00198 a[1] = n0->get_cell()->boundingBox().xmax(); 00199 b[1] = n1->get_cell()->boundingBox().xmax(); 00200 break; 00201 case 1: 00202 a[0] = n0->get_cell()->boundingBox().ymin(); 00203 b[0] = n1->get_cell()->boundingBox().ymin(); 00204 a[1] = n0->get_cell()->boundingBox().ymax(); 00205 b[1] = n1->get_cell()->boundingBox().ymax(); 00206 break; 00207 case 2: 00208 a[0] = n0->get_cell()->boundingBox().zmin(); 00209 b[0] = n1->get_cell()->boundingBox().zmin(); 00210 a[1] = n0->get_cell()->boundingBox().zmax(); 00211 b[1] = n1->get_cell()->boundingBox().zmax(); 00212 break; 00213 default: 00214 std::cerr << "Trying to take boudningbox along unknown axis!\n"; 00215 return false; 00216 } 00217 00218 // Now check the overlap of these intervals 00219 // Test taken from overlap method of boost::interval. 00220 00221 if ( strict ) 00222 return ( a[0] <= b[0] && b[0] < a[1] ) || 00223 ( b[0] <= a[0] && a[0] < b[1] ); 00224 else 00225 return ( a[0] <= b[0] && b[0] <= a[1] ) || 00226 ( b[0] <= a[0] && a[0] <= b[1] ); 00227 00228 } 00229 00230 template<class node_t> void EdgeListBuilder<node_t>::cellProc( const node_t* n ) { 00231 #ifdef ELB_DEBUG 00232 m_os << __FUNCTION__ << n->get_cell()->boundingBox() << std::endl; 00233 #endif 00234 00235 if ( n->is_leaf() ) 00236 return; 00237 00238 // Obvious, left and right share a face 00239 faceProc( n->left(), n->right(), n->split_dir() ); 00240 00241 // Obvious 00242 cellProc( n->left() ); 00243 cellProc( n->right() ); 00244 } 00245 00249 template<class node_t> bool EdgeListBuilder<node_t>::faceProc 00250 ( const node_t* left, const node_t* right, const int split_dir ) { 00251 00252 #ifdef ELB_DEBUG 00253 m_os << __FUNCTION__ << std::endl; 00254 // m_os << left->get_cell()->boundingBox() << std::endl; 00255 // m_os << right->get_cell()->boundingBox() << std::endl; 00256 #endif 00257 00258 // Check that boxes overlap 00259 //if ( !overlapAlongAxis( left, right, split_dir ) ) 00260 // return false; 00261 faceTuple_t f = {{left, right}}; 00262 if ( !commonFace( f ) ) 00263 return false; 00264 00265 if ( left->is_leaf() && right->is_leaf() ) { 00266 faceTuple_t f = {{ left, right }}; 00267 #ifdef ELB_DEBUG 00268 m_os << "Inserted\n"; 00269 m_os << left->get_cell()->boundingBox() << std::endl; 00270 m_os << right->get_cell()->boundingBox() << std::endl; 00271 #endif 00272 00273 #ifdef ELB_RUNTIME_CHECK 00274 if ( std::find( m_facelist.begin(), m_facelist.end(), f ) 00275 != m_facelist.end() ) { 00276 m_os << "Face already inserted!\n"; 00277 00278 exit(0); 00279 } 00280 #endif 00281 m_facelist.push_back( f ); 00282 m_facesplitdir.push_back( split_dir ); 00283 00284 return true; 00285 } 00286 00287 // We naively split into both children cells and test for 00288 // overlap at the start of the recursive step. 00289 if ( !left->is_leaf() ) { 00290 if ( left->split_dir() != split_dir ) { 00291 //edgeProc( left->left(), left->right(), right, split_dir ); 00292 faceProc( left->left(), right, split_dir ); 00293 faceProc( left->right(), right, split_dir ); 00294 } else { // left->split_dir() == split_dir 00295 faceProc( left->right(), right, split_dir ); 00296 } 00297 } else if ( !right->is_leaf() ) { 00298 if ( right->split_dir() != split_dir ) { 00299 //edgeProc( left, right->left(), right->right(), split_dir ); 00300 faceProc( left, right->left(), split_dir ); 00301 faceProc( left, right->right(), split_dir ); 00302 } else { // right->split_dir() == split_dir 00303 faceProc( left, right->left(), split_dir ); 00304 } 00305 } 00306 00307 return false; 00308 } 00309 00310 template<class node_t> void EdgeListBuilder<node_t>::edgeProc_helper 00311 ( const node_t* n0, const node_t* n1, const node_t* n2, const int split_dir ) { 00312 00313 if ( !n0->is_leaf() ) { 00314 if ( n0->split_dir() != split_dir ) { 00315 edgeProc( n0->left(), n1, n2, split_dir ); 00316 edgeProc( n0->right(), n1, n2, split_dir ); 00317 } else { // n0->split_dir() == split_dir 00318 //edgeProc( getChild( n0, getSplitDirParentType( n0, split_dir ) ), 00319 // n1, n2, split_dir ); 00320 } 00321 } 00322 } 00323 00324 template<class node_t> void EdgeListBuilder<node_t>::edgeProc 00325 ( const node_t* n0, const node_t* n1, const node_t* n2, const int split_dir ) { 00326 00327 #ifdef ELB_DEBUG 00328 m_os << __FUNCTION__ << std::endl; 00329 #endif 00330 00331 if ( n0->is_leaf() && n1->is_leaf() && n2->is_leaf() ) { 00332 edgeTuple_t n = {{ n0, n1, n2 }}; 00333 m_edgelist.push_back( n ); 00334 00335 #ifdef ELB_DEBUG 00336 m_os << "Found nodeTuple.\n"; 00337 for ( int i = 0; i < 3; ++i ) { 00338 m_os << n[i]->get_cell()->boundingBox() << std::endl; 00339 } 00340 m_os << std::endl; 00341 #endif 00342 } else { 00343 edgeProc_helper( n0, n1, n2, split_dir ); 00344 edgeProc_helper( n1, n2, n0, split_dir ); 00345 edgeProc_helper( n2, n0, n1, split_dir ); 00346 } 00347 } 00348 00349 }} // namespace mmx::shape 00350 00351 #endif // FIND_MINIMAL_EDGES_HPP
