1 //======================================================================= 2 // Copyright 1997, 1998, 1999, 2000 University of Notre Dame. 3 // Authors: Andrew Lumsdaine, Lie-Quan Lee, Jeremy G. Siek 4 // 5 // Distributed under the Boost Software License, Version 1.0. (See 6 // accompanying file LICENSE_1_0.txt or copy at 7 // http://www.boost.org/LICENSE_1_0.txt) 8 //======================================================================= 9 #ifndef BOOST_SELF_AVOIDING_WALK_HPP 10 #define BOOST_SELF_AVOIDING_WALK_HPP 11 12 /* 13 This file defines necessary components for SAW. 14 15 mesh language: (defined by myself to clearify what is what) 16 A triangle in mesh is called an triangle. 17 An edge in mesh is called an line. 18 A vertex in mesh is called a point. 19 20 A triangular mesh corresponds to a graph in which a vertex is a 21 triangle and an edge(u, v) stands for triangle u and triangle v 22 share an line. 23 24 After this point, a vertex always refers to vertex in graph, 25 therefore it is a traingle in mesh. 26 27 */ 28 29 #include <utility> 30 #include <boost/config.hpp> 31 #include <boost/graph/graph_traits.hpp> 32 #include <boost/property_map/property_map.hpp> 33 34 #define SAW_SENTINAL -1 35 36 namespace boost { 37 38 template <class T1, class T2, class T3> 39 struct triple { 40 T1 first; 41 T2 second; 42 T3 third; tripleboost::triple43 triple(const T1& a, const T2& b, const T3& c) : first(a), second(b), third(c) {} tripleboost::triple44 triple() : first(SAW_SENTINAL), second(SAW_SENTINAL), third(SAW_SENTINAL) {} 45 }; 46 47 typedef triple<int, int, int> Triple; 48 49 /* Define a vertex property which has a triangle inside. Triangle is 50 represented by a triple. */ 51 struct triangle_tag { enum { num = 100 }; }; 52 typedef property<triangle_tag,Triple> triangle_property; 53 54 /* Define an edge property with a line. A line is represented by a 55 pair. This is not required for SAW though. 56 */ 57 struct line_tag { enum { num = 101 }; }; 58 template <class T> struct line_property 59 : public property<line_tag, std::pair<T,T> > { }; 60 61 /*Precondition: Points in a Triangle are in order */ 62 template <class Triangle, class Line> get_sharing(const Triangle & a,const Triangle & b,Line & l)63 inline void get_sharing(const Triangle& a, const Triangle& b, Line& l) 64 { 65 l.first = SAW_SENTINAL; 66 l.second = SAW_SENTINAL; 67 68 if ( a.first == b.first ) { 69 l.first = a.first; 70 if ( a.second == b.second || a.second == b.third ) 71 l.second = a.second; 72 else if ( a.third == b.second || a.third == b.third ) 73 l.second = a.third; 74 75 } else if ( a.first == b.second ) { 76 l.first = a.first; 77 if ( a.second == b.third ) 78 l.second = a.second; 79 else if ( a.third == b.third ) 80 l.second = a.third; 81 82 } else if ( a.first == b.third ) { 83 l.first = a.first; 84 85 86 } else if ( a.second == b.first ) { 87 l.first = a.second; 88 if ( a.third == b.second || a.third == b.third ) 89 l.second = a.third; 90 91 } else if ( a.second == b.second ) { 92 l.first = a.second; 93 if ( a.third == b.third ) 94 l.second = a.third; 95 96 } else if ( a.second == b.third ) { 97 l.first = a.second; 98 99 100 } else if ( a.third == b.first 101 || a.third == b.second 102 || a.third == b.third ) 103 l.first = a.third; 104 105 /*Make it in order*/ 106 if ( l.first > l.second ) { 107 typename Line::first_type i = l.first; 108 l.first = l.second; 109 l.second = i; 110 } 111 112 } 113 114 template <class TriangleDecorator, class Vertex, class Line> 115 struct get_vertex_sharing { 116 typedef std::pair<Vertex, Line> Pair; get_vertex_sharingboost::get_vertex_sharing117 get_vertex_sharing(const TriangleDecorator& _td) : td(_td) {} operator ()boost::get_vertex_sharing118 inline Line operator()(const Vertex& u, const Vertex& v) const { 119 Line l; 120 get_sharing(td[u], td[v], l); 121 return l; 122 } operator ()boost::get_vertex_sharing123 inline Line operator()(const Pair& u, const Vertex& v) const { 124 Line l; 125 get_sharing(td[u.first], td[v], l); 126 return l; 127 } operator ()boost::get_vertex_sharing128 inline Line operator()(const Pair& u, const Pair& v) const { 129 Line l; 130 get_sharing(td[u.first], td[v.first], l); 131 return l; 132 } 133 TriangleDecorator td; 134 }; 135 136 /* HList has to be a handle of data holder so that pass-by-value is 137 * in right logic. 138 * 139 * The element of HList is a pair of vertex and line. (remember a 140 * line is a pair of two ints.). That indicates the walk w from 141 * current vertex is across line. (If the first of line is -1, it is 142 * a point though. 143 */ 144 template < class TriangleDecorator, class HList, class IteratorD> 145 class SAW_visitor 146 : public bfs_visitor<>, public dfs_visitor<> 147 { 148 typedef typename boost::property_traits<IteratorD>::value_type iter; 149 /*use boost shared_ptr*/ 150 typedef typename HList::element_type::value_type::second_type Line; 151 public: 152 153 typedef tree_edge_tag category; 154 SAW_visitor(TriangleDecorator _td,HList _hlist,IteratorD ia)155 inline SAW_visitor(TriangleDecorator _td, HList _hlist, IteratorD ia) 156 : td(_td), hlist(_hlist), iter_d(ia) {} 157 158 template <class Vertex, class Graph> start_vertex(Vertex v,Graph &)159 inline void start_vertex(Vertex v, Graph&) { 160 Line l1; 161 l1.first = SAW_SENTINAL; 162 l1.second = SAW_SENTINAL; 163 hlist->push_front(std::make_pair(v, l1)); 164 iter_d[v] = hlist->begin(); 165 } 166 167 /*Several symbols: 168 w(i): i-th triangle in walk w 169 w(i) |- w(i+1): w enter w(i+1) from w(i) over a line 170 w(i) ~> w(i+1): w enter w(i+1) from w(i) over a point 171 w(i) -> w(i+1): w enter w(i+1) from w(i) 172 w(i) ^ w(i+1): the line or point w go over from w(i) to w(i+1) 173 */ 174 template <class Edge, class Graph> tree_edge(Edge e,Graph & G)175 bool tree_edge(Edge e, Graph& G) { 176 using std::make_pair; 177 typedef typename boost::graph_traits<Graph>::vertex_descriptor Vertex; 178 Vertex tau = target(e, G); 179 Vertex i = source(e, G); 180 181 get_vertex_sharing<TriangleDecorator, Vertex, Line> get_sharing_line(td); 182 183 Line tau_i = get_sharing_line(tau, i); 184 185 iter w_end = hlist->end(); 186 187 iter w_i = iter_d[i]; 188 189 iter w_i_m_1 = w_i; 190 iter w_i_p_1 = w_i; 191 192 /*---------------------------------------------------------- 193 * true false 194 *========================================================== 195 *a w(i-1) |- w(i) w(i-1) ~> w(i) or w(i-1) is null 196 *---------------------------------------------------------- 197 *b w(i) |- w(i+1) w(i) ~> w(i+1) or no w(i+1) yet 198 *---------------------------------------------------------- 199 */ 200 201 bool a = false, b = false; 202 203 --w_i_m_1; 204 ++w_i_p_1; 205 b = ( w_i->second.first != SAW_SENTINAL ); 206 207 if ( w_i_m_1 != w_end ) { 208 a = ( w_i_m_1->second.first != SAW_SENTINAL ); 209 } 210 211 if ( a ) { 212 213 if ( b ) { 214 /*Case 1: 215 216 w(i-1) |- w(i) |- w(i+1) 217 */ 218 Line l1 = get_sharing_line(*w_i_m_1, tau); 219 220 iter w_i_m_2 = w_i_m_1; 221 --w_i_m_2; 222 223 bool c = true; 224 225 if ( w_i_m_2 != w_end ) { 226 c = w_i_m_2->second != l1; 227 } 228 229 if ( c ) { /* w(i-1) ^ tau != w(i-2) ^ w(i-1) */ 230 /*extension: w(i-1) -> tau |- w(i) */ 231 w_i_m_1->second = l1; 232 /*insert(pos, const T&) is to insert before pos*/ 233 iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i)); 234 235 } else { /* w(i-1) ^ tau == w(i-2) ^ w(i-1) */ 236 /*must be w(i-2) ~> w(i-1) */ 237 238 bool d = true; 239 //need to handle the case when w_i_p_1 is null 240 Line l3 = get_sharing_line(*w_i_p_1, tau); 241 if ( w_i_p_1 != w_end ) 242 d = w_i_p_1->second != l3; 243 if ( d ) { /* w(i+1) ^ tau != w(i+1) ^ w(i+2) */ 244 /*extension: w(i) |- tau -> w(i+1) */ 245 w_i->second = tau_i; 246 iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l3)); 247 } else { /* w(i+1) ^ tau == w(i+1) ^ w(i+2) */ 248 /*must be w(1+1) ~> w(i+2) */ 249 Line l5 = get_sharing_line(*w_i_m_1, *w_i_p_1); 250 if ( l5 != w_i_p_1->second ) { /* w(i-1) ^ w(i+1) != w(i+1) ^ w(i+2) */ 251 /*extension: w(i-2) -> tau |- w(i) |- w(i-1) -> w(i+1) */ 252 w_i_m_2->second = get_sharing_line(*w_i_m_2, tau); 253 iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i)); 254 w_i->second = w_i_m_1->second; 255 w_i_m_1->second = l5; 256 iter_d[w_i_m_1->first] = hlist->insert(w_i_p_1, *w_i_m_1); 257 hlist->erase(w_i_m_1); 258 } else { 259 /*mesh is tetrahedral*/ 260 // dont know what that means. 261 ; 262 } 263 } 264 265 } 266 } else { 267 /*Case 2: 268 269 w(i-1) |- w(i) ~> w(1+1) 270 */ 271 272 if ( w_i->second.second == tau_i.first 273 || w_i->second.second == tau_i.second ) { /*w(i) ^ w(i+1) < w(i) ^ tau*/ 274 /*extension: w(i) |- tau -> w(i+1) */ 275 w_i->second = tau_i; 276 Line l1 = get_sharing_line(*w_i_p_1, tau); 277 iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l1)); 278 } else { /*w(i) ^ w(i+1) !< w(i) ^ tau*/ 279 Line l1 = get_sharing_line(*w_i_m_1, tau); 280 bool c = true; 281 iter w_i_m_2 = w_i_m_1; 282 --w_i_m_2; 283 if ( w_i_m_2 != w_end ) 284 c = l1 != w_i_m_2->second; 285 if (c) { /*w(i-1) ^ tau != w(i-2) ^ w(i-1)*/ 286 /*extension: w(i-1) -> tau |- w(i)*/ 287 w_i_m_1->second = l1; 288 iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i)); 289 } else { /*w(i-1) ^ tau == w(i-2) ^ w(i-1)*/ 290 /*must be w(i-2)~>w(i-1)*/ 291 /*extension: w(i-2) -> tau |- w(i) |- w(i-1) -> w(i+1)*/ 292 w_i_m_2->second = get_sharing_line(*w_i_m_2, tau); 293 iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i)); 294 w_i->second = w_i_m_1->second; 295 w_i_m_1->second = get_sharing_line(*w_i_m_1, *w_i_p_1); 296 iter_d[w_i_m_1->first] = hlist->insert(w_i_p_1, *w_i_m_1); 297 hlist->erase(w_i_m_1); 298 } 299 300 } 301 302 } 303 304 } else { 305 306 if ( b ) { 307 /*Case 3: 308 309 w(i-1) ~> w(i) |- w(i+1) 310 */ 311 bool c = false; 312 if ( w_i_m_1 != w_end ) 313 c = ( w_i_m_1->second.second == tau_i.first) 314 || ( w_i_m_1->second.second == tau_i.second); 315 316 if ( c ) { /*w(i-1) ^ w(i) < w(i) ^ tau*/ 317 /* extension: w(i-1) -> tau |- w(i) */ 318 if ( w_i_m_1 != w_end ) 319 w_i_m_1->second = get_sharing_line(*w_i_m_1, tau); 320 iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i)); 321 } else { 322 bool d = true; 323 Line l1; 324 l1.first = SAW_SENTINAL; 325 l1.second = SAW_SENTINAL; 326 if ( w_i_p_1 != w_end ) { 327 l1 = get_sharing_line(*w_i_p_1, tau); 328 d = l1 != w_i_p_1->second; 329 } 330 if (d) { /*w(i+1) ^ tau != w(i+1) ^ w(i+2)*/ 331 /*extension: w(i) |- tau -> w(i+1) */ 332 w_i->second = tau_i; 333 iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l1)); 334 } else { 335 /*must be w(i+1) ~> w(i+2)*/ 336 /*extension: w(i-1) -> w(i+1) |- w(i) |- tau -> w(i+2) */ 337 iter w_i_p_2 = w_i_p_1; 338 ++w_i_p_2; 339 340 w_i_p_1->second = w_i->second; 341 iter_d[i] = hlist->insert(w_i_p_2, make_pair(i, tau_i)); 342 hlist->erase(w_i); 343 Line l2 = get_sharing_line(*w_i_p_2, tau); 344 iter_d[tau] = hlist->insert(w_i_p_2, make_pair(tau, l2)); 345 } 346 } 347 348 } else { 349 /*Case 4: 350 351 w(i-1) ~> w(i) ~> w(i+1) 352 353 */ 354 bool c = false; 355 if ( w_i_m_1 != w_end ) { 356 c = (w_i_m_1->second.second == tau_i.first) 357 || (w_i_m_1->second.second == tau_i.second); 358 } 359 if ( c ) { /*w(i-1) ^ w(i) < w(i) ^ tau */ 360 /*extension: w(i-1) -> tau |- w(i) */ 361 if ( w_i_m_1 != w_end ) 362 w_i_m_1->second = get_sharing_line(*w_i_m_1, tau); 363 iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i)); 364 } else { 365 /*extension: w(i) |- tau -> w(i+1) */ 366 w_i->second = tau_i; 367 Line l1; 368 l1.first = SAW_SENTINAL; 369 l1.second = SAW_SENTINAL; 370 if ( w_i_p_1 != w_end ) 371 l1 = get_sharing_line(*w_i_p_1, tau); 372 iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l1)); 373 } 374 } 375 376 } 377 378 return true; 379 } 380 381 protected: 382 TriangleDecorator td; /*a decorator for vertex*/ 383 HList hlist; 384 /*This must be a handle of list to record the SAW 385 The element type of the list is pair<Vertex, Line> 386 */ 387 388 IteratorD iter_d; 389 /*Problem statement: Need a fast access to w for triangle i. 390 *Possible solution: mantain an array to record. 391 iter_d[i] will return an iterator 392 which points to w(i), where i is a vertex 393 representing triangle i. 394 */ 395 }; 396 397 template <class Triangle, class HList, class Iterator> 398 inline 399 SAW_visitor<Triangle, HList, Iterator> visit_SAW(Triangle t,HList hl,Iterator i)400 visit_SAW(Triangle t, HList hl, Iterator i) { 401 return SAW_visitor<Triangle, HList, Iterator>(t, hl, i); 402 } 403 404 template <class Tri, class HList, class Iter> 405 inline 406 SAW_visitor< random_access_iterator_property_map<Tri*,Tri,Tri&>, 407 HList, random_access_iterator_property_map<Iter*,Iter,Iter&> > visit_SAW_ptr(Tri * t,HList hl,Iter * i)408 visit_SAW_ptr(Tri* t, HList hl, Iter* i) { 409 typedef random_access_iterator_property_map<Tri*,Tri,Tri&> TriD; 410 typedef random_access_iterator_property_map<Iter*,Iter,Iter&> IterD; 411 return SAW_visitor<TriD, HList, IterD>(t, hl, i); 412 } 413 414 // should also have combo's of pointers, and also const :( 415 416 } 417 418 #endif /*BOOST_SAW_H*/ 419