A.I.M.S algorithms


geodesic_algorithm_graph_base.h
Go to the documentation of this file.
1 #ifndef GEODESIC_ALGORITHM_GRAPH_BASE_010907
2 #define GEODESIC_ALGORITHM_GRAPH_BASE_010907
3 
6 #include <vector>
7 #include <set>
8 #include <assert.h>
9 
10 namespace geodesic{
11 
12 template<class Node>
14 {
15 public:
16  typedef Node* node_pointer;
17 
20  {};
21 
23 
24  void propagate(std::vector<SurfacePoint>& sources,
25  double max_propagation_distance = GEODESIC_INF, //propagation algorithm stops after reaching the certain distance from the source
26  std::vector<SurfacePoint>* stop_points = NULL); //or after ensuring that all the stop_points are covered
27 
28  void trace_back(SurfacePoint& destination, //trace back piecewise-linear path
29  std::vector<SurfacePoint>& path);
30 
31  void trace_back_with_index(SurfacePoint& destination,
32  std::vector<SurfacePoint>& path,std::vector<unsigned>& indexVertex);
33 
34  unsigned best_source(SurfacePoint& point, //quickly find what source this point belongs to and what is the distance to this source
35  double& best_source_distance);
36 
38  {
40 
41  double memory = m_nodes.size()*sizeof(Node);
42  std::cout << "uses about " << memory/1e6 << "Mb of memory" <<std::endl;
43  }
44 
45 protected:
46  unsigned node_index(vertex_pointer v) //gives index of the node that corresponds to this vertex
47  {
48  return v->id();
49  };
50 
51  void set_sources(std::vector<SurfacePoint>& sources)
52  {
53  m_sources = sources;
54  }
55 
56  node_pointer best_first_node(SurfacePoint& point, double& best_total_distance)
57  {
58  node_pointer best_node = NULL;
59  if(point.type() == VERTEX)
60  {
62  best_node = &m_nodes[node_index(v)];
63  best_total_distance = best_node->distance_from_source();
64  }
65  else
66  {
67  std::vector<node_pointer> possible_nodes;
68  list_nodes_visible_from_source(point.base_element(), possible_nodes);
69 
70  best_total_distance = GEODESIC_INF;
71  for(unsigned i=0; i<possible_nodes.size(); ++i)
72  {
73  node_pointer node = possible_nodes[i];
74 
75  double distance_from_dest = node->distance(&point);
76  if(node->distance_from_source() + distance_from_dest < best_total_distance)
77  {
78  best_total_distance = node->distance_from_source() + distance_from_dest;
79  best_node = node;
80  }
81  }
82  }
83 
84  //assert(best_node);
85  //assert(best_total_distance<GEODESIC_INF);
86  if(best_total_distance > m_propagation_distance_stopped) //result is unreliable
87  {
88  best_total_distance = GEODESIC_INF;
89  return NULL;
90  }
91  else
92  {
93  return best_node;
94  }
95  }; //quickly find what node will be the next one in geodesic path
96 
97  bool check_stop_conditions(unsigned& index); //check when propagation should stop
98 
100  std::vector<node_pointer>& storage) = 0; //list all nodes that belong to this mesh element
101 
102  virtual void list_nodes_visible_from_node(node_pointer node, //list all nodes that belong to this mesh element
103  std::vector<node_pointer>& storage,
104  std::vector<double>& distances,
105  double threshold_distance) = 0; //list only the nodes whose current distance is larger than the threshold
106 
107  std::vector<Node> m_nodes; //nodes of the graph
108 
109  typedef std::set<node_pointer, Node> queue_type;
111 
112  std::vector<SurfacePoint> m_sources; //for simplicity, we keep sources as they are
113 };
114 
115 template<class Node>
116 void GeodesicAlgorithmGraphBase<Node>::propagate(std::vector<SurfacePoint>& sources,
117  double max_propagation_distance, //propagation algorithm stops after reaching the certain distance from the source
118  std::vector<SurfacePoint>* stop_points) //or after ensuring that all the stop_points are covered
119 {
120  set_stop_conditions(stop_points, max_propagation_distance);
121  set_sources(sources);
122 
123  m_queue.clear();
124  m_propagation_distance_stopped = GEODESIC_INF;
125  for(unsigned i=0; i<m_nodes.size(); ++i)
126  {
127  m_nodes[i].clear();
128  }
129 
130  clock_t start = clock();
131 
132  std::vector<node_pointer> visible_nodes; //initialize vertices directly visible from sources
133  for(unsigned i=0; i<m_sources.size(); ++i)
134  {
135  SurfacePoint* source = &m_sources[i];
136  list_nodes_visible_from_source(source->base_element(),
137  visible_nodes);
138 
139  for(unsigned j=0; j<visible_nodes.size(); ++j)
140  {
141  node_pointer node = visible_nodes[j];
142  double distance = node->distance(source);
143  if(distance < node->distance_from_source())
144  {
145  node->distance_from_source() = distance;
146  node->source_index() = i;
147  node->previous() = NULL;
148  }
149  }
150  visible_nodes.clear();
151  }
152 
153  for(unsigned i=0; i<m_nodes.size(); ++i) //initialize the queue
154  {
155  if(m_nodes[i].distance_from_source() < GEODESIC_INF)
156  {
157  m_queue.insert(&m_nodes[i]);
158  }
159  }
160 
161  unsigned counter = 0;
162  unsigned satisfied_index = 0;
163 
164  std::vector<double> distances_between_nodes;
165  while(!m_queue.empty()) //main cycle
166  {
167  if(counter++ % 10 == 0) //check if we covered all required vertices
168  {
169  if (check_stop_conditions(satisfied_index))
170  {
171  break;
172  }
173  }
174 
175  node_pointer min_node = *m_queue.begin();
176  m_queue.erase(m_queue.begin());
177  assert(min_node->distance_from_source() < GEODESIC_INF);
178 
179  visible_nodes.clear();
180  distances_between_nodes.clear();
181  list_nodes_visible_from_node(min_node,
182  visible_nodes,
183  distances_between_nodes,
184  min_node->distance_from_source());
185 
186  for(unsigned i=0; i<visible_nodes.size(); ++i) //update all the adgecent vertices
187  {
188  node_pointer next_node = visible_nodes[i];
189 
190  if(next_node->distance_from_source() > min_node->distance_from_source() +
191  distances_between_nodes[i])
192  {
193  if(next_node->distance_from_source() < GEODESIC_INF) //remove it from the queue
194  {
195  typename queue_type::iterator iter = m_queue.find(next_node);
196  assert(iter != m_queue.end());
197  m_queue.erase(iter);
198  }
199  next_node->distance_from_source() = min_node->distance_from_source() +
200  distances_between_nodes[i];
201  next_node->source_index() = min_node->source_index();
202  next_node->previous() = min_node;
203  m_queue.insert(next_node);
204  }
205  }
206  }
207 
208  m_propagation_distance_stopped = m_queue.empty() ? GEODESIC_INF : (*m_queue.begin())->distance_from_source();
209  clock_t finish = clock();
210  m_time_consumed = (static_cast<double>(finish)-static_cast<double>(start))/CLOCKS_PER_SEC;
211  //std::cout << std::endl;
212 }
213 
214 template<class Node>
216 {
217  double queue_min_distance = (*m_queue.begin())->distance_from_source();
218 
219  if(queue_min_distance < m_max_propagation_distance)
220  {
221  return false;
222  }
223 
224  while(index < m_stop_vertices.size())
225  {
226  vertex_pointer v = m_stop_vertices[index].first;
227  Node& node = m_nodes[node_index(v)];
228  if(queue_min_distance < node.distance_from_source() + m_stop_vertices[index].second)
229  {
230  return false;
231  }
232  ++index;
233  }
234  return true;
235 }
236 
237 template<class Node>
238 inline void GeodesicAlgorithmGraphBase<Node>::trace_back(SurfacePoint& destination, //trace back piecewise-linear path
239  std::vector<SurfacePoint>& path)
240 {
241  path.clear();
242 
243  double total_path_length;
244  node_pointer node = best_first_node(destination, total_path_length);
245 
246  if(total_path_length>GEODESIC_INF/2.0) //unable to find the path
247  {
248  return;
249  }
250 
251  path.push_back(destination);
252 
253  if(node->distance(&destination) > 1e-50)
254  {
255  path.push_back(node->surface_point());
256  }
257 
258  while(node->previous()) //follow the path
259  {
260  node = node->previous();
261  path.push_back(node->surface_point());
262  }
263 
264  SurfacePoint& source = m_sources[node->source_index()]; //add source to the path if it is not already there
265  if(node->distance(&source) > 1e-50)
266  {
267  path.push_back(source);
268  }
269 }
270 
271 template<class Node>
272 inline void GeodesicAlgorithmGraphBase<Node>::trace_back_with_index(SurfacePoint& destination, //trace back piecewise-linear path
273  std::vector<SurfacePoint>& path,std::vector<unsigned>& indexVertex)
274 {
275  path.clear();
276  indexVertex.clear();
277 
278  double total_path_length;
279  node_pointer node = best_first_node(destination, total_path_length);
280 
281  if(total_path_length>GEODESIC_INF/2.0) //unable to find the path
282  {
283  return;
284  }
285 
286  path.push_back(destination);
287 
288  if(node->distance(&destination) > 1e-50)
289  {
290  path.push_back(node->surface_point());
291  }
292 
293  while(node->previous()) //follow the path
294  {
295  node = node->previous();
296  SurfacePoint s2 = node->surface_point();
297  SurfacePoint *source2 = &s2;
299  indexVertex.push_back((int)node_index(v));
300 
301  path.push_back(node->surface_point());
302  }
303 
304  SurfacePoint& source = m_sources[node->source_index()]; //add source to the path if it is not already there
305  if(node->distance(&source) > 1e-50)
306  {
307  path.push_back(source);
308  }
309 }
310 template<class Node>
311 inline unsigned GeodesicAlgorithmGraphBase<Node>::best_source(SurfacePoint& point, //quickly find what source this point belongs to and what is the distance to this source
312  double& best_source_distance)
313 {
314  node_pointer node = best_first_node(point, best_source_distance);
315  return node ? node->source_index() : 0;
316 };
317 
318 } //geodesic
319 
320 #endif //GEODESIC_ALGORITHM_GRAPH_BASE_010907
virtual void list_nodes_visible_from_node(node_pointer node, std::vector< node_pointer > &storage, std::vector< double > &distances, double threshold_distance)=0
unsigned best_source(SurfacePoint &point, double &best_source_distance)
Vertex * vertex_pointer
void trace_back(SurfacePoint &destination, std::vector< SurfacePoint > &path)
node_pointer best_first_node(SurfacePoint &point, double &best_total_distance)
void propagate(std::vector< SurfacePoint > &sources, double max_propagation_distance=GEODESIC_INF, std::vector< SurfacePoint > *stop_points=NULL)
void trace_back_with_index(SurfacePoint &destination, std::vector< SurfacePoint > &path, std::vector< unsigned > &indexVertex)
virtual void list_nodes_visible_from_source(MeshElementBase *p, std::vector< node_pointer > &storage)=0
double distance(double *v)
void set_sources(std::vector< SurfacePoint > &sources)