[857] | 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*/
|
---|