A.I.M.S algorithms


geodesic_algorithm_dijkstra_alternative.h
Go to the documentation of this file.
1 //Copyright (C) 2008 Danil Kirsanov, MIT License
2 #ifndef GEODESIC_ALGORITHM_DIJKSTRA_ALTERNATIVE_010506
3 #define GEODESIC_ALGORITHM_DIJKSTRA_ALTERNATIVE_010506
4 
7 #include <vector>
8 #include <set>
9 #include <assert.h>
10 
11 namespace geodesic{
12 
14 {
15  typedef DijkstraNode1* node_pointer;
16 public:
19 
20  double& distance_from_source(){return m_distance;};
21  node_pointer& previous(){return m_previous;};
22  unsigned& source_index(){return m_source_index;};
23  vertex_pointer& vertex(){return m_vertex;};
24 
25  void clear()
26  {
27  m_distance = GEODESIC_INF;
28  m_previous = NULL;
29  }
30 
31  bool operator()(node_pointer const s1, node_pointer const s2) const
32  {
33  return s1->distance_from_source() != s2->distance_from_source() ?
35  s1->vertex()->id() < s2->vertex()->id();
36  };
37 
38 private:
39  double m_distance; //distance to the closest source
40  unsigned m_source_index; //closest source index
41  node_pointer m_previous; //previous node in the geodesic path
42  vertex_pointer m_vertex; //correspoding vertex
43 };
44 
46 {
47 public:
49  typedef Node* node_pointer;
50 
53  m_nodes(mesh->vertices().size())
54  {
55  m_type = DIJKSTRA;
56  for(unsigned i=0; i<m_nodes.size(); ++i)
57  {
58  m_nodes[i].vertex() = &m_mesh->vertices()[i];
59  }
60  };
61 
63 
64  virtual void propagate(std::vector<SurfacePoint>& sources,
65  double max_propagation_distance = GEODESIC_INF, //propagation algorithm stops after reaching the certain distance from the source
66  std::vector<SurfacePoint>* stop_points = NULL); //or after ensuring that all the stop_points are covered
67 
68  virtual void trace_back(SurfacePoint& destination, //trace back piecewise-linear path
69  std::vector<SurfacePoint>& path);
70 
71  virtual unsigned best_source(SurfacePoint& point, //quickly find what source this point belongs to and what is the distance to this source
72  double& best_source_distance);
73 private:
74 
75  void set_sources(std::vector<SurfacePoint>& sources)
76  {
77  m_sources = sources;
78  }
79 
80  bool check_stop_conditions(unsigned& index);
81 
82  std::vector<Node> m_nodes; //nodes of the subdivision graph located on the vertices
83 
84  std::vector<SurfacePoint> m_sources; //for simplicity, we keep sources as they are
85 
86  typedef std::set<node_pointer, Node> queue_type;
87  queue_type m_queue;
88 };
89 
90 inline unsigned GeodesicAlgorithmDijkstraAlternative::best_source(SurfacePoint& point, //quickly find what source this point belongs to and what is the distance to this source
91  double& best_source_distance)
92 {
93  if(point.type() == VERTEX)
94  {
96  node_pointer node = &m_nodes[v->id()];
97  best_source_distance = node->distance_from_source();
98  return node->source_index();
99  }
100  else
101  {
102  std::vector<vertex_pointer> possible_vertices; //initialize vertices directly visible from sources
103  m_mesh->closest_vertices(&point, &possible_vertices);
104 
105  best_source_distance = GEODESIC_INF;
106  vertex_pointer min_vertex = NULL;
107  for(unsigned i=0; i<possible_vertices.size(); ++i)
108  {
109  vertex_pointer v = possible_vertices[i];
110 
111  double distance_from_source = m_nodes[v->id()].distance_from_source();
112  double distance_from_dest = v->distance(&point);
113  if(distance_from_source + distance_from_dest < best_source_distance)
114  {
115  best_source_distance = distance_from_source + distance_from_dest;
116  min_vertex = v;
117  }
118  }
119  assert(min_vertex);
120  node_pointer node = &m_nodes[min_vertex->id()];
121  return node->source_index();
122  }
123 }
124 
125 inline void GeodesicAlgorithmDijkstraAlternative::trace_back(SurfacePoint& destination, //trace back piecewise-linear path
126  std::vector<SurfacePoint>& path)
127 {
128  path.clear();
129  path.push_back(destination);
130 
131  if(destination.type() != VERTEX) //snap to the closest vertex
132  {
133  std::vector<vertex_pointer> possible_vertices; //initialize vertices directly visible from sources
134  m_mesh->closest_vertices(&destination, &possible_vertices);
135 
136  double min_distance = GEODESIC_INF;
137  vertex_pointer min_vertex = NULL;
138  for(unsigned i=0; i<possible_vertices.size(); ++i)
139  {
140  vertex_pointer v = possible_vertices[i];
141 
142  double distance_from_source = m_nodes[v->id()].distance_from_source();
143  double distance_from_dest = v->distance(&destination);
144  if(distance_from_source + distance_from_dest < min_distance)
145  {
146  min_distance = distance_from_source + distance_from_dest;
147  min_vertex = v;
148  }
149  }
150  assert(min_vertex);
151  path.push_back(SurfacePoint(min_vertex));
152  }
153 
154  node_pointer node = &m_nodes[path.back().base_element()->id()];
155  while(node->previous()) //follow the path
156  {
157  node = node->previous();
158  path.push_back(SurfacePoint(node->vertex()));
159  }
160 
161  SurfacePoint& source = m_sources[node->source_index()]; //add source to the path if it is not already there
162  if(source.type() != VERTEX ||
163  source.base_element()->id() != node->vertex()->id())
164  {
165  path.push_back(source);
166  }
167 }
168 
169 inline void GeodesicAlgorithmDijkstraAlternative::propagate(std::vector<SurfacePoint>& sources,
170  double max_propagation_distance, //propagation algorithm stops after reaching the certain distance from the source
171  std::vector<SurfacePoint>* stop_points) //or after ensuring that all the stop_points are covered
172 {
173  set_stop_conditions(stop_points, max_propagation_distance);
174  set_sources(sources);
175 
176  m_queue.clear(); //clear everything
177  for(unsigned i=0; i<m_nodes.size(); ++i)
178  {
179  m_nodes[i].clear();
180  }
181 
182  clock_t start = clock();
183 
184  std::vector<vertex_pointer> visible_vertices; //initialize vertices directly visible from sources
185  for(unsigned i=0; i<m_sources.size(); ++i)
186  {
187  SurfacePoint* s = &m_sources[i];
188  m_mesh->closest_vertices(s, &visible_vertices);
189  for(unsigned j=0; j<visible_vertices.size(); ++j)
190  {
191  vertex_pointer v = visible_vertices[j];
192  double distance = s->distance(v);
193 
194  node_pointer node = &m_nodes[v->id()];
195  if(distance < node->distance_from_source())
196  {
197  node->distance_from_source() = distance;
198  node->source_index() = i;
199  node->previous() = NULL;
200  }
201  }
202  visible_vertices.clear();
203  }
204 
205  for(unsigned i=0; i<m_nodes.size(); ++i) //initialize the queue
206  {
207  if(m_nodes[i].distance_from_source() < GEODESIC_INF)
208  {
209  m_queue.insert(&m_nodes[i]);
210  }
211  }
212 
213  unsigned counter = 0;
214  unsigned satisfied_index = 0;
215 
216  while(!m_queue.empty()) //main cycle
217  {
218  if(counter++ % 10 == 0) //check if we covered all required vertices
219  {
220  if (check_stop_conditions(satisfied_index))
221  {
222  break;
223  }
224  //std::cout << counter << " " << m_queue.size() << " " << (*m_queue.begin())->distance_from_source()<< std::endl;
225  }
226 
227  node_pointer min_node = *m_queue.begin();
228  m_queue.erase(m_queue.begin());
229 
230  vertex_pointer v = min_node->vertex();
231 
232  for(unsigned i=0; i<v->adjacent_edges().size(); ++i) //update all the adgecent vertices
233  {
234  edge_pointer e = v->adjacent_edges()[i];
235  vertex_pointer next_v = e->opposite_vertex(v);
236  node_pointer next_node = &m_nodes[next_v->id()];
237  //double current_distance =
238  if(next_node->distance_from_source() > min_node->distance_from_source() + e->length())
239  {
240  if(next_node->distance_from_source() < GEODESIC_INF) //remove it from the queue
241  {
242  queue_type::iterator iter = m_queue.find(next_node);
243  assert(iter != m_queue.end());
244  m_queue.erase(iter);
245  }
246  next_node->distance_from_source() = min_node->distance_from_source() + e->length();
247  next_node->source_index() = min_node->source_index();
248  next_node->previous() = min_node;
249  m_queue.insert(next_node);
250  }
251  }
252  }
253  //std::cout << std::endl;
254 
255  clock_t finish = clock();
256  m_time_consumed = (static_cast<double>(finish)-static_cast<double>(start))/CLOCKS_PER_SEC;
257 }
258 
259 inline bool GeodesicAlgorithmDijkstraAlternative::check_stop_conditions(unsigned& index)
260 {
261  double queue_min_distance = (*m_queue.begin())->distance_from_source();
262 
263  if(queue_min_distance < m_max_propagation_distance)
264  {
265  return false;
266  }
267 
268  while(index < m_stop_vertices.size())
269  {
270  vertex_pointer v = m_stop_vertices[index].first;
271  Node& node = m_nodes[v->id()];
272  if(queue_min_distance < node.distance_from_source() + m_stop_vertices[index].second)
273  {
274  return false;
275  }
276  ++index;
277  }
278  return true;
279 }
280 
281 } //geodesic
282 
283 #endif //GEODESIC_ALGORITHM_DIJKSTRA_ALTERNATIVE_010506
virtual unsigned best_source(SurfacePoint &point, double &best_source_distance)
edge_pointer_vector & adjacent_edges()
void set_stop_conditions(std::vector< SurfacePoint > *stop_points, double stop_distance)
std::vector< stop_vertex_with_distace_type > m_stop_vertices
bool operator()(node_pointer const s1, node_pointer const s2) const
Vertex * vertex_pointer
std::vector< Vertex > & vertices()
Definition: geodesic_mesh.h:43
virtual void trace_back(SurfacePoint &destination, std::vector< SurfacePoint > &path)
double distance(double *v)
vertex_pointer opposite_vertex(vertex_pointer v)
virtual void propagate(std::vector< SurfacePoint > &sources, double max_propagation_distance=GEODESIC_INF, std::vector< SurfacePoint > *stop_points=NULL)
unsigned closest_vertices(SurfacePoint *p, std::vector< vertex_pointer > *storage=NULL)
Definition: geodesic_mesh.h:68