aimsalgo  5.1.2
Neuroimaging image processing
geodesic_algorithm_base.h
Go to the documentation of this file.
1 //Copyright (C) 2008 Danil Kirsanov, MIT License
2 
3 #ifndef GEODESIC_ALGORITHM_BASE_122806
4 #define GEODESIC_ALGORITHM_BASE_122806
5 
6 #include "geodesic_mesh.h"
8 #include <iostream>
9 #include <ctime>
10 
11 namespace geodesic{
12 
14 {
15 public:
17  {
22  };
23 
27  m_mesh(mesh)
28  {};
29 
31 
32  virtual void propagate(std::vector<SurfacePoint>& sources,
33  double max_propagation_distance = GEODESIC_INF, //propagation algorithm stops after reaching the certain distance from the source
34  std::vector<SurfacePoint>* stop_points = NULL) = 0; //or after ensuring that all the stop_points are covered
35 
36  virtual void trace_back(SurfacePoint& destination, //trace back piecewise-linear path
37  std::vector<SurfacePoint>& path) = 0;
38 
39  virtual void trace_back_with_index(SurfacePoint& destination,
40  std::vector<SurfacePoint>& path,std::vector<unsigned>& indexVertex) = 0;
41 
42  void geodesic(SurfacePoint& source,
43  SurfacePoint& destination,
44  std::vector<SurfacePoint>& path); //lazy people can find geodesic path with one function call
45 
46  void geodesic(std::vector<SurfacePoint>& sources,
47  std::vector<SurfacePoint>& destinations,
48  std::vector<std::vector<SurfacePoint> >& paths); //lazy people can find geodesic paths with one function call
49 
50  void geodesic(SurfacePoint& source,
51  SurfacePoint& destination,
52  std::vector<SurfacePoint>& path,std::vector<unsigned>& indexVertex);
53 
54  void geodesic(SurfacePoint& sources,
55  std::vector<SurfacePoint>& destinations,
56  std::vector<std::vector<SurfacePoint> >& paths,
57  std::vector<std::vector<unsigned> >& indexVertex);
58 
59  double length(std::vector<SurfacePoint>& path);
60  void print_info_about_path(std::vector<SurfacePoint>& path);
61 
62  virtual unsigned best_source(SurfacePoint& point, //after propagation step is done, quickly find what source this point belongs to and what is the distance to this source
63  double& best_source_distance) = 0;
64 
65  virtual void print_statistics() //print info about timing and memory usage in the propagation step of the algorithm
66  {
67  std::cout << "propagation step took " << m_time_consumed << " seconds " << std::endl;
68  };
69 
71 
72  virtual std::string name();
73 
74  geodesic::Mesh* mesh(){return m_mesh;};
75 protected:
76 
77  void set_stop_conditions(std::vector<SurfacePoint>* stop_points,
78  double stop_distance);
79  double stop_distance()
80  {
82  }
83 
84  AlgorithmType m_type; // type of the algorithm
85 
86  typedef std::pair<vertex_pointer, double> stop_vertex_with_distace_type;
87  std::vector<stop_vertex_with_distace_type> m_stop_vertices; // algorithm stops propagation after covering certain vertices
88  double m_max_propagation_distance; // or reaching the certain distance
89 
91 
92  double m_time_consumed; //how much time does the propagation step takes
93  double m_propagation_distance_stopped; //at what distance (if any) the propagation algorithm stopped
94 };
95 
96 inline double GeodesicAlgorithmBase::length(std::vector<SurfacePoint>& path)
97 {
98  double length = 0;
99  if(!path.empty())
100  {
101  for(unsigned i=0; i<path.size()-1; ++i)
102  {
103  length += path[i].distance(&path[i+1]);
104  }
105  }
106  return length;
107 }
108 
109 inline void GeodesicAlgorithmBase::print_info_about_path(std::vector<SurfacePoint>& path)
110 {
111  std::cout << "number of the points in the path = " << path.size()
112  << ", length of the path = " << length(path)
113  << std::endl;
114 }
115 
116 inline std::string GeodesicAlgorithmBase::name()
117 {
118  switch(m_type)
119  {
120  case EXACT:
121  return "exact";
122  case DIJKSTRA:
123  return "dijkstra";
124  case SUBDIVISION:
125  return "subdivision";
126  default:
127  case UNDEFINED_ALGORITHM:
128  return "undefined";
129  }
130 }
131 
133  SurfacePoint& destination,
134  std::vector<SurfacePoint>& path) //lazy people can find geodesic path with one function call
135 {
136  std::vector<SurfacePoint> sources(1, source);
137  std::vector<SurfacePoint> stop_points(1, destination);
138  double const max_propagation_distance = GEODESIC_INF;
139 
140  propagate(sources,
141  max_propagation_distance,
142  &stop_points);
143 
144  trace_back(destination, path);
145 }
146 
148  SurfacePoint& destination,
149  std::vector<SurfacePoint>& path,std::vector<unsigned>& indexVertex) //lazy people can find geodesic path with one function call
150 {
151  std::vector<SurfacePoint> sources(1, source);
152  std::vector<SurfacePoint> stop_points(1, destination);
153  double const max_propagation_distance = GEODESIC_INF;
154 
155  propagate(sources,
156  max_propagation_distance,
157  &stop_points);
158 
159  trace_back_with_index(destination, path,indexVertex);
160 }
161 
162 inline void GeodesicAlgorithmBase::geodesic(std::vector<SurfacePoint>& sources,
163  std::vector<SurfacePoint>& destinations,
164  std::vector<std::vector<SurfacePoint> >& paths) //lazy people can find geodesic paths with one function call
165 {
166  double const max_propagation_distance = GEODESIC_INF;
167 
168  propagate(sources,
169  max_propagation_distance,
170  &destinations); //we use desinations as stop points
171 
172  paths.resize(destinations.size());
173 
174  for(unsigned i=0; i<paths.size(); ++i)
175  {
176  trace_back(destinations[i], paths[i]);
177  }
178 }
179 
181  std::vector<SurfacePoint>& destinations, std::vector<std::vector<SurfacePoint> >& paths,
182  std::vector<std::vector<unsigned> >& indexVertex) //lazy people can find geodesic paths with one function call
183 {
184  double const max_propagation_distance = GEODESIC_INF;
185 
186  std::vector<SurfacePoint> sources(1, source);
187 
188  propagate(sources,
189  max_propagation_distance,
190  &destinations); //we use desinations as stop points
191 
192  paths.resize(destinations.size());
193  indexVertex.resize(destinations.size());
194 
195  for(unsigned i=0; i<paths.size(); ++i)
196  {
197  trace_back_with_index(destinations[i], paths[i],indexVertex[i]);
198  }
199 }
200 
201 inline void GeodesicAlgorithmBase::set_stop_conditions(std::vector<SurfacePoint>* stop_points,
202  double stop_distance)
203 {
205 
206  if(!stop_points)
207  {
208  m_stop_vertices.clear();
209  return;
210  }
211 
212  m_stop_vertices.resize(stop_points->size());
213 
214  std::vector<vertex_pointer> possible_vertices;
215  for(unsigned i = 0; i < stop_points->size(); ++i)
216  {
217  SurfacePoint* point = &(*stop_points)[i];
218 
219  possible_vertices.clear();
220  m_mesh->closest_vertices(point, &possible_vertices);
221 
222  vertex_pointer closest_vertex = NULL;
223  double min_distance = 1e100;
224  for(unsigned j = 0; j < possible_vertices.size(); ++j)
225  {
226  double distance = point->distance(possible_vertices[j]);
227  if(distance < min_distance)
228  {
229  min_distance = distance;
230  closest_vertex = possible_vertices[j];
231  }
232  }
233  assert(closest_vertex);
234 
235  m_stop_vertices[i].first = closest_vertex;
236  m_stop_vertices[i].second = min_distance;
237  }
238 }
239 
240 }//geodesic
241 
242 #endif //GEODESIC_ALGORITHM_BASE_122806
GeodesicAlgorithmBase(geodesic::Mesh *mesh)
virtual void trace_back_with_index(SurfacePoint &destination, std::vector< SurfacePoint > &path, std::vector< unsigned > &indexVertex)=0
virtual void trace_back(SurfacePoint &destination, std::vector< SurfacePoint > &path)=0
virtual unsigned best_source(SurfacePoint &point, double &best_source_distance)=0
void set_stop_conditions(std::vector< SurfacePoint > *stop_points, double stop_distance)
void print_info_about_path(std::vector< SurfacePoint > &path)
double length(std::vector< SurfacePoint > &path)
void geodesic(SurfacePoint &source, SurfacePoint &destination, std::vector< SurfacePoint > &path)
std::vector< stop_vertex_with_distace_type > m_stop_vertices
virtual void propagate(std::vector< SurfacePoint > &sources, double max_propagation_distance=GEODESIC_INF, std::vector< SurfacePoint > *stop_points=NULL)=0
std::pair< vertex_pointer, double > stop_vertex_with_distace_type
unsigned closest_vertices(SurfacePoint *p, std::vector< vertex_pointer > *storage=NULL)
Definition: geodesic_mesh.h:68
double distance(double *v)