self_avoiding_walk.hpp 15 KB

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