Chris@16: //=======================================================================
Chris@16: // Copyright 1997, 1998, 1999, 2000 University of Notre Dame.
Chris@16: // Authors: Andrew Lumsdaine, Lie-Quan Lee, Jeremy G. Siek
Chris@16: //
Chris@16: // Distributed under the Boost Software License, Version 1.0. (See
Chris@16: // accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: //=======================================================================
Chris@16: #ifndef BOOST_SELF_AVOIDING_WALK_HPP
Chris@16: #define BOOST_SELF_AVOIDING_WALK_HPP
Chris@16: 
Chris@16: /*
Chris@16:   This file defines necessary components for SAW. 
Chris@16: 
Chris@16:   mesh language: (defined by myself to clearify what is what)
Chris@16:   A triangle in mesh is called an triangle.
Chris@16:   An edge in mesh is called an line. 
Chris@16:   A vertex in mesh is called a point.
Chris@16: 
Chris@16:   A triangular mesh corresponds to a graph in which a vertex is a
Chris@16:   triangle and an edge(u, v) stands for triangle u and triangle v 
Chris@16:   share an line.
Chris@16: 
Chris@16:   After this point, a vertex always refers to vertex in graph,
Chris@16:   therefore it is a traingle in mesh. 
Chris@16: 
Chris@16:  */
Chris@16: 
Chris@16: #include <utility>
Chris@16: #include <boost/config.hpp>
Chris@16: #include <boost/graph/graph_traits.hpp>
Chris@16: #include <boost/property_map/property_map.hpp>
Chris@16: 
Chris@16: #define SAW_SENTINAL -1
Chris@16: 
Chris@16: namespace boost {
Chris@16: 
Chris@16:   template <class T1, class T2, class T3>
Chris@16:   struct triple {
Chris@16:     T1 first;
Chris@16:     T2 second;
Chris@16:     T3 third;
Chris@16:     triple(const T1& a, const T2& b, const T3& c) : first(a), second(b), third(c) {}
Chris@16:     triple() : first(SAW_SENTINAL), second(SAW_SENTINAL), third(SAW_SENTINAL) {}
Chris@16:   };
Chris@16:  
Chris@16:   typedef triple<int, int, int> Triple;
Chris@16: 
Chris@16:   /* Define a vertex property which has a triangle inside. Triangle is
Chris@16:     represented by a triple.  */
Chris@16:   struct triangle_tag { enum { num = 100 }; };
Chris@16:   typedef property<triangle_tag,Triple> triangle_property;
Chris@16: 
Chris@16:   /* Define an edge property with a line. A line is represented by a
Chris@16:     pair.  This is not required for SAW though.
Chris@16:   */
Chris@16:   struct line_tag { enum { num = 101 }; };
Chris@16:   template <class T> struct line_property
Chris@16:     : public property<line_tag, std::pair<T,T> > { };
Chris@16: 
Chris@16:   /*Precondition: Points in a Triangle are in order */
Chris@16:   template <class Triangle, class Line>
Chris@16:   inline void get_sharing(const Triangle& a, const Triangle& b, Line& l)
Chris@16:   {
Chris@16:     l.first = SAW_SENTINAL;
Chris@16:     l.second = SAW_SENTINAL;
Chris@16: 
Chris@16:     if ( a.first == b.first ) {
Chris@16:       l.first = a.first;
Chris@16:       if ( a.second == b.second || a.second == b.third )
Chris@16:         l.second = a.second;
Chris@16:       else if ( a.third == b.second || a.third == b.third )
Chris@16:         l.second = a.third; 
Chris@16: 
Chris@16:     }  else if ( a.first == b.second ) {
Chris@16:       l.first = a.first;
Chris@16:       if ( a.second == b.third )
Chris@16:         l.second = a.second;
Chris@16:       else if ( a.third == b.third )
Chris@16:         l.second = a.third; 
Chris@16: 
Chris@16:     } else if ( a.first == b.third ) {
Chris@16:       l.first = a.first;
Chris@16: 
Chris@16: 
Chris@16:     } else if ( a.second == b.first ) {
Chris@16:       l.first = a.second;
Chris@16:       if ( a.third == b.second || a.third == b.third )
Chris@16:         l.second = a.third;
Chris@16: 
Chris@16:     } else if ( a.second == b.second ) {
Chris@16:       l.first = a.second;
Chris@16:       if ( a.third == b.third ) 
Chris@16:         l.second = a.third;
Chris@16: 
Chris@16:     } else if ( a.second == b.third ) {
Chris@16:       l.first = a.second;
Chris@16: 
Chris@16: 
Chris@16:     } else if ( a.third == b.first 
Chris@16:                 || a.third == b.second  
Chris@16:                 || a.third == b.third )
Chris@16:       l.first = a.third; 
Chris@16: 
Chris@16:     /*Make it in order*/
Chris@16:     if ( l.first > l.second ) {
Chris@16:       typename Line::first_type i = l.first;
Chris@16:       l.first = l.second;
Chris@16:       l.second = i;
Chris@16:     }
Chris@16: 
Chris@16:   }
Chris@16: 
Chris@16:   template <class TriangleDecorator, class Vertex, class Line>
Chris@16:   struct get_vertex_sharing {
Chris@16:     typedef std::pair<Vertex, Line> Pair;
Chris@16:     get_vertex_sharing(const TriangleDecorator& _td) : td(_td) {}
Chris@16:     inline Line operator()(const Vertex& u, const Vertex& v) const {
Chris@16:       Line l;
Chris@16:       get_sharing(td[u], td[v], l);
Chris@16:       return l;
Chris@16:     }
Chris@16:     inline Line operator()(const Pair& u, const Vertex& v) const {
Chris@16:       Line l;
Chris@16:       get_sharing(td[u.first], td[v], l);
Chris@16:       return l;
Chris@16:     }
Chris@16:     inline Line operator()(const Pair& u, const Pair& v) const {
Chris@16:       Line l;
Chris@16:       get_sharing(td[u.first], td[v.first], l);
Chris@16:       return l;
Chris@16:     }
Chris@16:     TriangleDecorator td;
Chris@16:   };
Chris@16:   
Chris@16:   /* HList has to be a handle of data holder so that pass-by-value is
Chris@16:    * in right logic.
Chris@16:    *
Chris@16:    * The element of HList is a pair of vertex and line. (remember a
Chris@16:    * line is a pair of two ints.). That indicates the walk w from
Chris@16:    * current vertex is across line. (If the first of line is -1, it is
Chris@16:    * a point though.
Chris@16:    */
Chris@16:   template < class TriangleDecorator, class HList, class IteratorD>
Chris@16:   class SAW_visitor
Chris@16:     : public bfs_visitor<>, public dfs_visitor<>
Chris@16:   {
Chris@16:     typedef typename boost::property_traits<IteratorD>::value_type iter;
Chris@16:     /*use boost shared_ptr*/
Chris@16:     typedef typename HList::element_type::value_type::second_type     Line;
Chris@16:   public:
Chris@16: 
Chris@16:     typedef tree_edge_tag category;
Chris@16: 
Chris@16:     inline SAW_visitor(TriangleDecorator _td, HList _hlist, IteratorD ia)
Chris@16:       : td(_td), hlist(_hlist), iter_d(ia) {}
Chris@16: 
Chris@16:     template <class Vertex, class Graph>
Chris@16:     inline void start_vertex(Vertex v, Graph&) {
Chris@16:       Line l1;
Chris@16:       l1.first = SAW_SENTINAL;
Chris@16:       l1.second = SAW_SENTINAL;
Chris@16:       hlist->push_front(std::make_pair(v, l1));
Chris@16:       iter_d[v] = hlist->begin();
Chris@16:     }
Chris@16: 
Chris@16:     /*Several symbols:
Chris@16:       w(i): i-th triangle in walk w
Chris@16:       w(i) |- w(i+1): w enter w(i+1) from w(i) over a line
Chris@16:       w(i) ~> w(i+1): w enter w(i+1) from w(i) over a point
Chris@16:       w(i) -> w(i+1): w enter w(i+1) from w(i)
Chris@16:       w(i) ^ w(i+1): the line or point w go over from w(i) to w(i+1)
Chris@16:     */
Chris@16:     template <class Edge, class Graph>
Chris@16:     bool tree_edge(Edge e, Graph& G) {
Chris@16:           using std::make_pair;
Chris@16:       typedef typename boost::graph_traits<Graph>::vertex_descriptor Vertex;
Chris@16:       Vertex tau = target(e, G);
Chris@16:       Vertex i   = source(e, G); 
Chris@16: 
Chris@16:       get_vertex_sharing<TriangleDecorator, Vertex, Line> get_sharing_line(td);
Chris@16:       
Chris@16:       Line tau_i = get_sharing_line(tau, i);
Chris@16:       
Chris@16:       iter w_end = hlist->end();
Chris@16: 
Chris@16:       iter w_i = iter_d[i];
Chris@16: 
Chris@16:       iter w_i_m_1 = w_i;
Chris@16:       iter w_i_p_1 = w_i;
Chris@16: 
Chris@16:       /*----------------------------------------------------------
Chris@16:        *             true             false
Chris@16:        *==========================================================
Chris@16:        *a       w(i-1) |- w(i)    w(i-1) ~> w(i) or w(i-1) is null
Chris@16:        *----------------------------------------------------------
Chris@16:        *b       w(i) |- w(i+1)    w(i) ~> w(i+1) or no w(i+1) yet
Chris@16:        *----------------------------------------------------------
Chris@16:        */
Chris@16:       
Chris@16:       bool a = false, b = false;
Chris@16: 
Chris@16:       --w_i_m_1;
Chris@16:       ++w_i_p_1;
Chris@16:       b = ( w_i->second.first != SAW_SENTINAL );
Chris@16: 
Chris@16:       if ( w_i_m_1 != w_end ) {
Chris@16:         a = ( w_i_m_1->second.first != SAW_SENTINAL );
Chris@16:       }      
Chris@16:       
Chris@16:       if ( a ) {
Chris@16:         
Chris@16:         if ( b ) {
Chris@16:           /*Case 1: 
Chris@16:             
Chris@16:             w(i-1) |- w(i) |- w(i+1)
Chris@16:           */
Chris@16:           Line l1 = get_sharing_line(*w_i_m_1, tau);
Chris@16: 
Chris@16:           iter w_i_m_2 = w_i_m_1;
Chris@16:           --w_i_m_2;
Chris@16: 
Chris@16:           bool c = true;
Chris@16:           
Chris@16:           if ( w_i_m_2 != w_end ) {
Chris@16:             c = w_i_m_2->second != l1;
Chris@16:           }
Chris@16: 
Chris@16:           if ( c ) {  /* w(i-1) ^ tau != w(i-2) ^ w(i-1)  */
Chris@16:             /*extension: w(i-1) -> tau |- w(i) */
Chris@16:             w_i_m_1->second = l1;
Chris@16:             /*insert(pos, const T&) is to insert before pos*/
Chris@16:             iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));  
Chris@16:             
Chris@16:           } else {  /* w(i-1) ^ tau == w(i-2) ^ w(i-1)  */
Chris@16:             /*must be w(i-2) ~> w(i-1) */
Chris@16:             
Chris@16:             bool d = true;
Chris@16:             //need to handle the case when w_i_p_1 is null
Chris@16:             Line l3 = get_sharing_line(*w_i_p_1, tau);
Chris@16:             if ( w_i_p_1 != w_end )
Chris@16:               d = w_i_p_1->second != l3;
Chris@16:             if ( d ) { /* w(i+1) ^ tau != w(i+1) ^ w(i+2) */
Chris@16:               /*extension: w(i) |- tau -> w(i+1) */
Chris@16:               w_i->second = tau_i;
Chris@16:               iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l3));
Chris@16:             } else { /* w(i+1) ^ tau == w(i+1) ^ w(i+2) */
Chris@16:               /*must be w(1+1) ~> w(i+2) */
Chris@16:               Line l5 = get_sharing_line(*w_i_m_1, *w_i_p_1);
Chris@16:               if ( l5 != w_i_p_1->second ) { /* w(i-1) ^ w(i+1) != w(i+1) ^ w(i+2) */
Chris@16:                 /*extension: w(i-2) -> tau |- w(i) |- w(i-1) -> w(i+1) */
Chris@16:                 w_i_m_2->second = get_sharing_line(*w_i_m_2, tau);
Chris@16:                 iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
Chris@16:                 w_i->second = w_i_m_1->second;
Chris@16:                 w_i_m_1->second = l5;
Chris@16:                 iter_d[w_i_m_1->first] = hlist->insert(w_i_p_1, *w_i_m_1);
Chris@16:                 hlist->erase(w_i_m_1);
Chris@16:               } else {
Chris@16:                 /*mesh is tetrahedral*/
Chris@16:                 // dont know what that means.
Chris@16:                 ;
Chris@16:               }
Chris@16:             }
Chris@16:             
Chris@16:           }
Chris@16:         } else { 
Chris@16:           /*Case 2:
Chris@16:             
Chris@16:             w(i-1) |- w(i) ~> w(1+1)
Chris@16:           */
Chris@16:           
Chris@16:           if ( w_i->second.second == tau_i.first 
Chris@16:                || w_i->second.second == tau_i.second ) {  /*w(i) ^ w(i+1) < w(i) ^ tau*/
Chris@16:             /*extension: w(i) |- tau -> w(i+1) */
Chris@16:             w_i->second = tau_i;
Chris@16:             Line l1 = get_sharing_line(*w_i_p_1, tau);
Chris@16:             iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l1));
Chris@16:           } else { /*w(i) ^ w(i+1) !< w(i) ^ tau*/
Chris@16:             Line l1 = get_sharing_line(*w_i_m_1, tau);
Chris@16:             bool c = true;
Chris@16:             iter w_i_m_2 = w_i_m_1;
Chris@16:             --w_i_m_2;
Chris@16:             if ( w_i_m_2 != w_end )
Chris@16:               c =  l1 != w_i_m_2->second;
Chris@16:             if (c) { /*w(i-1) ^ tau != w(i-2) ^ w(i-1)*/
Chris@16:               /*extension: w(i-1) -> tau |- w(i)*/
Chris@16:               w_i_m_1->second = l1;
Chris@16:               iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
Chris@16:             } else { /*w(i-1) ^ tau == w(i-2) ^ w(i-1)*/
Chris@16:               /*must be w(i-2)~>w(i-1)*/
Chris@16:               /*extension: w(i-2) -> tau |- w(i) |- w(i-1) -> w(i+1)*/
Chris@16:               w_i_m_2->second = get_sharing_line(*w_i_m_2, tau);
Chris@16:               iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
Chris@16:               w_i->second = w_i_m_1->second;
Chris@16:               w_i_m_1->second = get_sharing_line(*w_i_m_1, *w_i_p_1);
Chris@16:               iter_d[w_i_m_1->first] = hlist->insert(w_i_p_1, *w_i_m_1);
Chris@16:               hlist->erase(w_i_m_1);
Chris@16:             }
Chris@16:             
Chris@16:           }    
Chris@16:           
Chris@16:         }
Chris@16:         
Chris@16:       } else {
Chris@16:         
Chris@16:         if ( b ) {
Chris@16:           /*Case 3:
Chris@16:             
Chris@16:             w(i-1) ~> w(i) |- w(i+1)
Chris@16:           */
Chris@16:           bool c = false;
Chris@16:           if ( w_i_m_1 != w_end )
Chris@16:             c = ( w_i_m_1->second.second == tau_i.first)
Chris@16:               || ( w_i_m_1->second.second == tau_i.second);
Chris@16:           
Chris@16:           if ( c ) {  /*w(i-1) ^ w(i) < w(i) ^ tau*/
Chris@16:             /* extension: w(i-1) -> tau |- w(i) */
Chris@16:             if ( w_i_m_1 != w_end )
Chris@16:               w_i_m_1->second = get_sharing_line(*w_i_m_1, tau);
Chris@16:             iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
Chris@16:           } else {
Chris@16:             bool d = true;
Chris@16:             Line l1;
Chris@16:             l1.first = SAW_SENTINAL;
Chris@16:             l1.second = SAW_SENTINAL;
Chris@16:             if ( w_i_p_1 != w_end ) {
Chris@16:               l1 = get_sharing_line(*w_i_p_1, tau);
Chris@16:               d = l1 != w_i_p_1->second;
Chris@16:             }
Chris@16:             if (d) { /*w(i+1) ^ tau != w(i+1) ^ w(i+2)*/
Chris@16:               /*extension: w(i) |- tau -> w(i+1) */
Chris@16:               w_i->second = tau_i;
Chris@16:               iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l1));
Chris@16:             } else {
Chris@16:               /*must be w(i+1) ~> w(i+2)*/
Chris@16:               /*extension: w(i-1) -> w(i+1) |- w(i) |- tau -> w(i+2) */
Chris@16:               iter w_i_p_2 = w_i_p_1;
Chris@16:               ++w_i_p_2;
Chris@16:               
Chris@16:               w_i_p_1->second = w_i->second;
Chris@16:               iter_d[i] = hlist->insert(w_i_p_2, make_pair(i, tau_i));
Chris@16:               hlist->erase(w_i);
Chris@16:               Line l2 = get_sharing_line(*w_i_p_2, tau);
Chris@16:               iter_d[tau] = hlist->insert(w_i_p_2, make_pair(tau, l2));
Chris@16:             }
Chris@16:           }
Chris@16:           
Chris@16:         } else {
Chris@16:           /*Case 4:
Chris@16:             
Chris@16:             w(i-1) ~> w(i) ~> w(i+1)
Chris@16:             
Chris@16:           */
Chris@16:           bool c = false;
Chris@16:           if ( w_i_m_1 != w_end ) {
Chris@16:             c = (w_i_m_1->second.second == tau_i.first) 
Chris@16:               || (w_i_m_1->second.second == tau_i.second);
Chris@16:           }
Chris@16:           if ( c ) {   /*w(i-1) ^ w(i) < w(i) ^ tau */
Chris@16:             /*extension: w(i-1) -> tau |- w(i) */
Chris@16:             if ( w_i_m_1 != w_end )
Chris@16:               w_i_m_1->second = get_sharing_line(*w_i_m_1, tau);
Chris@16:             iter_d[tau] = hlist->insert(w_i, make_pair(tau, tau_i));
Chris@16:           } else { 
Chris@16:             /*extension: w(i) |- tau -> w(i+1) */
Chris@16:             w_i->second = tau_i;
Chris@16:             Line l1;
Chris@16:             l1.first = SAW_SENTINAL;
Chris@16:             l1.second = SAW_SENTINAL;
Chris@16:             if ( w_i_p_1 != w_end ) 
Chris@16:               l1 = get_sharing_line(*w_i_p_1, tau);
Chris@16:             iter_d[tau] = hlist->insert(w_i_p_1, make_pair(tau, l1));
Chris@16:           }
Chris@16:         }
Chris@16:         
Chris@16:       }
Chris@16: 
Chris@16:       return true;
Chris@16:     }
Chris@16:     
Chris@16:   protected:
Chris@16:     TriangleDecorator td; /*a decorator for vertex*/
Chris@16:     HList             hlist;
Chris@16:     /*This must be a handle of list to record the SAW
Chris@16:       The element type of the list is pair<Vertex, Line>
Chris@16:      */
Chris@16:    
Chris@16:     IteratorD         iter_d; 
Chris@16:     /*Problem statement: Need a fast access to w for triangle i.
Chris@16:      *Possible solution: mantain an array to record. 
Chris@16:      iter_d[i] will return an iterator 
Chris@16:      which points to w(i), where i is a vertex
Chris@16:      representing triangle i.
Chris@16:     */
Chris@16:   };
Chris@16: 
Chris@16:   template <class Triangle, class HList, class Iterator>
Chris@16:   inline 
Chris@16:   SAW_visitor<Triangle, HList, Iterator>
Chris@16:   visit_SAW(Triangle t, HList hl, Iterator i) {
Chris@16:     return SAW_visitor<Triangle, HList, Iterator>(t, hl, i);
Chris@16:   }
Chris@16: 
Chris@16:   template <class Tri, class HList, class Iter>
Chris@16:   inline 
Chris@16:   SAW_visitor< random_access_iterator_property_map<Tri*,Tri,Tri&>,
Chris@16:     HList, random_access_iterator_property_map<Iter*,Iter,Iter&> >
Chris@16:   visit_SAW_ptr(Tri* t, HList hl, Iter* i) {
Chris@16:     typedef random_access_iterator_property_map<Tri*,Tri,Tri&> TriD;
Chris@16:     typedef random_access_iterator_property_map<Iter*,Iter,Iter&> IterD;
Chris@16:     return SAW_visitor<TriD, HList, IterD>(t, hl, i);
Chris@16:   }
Chris@16: 
Chris@16:   // should also have combo's of pointers, and also const :(
Chris@16:   
Chris@16: }
Chris@16: 
Chris@16: #endif /*BOOST_SAW_H*/