self_avoiding_walk.hpp 14 KB

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