graph/detail/self_avoiding_walk.hpp

0001 //=======================================================================
0002 // Copyright 1997, 1998, 1999, 2000 University of Notre Dame.
0003 // Authors: Andrew Lumsdaine, Lie-Quan Lee, Jeremy G. Siek
0004 //
0005 // Distributed under the Boost Software License, Version 1.0. (See
0006 // accompanying file LICENSE_1_0.txt or copy at
0007 // http://www.boost.org/LICENSE_1_0.txt)
0008 //=======================================================================
0009 #ifndef BOOST_SELF_AVOIDING_WALK_HPP
0010 #define BOOST_SELF_AVOIDING_WALK_HPP
0011
0012 /*
0013   This file defines necessary components for SAW.
0014
0015   mesh language: (defined by myself to clearify what is what)
0016   A triangle in mesh is called an triangle.
0017   An edge in mesh is called an line.
0018   A vertex in mesh is called a point.
0019
0020   A triangular mesh corresponds to a graph in which a vertex is a
0021   triangle and an edge(u, v) stands for triangle u and triangle v
0022   share an line.
0023
0024   After this point, a vertex always refers to vertex in graph,
0025   therefore it is a traingle in mesh.
0026
0027  */
0028
0029 #include <utility>
0030 #include <boost/config.hpp>
0031 #include <boost/graph/graph_traits.hpp>
0032 #include <boost/property_map/property_map.hpp>
0033
0034 #define SAW_SENTINAL -1
0035
0036 namespace boost
0037 {
0038
0039 template < class T1, class T2, class T3 > struct triple
0040 {
0041     T1 first;
0042     T2 second;
0043     T3 third;
0044     triple(const T1& a, const T2& b, const T3& c)
0045     : first(a), second(b), third(c)
0046     {
0047     }
0048     triple() : first(SAW_SENTINAL), second(SAW_SENTINAL), third(SAW_SENTINAL) {}
0049 };
0050
0051 typedef triple< int, int, int > Triple;
0052
0053 /* Define a vertex property which has a triangle inside. Triangle is
0054   represented by a triple.  */
0055 struct triangle_tag
0056 {
0057     enum
0058     {
0059         num = 100
0060     };
0061 };
0062 typedef property< triangle_tag, Triple > triangle_property;
0063
0064 /* Define an edge property with a line. A line is represented by a
0065   pair.  This is not required for SAW though.
0066 */
0067 struct line_tag
0068 {
0069     enum
0070     {
0071         num = 101
0072     };
0073 };
0074 template < class T >
0075 struct line_property : public property< line_tag, std::pair< T, T > >
0076 {
0077 };
0078
0079 /*Precondition: Points in a Triangle are in order */
0080 template < class Triangle, class Line >
0081 inline void get_sharing(const Triangle& a, const Triangle& b, Line& l)
0082 {
0083     l.first = SAW_SENTINAL;
0084     l.second = SAW_SENTINAL;
0085
0086     if (a.first == b.first)
0087     {
0088         l.first = a.first;
0089         if (a.second == b.second || a.second == b.third)
0090             l.second = a.second;
0091         else if (a.third == b.second || a.third == b.third)
0092             l.second = a.third;
0093     }
0094     else if (a.first == b.second)
0095     {
0096         l.first = a.first;
0097         if (a.second == b.third)
0098             l.second = a.second;
0099         else if (a.third == b.third)
0100             l.second = a.third;
0101     }
0102     else if (a.first == b.third)
0103     {
0104         l.first = a.first;
0105     }
0106     else if (a.second == b.first)
0107     {
0108         l.first = a.second;
0109         if (a.third == b.second || a.third == b.third)
0110             l.second = a.third;
0111     }
0112     else if (a.second == b.second)
0113     {
0114         l.first = a.second;
0115         if (a.third == b.third)
0116             l.second = a.third;
0117     }
0118     else if (a.second == b.third)
0119     {
0120         l.first = a.second;
0121     }
0122     else if (a.third == b.first || a.third == b.second || a.third == b.third)
0123         l.first = a.third;
0124
0125     /*Make it in order*/
0126     if (l.first > l.second)
0127     {
0128         typename Line::first_type i = l.first;
0129         l.first = l.second;
0130         l.second = i;
0131     }
0132 }
0133
0134 template < class TriangleDecorator, class Vertex, class Line >
0135 struct get_vertex_sharing
0136 {
0137     typedef std::pair< Vertex, Line > Pair;
0138     get_vertex_sharing(const TriangleDecorator& _td) : td(_td) {}
0139     inline Line operator()(const Vertex& u, const Vertex& v) const
0140     {
0141         Line l;
0142         get_sharing(td[u], td[v], l);
0143         return l;
0144     }
0145     inline Line operator()(const Pair& u, const Vertex& v) const
0146     {
0147         Line l;
0148         get_sharing(td[u.first], td[v], l);
0149         return l;
0150     }
0151     inline Line operator()(const Pair& u, const Pair& v) const
0152     {
0153         Line l;
0154         get_sharing(td[u.first], td[v.first], l);
0155         return l;
0156     }
0157     TriangleDecorator td;
0158 };
0159
0160 /* HList has to be a handle of data holder so that pass-by-value is
0161  * in right logic.
0162  *
0163  * The element of HList is a pair of vertex and line. (remember a
0164  * line is a pair of two ints.). That indicates the walk w from
0165  * current vertex is across line. (If the first of line is -1, it is
0166  * a point though.
0167  */
0168 template < class TriangleDecorator, class HList, class IteratorD >
0169 class SAW_visitor : public bfs_visitor<>, public dfs_visitor<>
0170 {
0171     typedef typename boost::property_traits< IteratorD >::value_type iter;
0172     /*use boost shared_ptr*/
0173     typedef typename HList::element_type::value_type::second_type Line;
0174
0175 public:
0176     typedef tree_edge_tag category;
0177
0178     inline SAW_visitor(TriangleDecorator _td, HList _hlist, IteratorD ia)
0179     : td(_td), hlist(_hlist), iter_d(ia)
0180     {
0181     }
0182
0183     template < class Vertex, class Graph >
0184     inline void start_vertex(Vertex v, Graph&)
0185     {
0186         Line l1;
0187         l1.first = SAW_SENTINAL;
0188         l1.second = SAW_SENTINAL;
0189         hlist->push_front(std::make_pair(v, l1));
0190         iter_d[v] = hlist->begin();
0191     }
0192
0193     /*Several symbols:
0194       w(i): i-th triangle in walk w
0195       w(i) |- w(i+1): w enter w(i+1) from w(i) over a line
0196       w(i) ~> w(i+1): w enter w(i+1) from w(i) over a point
0197       w(i) -> w(i+1): w enter w(i+1) from w(i)
0198       w(i) ^ w(i+1): the line or point w go over from w(i) to w(i+1)
0199     */
0200     template < class Edge, class Graph > bool tree_edge(Edge e, Graph& G)
0201     {
0202         using std::make_pair;
0203         typedef typename boost::graph_traits< Graph >::vertex_descriptor Vertex;
0204         Vertex tau = target(e, G);
0205         Vertex i = source(e, G);
0206
0207         get_vertex_sharing< TriangleDecorator, Vertex, Line > get_sharing_line(
0208             td);
0209
0210         Line tau_i = get_sharing_line(tau, i);
0211
0212         iter w_end = hlist->end();
0213
0214         iter w_i = iter_d[i];
0215
0216         iter w_i_m_1 = w_i;
0217         iter w_i_p_1 = w_i;
0218
0219         /*----------------------------------------------------------
0220          *             true             false
0221          *==========================================================
0222          *a       w(i-1) |- w(i)    w(i-1) ~> w(i) or w(i-1) is null
0223          *----------------------------------------------------------
0224          *b       w(i) |- w(i+1)    w(i) ~> w(i+1) or no w(i+1) yet
0225          *----------------------------------------------------------
0226          */
0227
0228         bool a = false, b = false;
0229
0230         --w_i_m_1;
0231         ++w_i_p_1;
0232         b = (w_i->second.first != SAW_SENTINAL);
0233
0234         if (w_i_m_1 != w_end)
0235         {
0236             a = (w_i_m_1->second.first != SAW_SENTINAL);
0237         }
0238
0239         if (a)
0240         {
0241
0242             if (b)
0243             {
0244                 /*Case 1:
0245
0246                   w(i-1) |- w(i) |- w(i+1)
0247                 */
0248                 Line l1 = get_sharing_line(*w_i_m_1, tau);
0249
0250                 iter w_i_m_2 = w_i_m_1;
0251                 --w_i_m_2;
0252
0253                 bool c = true;
0254
0255                 if (w_i_m_2 != w_end)
0256                 {
0257                     c = w_i_m_2->second != l1;
0258                 }
0259
0260                 if (c)
0261                 { /* w(i-1) ^ tau != w(i-2) ^ w(i-1)  */
0262                     /*extension: w(i-1) -> tau |- w(i) */
0263                     w_i_m_1->second = l1;
0264                     /*insert(pos, const T&) is to insert before pos*/
0265                     iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
0266                 }
0267                 else
0268                 { /* w(i-1) ^ tau == w(i-2) ^ w(i-1)  */
0269                     /*must be w(i-2) ~> w(i-1) */
0270
0271                     bool d = true;
0272                     // need to handle the case when w_i_p_1 is null
0273                     Line l3 = get_sharing_line(*w_i_p_1, tau);
0274                     if (w_i_p_1 != w_end)
0275                         d = w_i_p_1->second != l3;
0276                     if (d)
0277                     { /* w(i+1) ^ tau != w(i+1) ^ w(i+2) */
0278                         /*extension: w(i) |- tau -> w(i+1) */
0279                         w_i->second = tau_i;
0280                         iter_d[tau]
0281                             = hlist->insert(w_i_p_1, make_pair(tau, l3));
0282                     }
0283                     else
0284                     { /* w(i+1) ^ tau == w(i+1) ^ w(i+2) */
0285                         /*must be w(1+1) ~> w(i+2) */
0286                         Line l5 = get_sharing_line(*w_i_m_1, *w_i_p_1);
0287                         if (l5 != w_i_p_1->second)
0288                         { /* w(i-1) ^ w(i+1) != w(i+1) ^ w(i+2) */
0289                             /*extension: w(i-2) -> tau |- w(i) |- w(i-1) ->
0290                              * w(i+1) */
0291                             w_i_m_2->second = get_sharing_line(*w_i_m_2, tau);
0292                             iter_d[tau]
0293                                 = hlist->insert(w_i, make_pair(tau, tau_i));
0294                             w_i->second = w_i_m_1->second;
0295                             w_i_m_1->second = l5;
0296                             iter_d[w_i_m_1->first]
0297                                 = hlist->insert(w_i_p_1, *w_i_m_1);
0298                             hlist->erase(w_i_m_1);
0299                         }
0300                         else
0301                         {
0302                             /*mesh is tetrahedral*/
0303                             // dont know what that means.
0304                             ;
0305                         }
0306                     }
0307                 }
0308             }
0309             else
0310             {
0311                 /*Case 2:
0312
0313                   w(i-1) |- w(i) ~> w(1+1)
0314                 */
0315
0316                 if (w_i->second.second == tau_i.first
0317                     || w_i->second.second == tau_i.second)
0318                 { /*w(i) ^ w(i+1) < w(i) ^ tau*/
0319                     /*extension: w(i) |- tau -> w(i+1) */
0320                     w_i->second = tau_i;
0321                     Line l1 = get_sharing_line(*w_i_p_1, tau);
0322                     iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l1));
0323                 }
0324                 else
0325                 { /*w(i) ^ w(i+1) !< w(i) ^ tau*/
0326                     Line l1 = get_sharing_line(*w_i_m_1, tau);
0327                     bool c = true;
0328                     iter w_i_m_2 = w_i_m_1;
0329                     --w_i_m_2;
0330                     if (w_i_m_2 != w_end)
0331                         c = l1 != w_i_m_2->second;
0332                     if (c)
0333                     { /*w(i-1) ^ tau != w(i-2) ^ w(i-1)*/
0334                         /*extension: w(i-1) -> tau |- w(i)*/
0335                         w_i_m_1->second = l1;
0336                         iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
0337                     }
0338                     else
0339                     { /*w(i-1) ^ tau == w(i-2) ^ w(i-1)*/
0340                         /*must be w(i-2)~>w(i-1)*/
0341                         /*extension: w(i-2) -> tau |- w(i) |- w(i-1) -> w(i+1)*/
0342                         w_i_m_2->second = get_sharing_line(*w_i_m_2, tau);
0343                         iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
0344                         w_i->second = w_i_m_1->second;
0345                         w_i_m_1->second = get_sharing_line(*w_i_m_1, *w_i_p_1);
0346                         iter_d[w_i_m_1->first]
0347                             = hlist->insert(w_i_p_1, *w_i_m_1);
0348                         hlist->erase(w_i_m_1);
0349                     }
0350                 }
0351             }
0352         }
0353         else
0354         {
0355
0356             if (b)
0357             {
0358                 /*Case 3:
0359
0360                   w(i-1) ~> w(i) |- w(i+1)
0361                 */
0362                 bool c = false;
0363                 if (w_i_m_1 != w_end)
0364                     c = (w_i_m_1->second.second == tau_i.first)
0365                         || (w_i_m_1->second.second == tau_i.second);
0366
0367                 if (c)
0368                 { /*w(i-1) ^ w(i) < w(i) ^ tau*/
0369                     /* extension: w(i-1) -> tau |- w(i) */
0370                     if (w_i_m_1 != w_end)
0371                         w_i_m_1->second = get_sharing_line(*w_i_m_1, tau);
0372                     iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
0373                 }
0374                 else
0375                 {
0376                     bool d = true;
0377                     Line l1;
0378                     l1.first = SAW_SENTINAL;
0379                     l1.second = SAW_SENTINAL;
0380                     if (w_i_p_1 != w_end)
0381                     {
0382                         l1 = get_sharing_line(*w_i_p_1, tau);
0383                         d = l1 != w_i_p_1->second;
0384                     }
0385                     if (d)
0386                     { /*w(i+1) ^ tau != w(i+1) ^ w(i+2)*/
0387                         /*extension: w(i) |- tau -> w(i+1) */
0388                         w_i->second = tau_i;
0389                         iter_d[tau]
0390                             = hlist->insert(w_i_p_1, make_pair(tau, l1));
0391                     }
0392                     else
0393                     {
0394                         /*must be w(i+1) ~> w(i+2)*/
0395                         /*extension: w(i-1) -> w(i+1) |- w(i) |- tau -> w(i+2)
0396                          */
0397                         iter w_i_p_2 = w_i_p_1;
0398                         ++w_i_p_2;
0399
0400                         w_i_p_1->second = w_i->second;
0401                         iter_d[i] = hlist->insert(w_i_p_2, make_pair(i, tau_i));
0402                         hlist->erase(w_i);
0403                         Line l2 = get_sharing_line(*w_i_p_2, tau);
0404                         iter_d[tau]
0405                             = hlist->insert(w_i_p_2, make_pair(tau, l2));
0406                     }
0407                 }
0408             }
0409             else
0410             {
0411                 /*Case 4:
0412
0413                   w(i-1) ~> w(i) ~> w(i+1)
0414
0415                 */
0416                 bool c = false;
0417                 if (w_i_m_1 != w_end)
0418                 {
0419                     c = (w_i_m_1->second.second == tau_i.first)
0420                         || (w_i_m_1->second.second == tau_i.second);
0421                 }
0422                 if (c)
0423                 { /*w(i-1) ^ w(i) < w(i) ^ tau */
0424                     /*extension: w(i-1) -> tau |- w(i) */
0425                     if (w_i_m_1 != w_end)
0426                         w_i_m_1->second = get_sharing_line(*w_i_m_1, tau);
0427                     iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
0428                 }
0429                 else
0430                 {
0431                     /*extension: w(i) |- tau -> w(i+1) */
0432                     w_i->second = tau_i;
0433                     Line l1;
0434                     l1.first = SAW_SENTINAL;
0435                     l1.second = SAW_SENTINAL;
0436                     if (w_i_p_1 != w_end)
0437                         l1 = get_sharing_line(*w_i_p_1, tau);
0438                     iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l1));
0439                 }
0440             }
0441         }
0442
0443         return true;
0444     }
0445
0446 protected:
0447     TriangleDecorator td; /*a decorator for vertex*/
0448     HList hlist;
0449     /*This must be a handle of list to record the SAW
0450       The element type of the list is pair<Vertex, Line>
0451      */
0452
0453     IteratorD iter_d;
0454     /*Problem statement: Need a fast access to w for triangle i.
0455      *Possible solution: mantain an array to record.
0456      iter_d[i] will return an iterator
0457      which points to w(i), where i is a vertex
0458      representing triangle i.
0459     */
0460 };
0461
0462 template < class Triangle, class HList, class Iterator >
0463 inline SAW_visitor< Triangle, HList, Iterator > visit_SAW(
0464     Triangle t, HList hl, Iterator i)
0465 {
0466     return SAW_visitor< Triangle, HList, Iterator >(t, hl, i);
0467 }
0468
0469 template < class Tri, class HList, class Iter >
0470 inline SAW_visitor< random_access_iterator_property_map< Tri*, Tri, Tri& >,
0471     HList, random_access_iterator_property_map< Iter*, Iter, Iter& > >
0472 visit_SAW_ptr(Tri* t, HList hl, Iter* i)
0473 {
0474     typedef random_access_iterator_property_map< Tri*, Tri, Tri& > TriD;
0475     typedef random_access_iterator_property_map< Iter*, Iter, Iter& > IterD;
0476     return SAW_visitor< TriD, HList, IterD >(t, hl, i);
0477 }
0478
0479 // should also have combo's of pointers, and also const :(
0480
0481 }
0482
0483 #endif /*BOOST_SAW_H*/