| 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.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; |
|---|
| 43 | triple(const T1& a, const T2& b, const T3& c) : first(a), second(b), third(c) {} |
|---|
| 44 | 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> |
|---|
| 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; |
|---|
| 117 | get_vertex_sharing(const TriangleDecorator& _td) : td(_td) {} |
|---|
| 118 | 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 | } |
|---|
| 123 | 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 | } |
|---|
| 128 | 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 | |
|---|
| 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> |
|---|
| 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> |
|---|
| 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> |
|---|
| 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&> > |
|---|
| 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*/ |
|---|