A.I.M.S algorithms


geodesic_algorithm_exact_elements.h
Go to the documentation of this file.
1 //Copyright (C) 2008 Danil Kirsanov, MIT License
2 #ifndef GEODESIC_ALGORITHM_EXACT_ELEMENTS_20071231
3 #define GEODESIC_ALGORITHM_EXACT_ELEMENTS_20071231
4 
5 #include "geodesic_memory.h"
7 #include <vector>
8 #include <cmath>
9 #include <assert.h>
10 #include <algorithm>
11 
12 namespace geodesic{
13 
14 class Interval;
16 typedef Interval* interval_pointer;
18 
19 class Interval //interval of the edge
20 {
21 public:
22 
23  Interval(){};
24  ~Interval(){};
25 
27  {
32  };
33 
34  double signal(double x) //geodesic distance function at point x
35  {
36  assert(x>=0.0 && x <= m_edge->length());
37 
38  if(m_d == GEODESIC_INF)
39  {
40  return GEODESIC_INF;
41  }
42  else
43  {
44  double dx = x - m_pseudo_x;
45  if(m_pseudo_y == 0.0)
46  {
47  return m_d + std::abs(dx);
48  }
49  else
50  {
51  return m_d + sqrt(dx*dx + m_pseudo_y*m_pseudo_y);
52  }
53  }
54  }
55 
56  double max_distance(double end)
57  {
58  if(m_d == GEODESIC_INF)
59  {
60  return GEODESIC_INF;
61  }
62  else
63  {
64  double a = std::abs(m_start - m_pseudo_x);
65  double b = std::abs(end - m_pseudo_x);
66 
67  return a > b ? m_d + sqrt(a*a + m_pseudo_y*m_pseudo_y):
68  m_d + sqrt(b*b + m_pseudo_y*m_pseudo_y);
69  }
70  }
71 
72  void compute_min_distance(double stop) //compute min, given c,d theta, start, end.
73  {
74  assert(stop > m_start);
75 
76  if(m_d == GEODESIC_INF)
77  {
79  }
80  else if(m_start > m_pseudo_x)
81  {
82  m_min = signal(m_start);
83  }
84  else if(stop < m_pseudo_x)
85  {
86  m_min = signal(stop);
87  }
88  else
89  {
90  assert(m_pseudo_y<=0);
91  m_min = m_d - m_pseudo_y;
92  }
93  }
94  //compare two intervals in the queue
95  bool operator()(interval_pointer const x, interval_pointer const y) const
96  {
97  if(x->min() != y->min())
98  {
99  return x->min() < y->min();
100  }
101  else if(x->start() != y->start())
102  {
103  return x->start() < y->start();
104  }
105  else
106  {
107  return x->edge()->id() < y->edge()->id();
108  }
109  }
110 
111  double stop() //return the endpoint of the interval
112  {
113  return m_next ? m_next->start() : m_edge->length();
114  }
115 
116  double hypotenuse(double a, double b)
117  {
118  return sqrt(a*a + b*b);
119  }
120 
121  void find_closest_point(double const x,
122  double const y,
123  double& offset,
124  double& distance); //find the point on the interval that is closest to the point (alpha, s)
125 
126  double& start(){return m_start;};
127  double& d(){return m_d;};
128  double& pseudo_x(){return m_pseudo_x;};
129  double& pseudo_y(){return m_pseudo_y;};
130  double& min(){return m_min;};
132  edge_pointer& edge(){return m_edge;};
135  unsigned& source_index(){return m_source_index;};
136 
138  SurfacePoint* point = NULL,
139  unsigned source_index = 0);
140 
141 protected:
142  double m_start; //initial point of the interval on the edge
143  double m_d; //distance from the source to the pseudo-source
144  double m_pseudo_x; //coordinates of the pseudo-source in the local coordinate system
145  double m_pseudo_y; //y-coordinate should be always negative
146  double m_min; //minimum distance on the interval
147 
148  interval_pointer m_next; //pointer to the next interval in the list
149  edge_pointer m_edge; //edge that the interval belongs to
150  unsigned m_source_index; //the source it belongs to
151  DirectionType m_direction; //where the interval is coming from
152 };
153 
154 struct IntervalWithStop : public Interval
155 {
156 public:
157  double& stop(){return m_stop;};
158 protected:
159  double m_stop;
160 };
161 
162 class IntervalList //list of the of intervals of the given edge
163 {
164 public:
165  IntervalList(){m_first = NULL;};
167 
168  void clear()
169  {
170  m_first = NULL;
171  };
172 
174  {
175  m_edge = e;
176  m_first = NULL;
177  };
178 
179  interval_pointer covering_interval(double offset) //returns the interval that covers the offset
180  {
181  assert(offset >= 0.0 && offset <= m_edge->length());
182 
183  interval_pointer p = m_first;
184  while(p && p->stop() < offset)
185  {
186  p = p->next();
187  }
188 
189  return p;// && p->start() <= offset ? p : NULL;
190  };
191 
193  double& offset,
194  double& distance,
195  interval_pointer& interval)
196  {
197  interval_pointer p = m_first;
198  distance = GEODESIC_INF;
199  interval = NULL;
200 
201  double x,y;
202  m_edge->local_coordinates(point, x, y);
203 
204  while(p)
205  {
206  if(p->min()<GEODESIC_INF)
207  {
208  double o, d;
209  p->find_closest_point(x, y, o, d);
210  if(d < distance)
211  {
212  distance = d;
213  offset = o;
214  interval = p;
215  }
216  }
217  p = p->next();
218  }
219  };
220 
222  {
223  interval_pointer p = m_first;
224  unsigned count = 0;
225  while(p)
226  {
227  ++count;
228  p = p->next();
229  }
230  return count;
231  }
232 
234  {
235  interval_pointer p = m_first;
236  if(p)
237  {
238  while(p->next())
239  {
240  p = p->next();
241  }
242  }
243  return p;
244  }
245 
246  double signal(double x)
247  {
248  interval_pointer interval = covering_interval(x);
249 
250  return interval ? interval->signal(x) : GEODESIC_INF;
251  }
252 
253  interval_pointer& first(){return m_first;};
254  edge_pointer& edge(){return m_edge;};
255 private:
256  interval_pointer m_first; //pointer to the first member of the list
257  edge_pointer m_edge; //edge that owns this list
258 };
259 
261 {
262 public:
263  unsigned index(){return m_index;};
264 
265  void initialize(SurfacePoint& p, unsigned index)
266  {
268  m_index = index;
269  }
270 
271  bool operator()(SurfacePointWithIndex* x, SurfacePointWithIndex* y) const //used for sorting
272  {
273  assert(x->type() != UNDEFINED_POINT && y->type() !=UNDEFINED_POINT);
274 
275  if(x->type() != y->type())
276  {
277  return x->type() < y->type();
278  }
279  else
280  {
281  return x->base_element()->id() < y->base_element()->id();
282  }
283  }
284 
285 private:
286  unsigned m_index;
287 };
288 
289 class SortedSources : public std::vector<SurfacePointWithIndex>
290 {
291 private:
292  typedef std::vector<SurfacePointWithIndex*> sorted_vector_type;
293 public:
294  typedef sorted_vector_type::iterator sorted_iterator;
295  typedef std::pair<sorted_iterator, sorted_iterator> sorted_iterator_pair;
296 
298  {
299  m_search_dummy.base_element() = mesh_element;
300 
301  return equal_range(m_sorted.begin(),
302  m_sorted.end(),
303  &m_search_dummy,
304  m_compare_less);
305  }
306 
307  void initialize(std::vector<SurfacePoint>& sources) //we initialize the sources by copie
308  {
309  resize(sources.size());
310  m_sorted.resize(sources.size());
311  for(unsigned i=0; i<sources.size(); ++i)
312  {
313  SurfacePointWithIndex& p = *(begin() + i);
314 
315  p.initialize(sources[i],i);
316  m_sorted[i] = &p;
317  }
318 
319  std::sort(m_sorted.begin(), m_sorted.end(), m_compare_less);
320  };
321 
323  {
324  assert(i < size());
325  return *(begin() + i);
326  }
327 
328 private:
329  sorted_vector_type m_sorted;
330  SurfacePointWithIndex m_search_dummy; //used as a search template
331  SurfacePointWithIndex m_compare_less; //used as a compare functor
332 };
333 
334 
335 inline void Interval::find_closest_point(double const rs,
336  double const hs,
337  double& r,
338  double& d_out) //find the point on the interval that is closest to the point (alpha, s)
339  {
340  if(m_d == GEODESIC_INF)
341  {
342  r = GEODESIC_INF;
343  d_out = GEODESIC_INF;
344  return;
345  }
346 
347  double hc = -m_pseudo_y;
348  double rc = m_pseudo_x;
349  double end = stop();
350 
351  double local_epsilon = SMALLEST_INTERVAL_RATIO*m_edge->length();
352  if(std::abs(hs+hc) < local_epsilon)
353  {
354  if(rs<=m_start)
355  {
356  r = m_start;
357  d_out = signal(m_start) + std::abs(rs - m_start);
358  }
359  else if(rs>=end)
360  {
361  r = end;
362  d_out = signal(end) + fabs(end - rs);
363  }
364  else
365  {
366  r = rs;
367  d_out = signal(rs);
368  }
369  }
370  else
371  {
372  double ri = (rs*hc + hs*rc)/(hs+hc);
373 
374  if(ri<m_start)
375  {
376  r = m_start;
377  d_out = signal(m_start) + hypotenuse(m_start - rs, hs);
378  }
379  else if(ri>end)
380  {
381  r = end;
382  d_out = signal(end) + hypotenuse(end - rs, hs);
383  }
384  else
385  {
386  r = ri;
387  d_out = m_d + hypotenuse(rc - rs, hc + hs);
388  }
389  }
390  }
391 
392 
394  SurfacePoint* source,
395  unsigned source_index)
396 {
397  m_next = NULL;
398  //m_geodesic_previous = NULL;
400  m_edge = edge;
402 
403  m_start = 0.0;
404  //m_stop = edge->length();
405  if(!source)
406  {
407  m_d = GEODESIC_INF;
409  return;
410  }
411  m_d = 0;
412 
413  if(source->base_element()->type() == VERTEX)
414  {
415  if(source->base_element()->id() == edge->v0()->id())
416  {
417  m_pseudo_x = 0.0;
418  m_pseudo_y = 0.0;
419  m_min = 0.0;
420  return;
421  }
422  else if(source->base_element()->id() == edge->v1()->id())
423  {
424  m_pseudo_x = stop();
425  m_pseudo_y = 0.0;
426  m_min = 0.0;
427  return;
428  }
429  }
430 
431  edge->local_coordinates(source, m_pseudo_x, m_pseudo_y);
433 
435 }
436 
437 } //geodesic
438 
439 #endif //GEODESIC_ALGORITHM_EXACT_ELEMENTS_20071231
bool operator()(interval_pointer const x, interval_pointer const y) const
bool operator()(SurfacePointWithIndex *x, SurfacePointWithIndex *y) const
void local_coordinates(Point3D *point, double &x, double &y)
std::pair< sorted_iterator, sorted_iterator > sorted_iterator_pair
void initialize(edge_pointer edge, SurfacePoint *point=NULL, unsigned source_index=0)
vertex_pointer v0()
void initialize(SurfacePoint const &p)
sorted_vector_type::iterator sorted_iterator
interval_pointer covering_interval(double offset)
void find_closest_point(double const x, double const y, double &offset, double &distance)
void initialize(SurfacePoint &p, unsigned index)
vertex_pointer v1()
double hypotenuse(double a, double b)
void initialize(std::vector< SurfacePoint > &sources)
sorted_iterator_pair sources(base_pointer mesh_element)
SurfacePointWithIndex & operator[](unsigned i)
void find_closest_point(SurfacePoint *point, double &offset, double &distance, interval_pointer &interval)