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*/
|
---|