geodesic_algorithm_exact_elements.h

Go to the documentation of this file.

//Copyright (C) 2008 Danil Kirsanov, MIT License
#ifndef GEODESIC_ALGORITHM_EXACT_ELEMENTS_20071231
#define GEODESIC_ALGORITHM_EXACT_ELEMENTS_20071231

#include "geodesic_memory.h"
#include "geodesic_mesh_elements.h"
#include <vector>
#include <cmath>
#include <assert.h>
#include <algorithm>

namespace geodesic{

class Interval;
class IntervalList;
typedef Interval* interval_pointer;
typedef IntervalList* list_pointer;

class Interval            //interval of the edge
{
public:

  Interval(){};
  ~Interval(){};

  enum DirectionType
    {
        FROM_FACE_0,
    FROM_FACE_1,
    FROM_SOURCE,
    UNDEFINED_DIRECTION
    };

  double signal(double x)   //geodesic distance function at point x
  {
    assert(x>=0.0 && x <= m_edge->length());

    if(m_d == GEODESIC_INF)
    {
      return GEODESIC_INF;
    }
    else
    {
      double dx = x - m_pseudo_x;
      if(m_pseudo_y == 0.0)
      {
        return m_d + std::abs(dx);
      }
      else
      {
        return m_d + sqrt(dx*dx + m_pseudo_y*m_pseudo_y);
      }
    }
  }

  double max_distance(double end)
  {
    if(m_d == GEODESIC_INF)
    {
      return GEODESIC_INF;
    }
    else
    {
      double a = std::abs(m_start - m_pseudo_x);
      double b = std::abs(end - m_pseudo_x);

      return a > b ? m_d + sqrt(a*a + m_pseudo_y*m_pseudo_y):
               m_d + sqrt(b*b + m_pseudo_y*m_pseudo_y);
    }
  }

  void compute_min_distance(double stop)      //compute min, given c,d theta, start, end.
  {
    assert(stop > m_start);

    if(m_d == GEODESIC_INF)
    {
      m_min = GEODESIC_INF;
    }
    else if(m_start > m_pseudo_x)
    {
      m_min = signal(m_start);
    }
    else if(stop < m_pseudo_x)
    {
      m_min = signal(stop);
    }
    else
    {
      assert(m_pseudo_y<=0);
      m_min = m_d - m_pseudo_y;
    }
  }
      //compare two intervals in the queue
  bool operator()(interval_pointer const x, interval_pointer const y) const
  {
    if(x->min() != y->min())
    {
      return x->min() < y->min();
    }
    else if(x->start() != y->start())
    {
      return x->start() < y->start();
    }
    else
    {
      return x->edge()->id() < y->edge()->id();
    }
  }

  double stop()   //return the endpoint of the interval
  {
    return m_next ? m_next->start() : m_edge->length();
  }

  double hypotenuse(double a, double b)
  {
    return sqrt(a*a + b*b);
  }

  void find_closest_point(double const x,
                double const y,
                double& offset,
                double& distance);      //find the point on the interval that is closest to the point (alpha, s)

  double& start(){return m_start;};
  double& d(){return m_d;};
  double& pseudo_x(){return m_pseudo_x;};
  double& pseudo_y(){return m_pseudo_y;};
  double& min(){return m_min;};
  interval_pointer& next(){return m_next;};
  edge_pointer& edge(){return m_edge;};
  DirectionType& direction(){return m_direction;};
  bool visible_from_source(){return m_direction == FROM_SOURCE;};
  unsigned& source_index(){return m_source_index;};

  void initialize(edge_pointer edge,
          SurfacePoint* point = NULL,
          unsigned source_index = 0);

protected:
  double m_start;           //initial point of the interval on the edge
  double m_d;             //distance from the source to the pseudo-source
  double m_pseudo_x;          //coordinates of the pseudo-source in the local coordinate system
  double m_pseudo_y;          //y-coordinate should be always negative
  double m_min;           //minimum distance on the interval

  interval_pointer m_next;      //pointer to the next interval in the list
  edge_pointer m_edge;        //edge that the interval belongs to
  unsigned m_source_index;      //the source it belongs to
  DirectionType m_direction;      //where the interval is coming from
};

struct IntervalWithStop : public Interval
{
public:
  double& stop(){return m_stop;};
protected:
  double m_stop;
};

class IntervalList            //list of the of intervals of the given edge
{
public:
  IntervalList(){m_first = NULL;};
  ~IntervalList(){};

  void clear()
  {
    m_first = NULL;
  };

  void initialize(edge_pointer e)
  {
    m_edge = e;
    m_first = NULL;
  };

  interval_pointer covering_interval(double offset)     //returns the interval that covers the offset
  {
    assert(offset >= 0.0 && offset <= m_edge->length());

    interval_pointer p = m_first;
    while(p && p->stop() < offset)
    {
      p = p->next();
    }

    return p;// && p->start() <= offset ? p : NULL;
  };

  void find_closest_point(SurfacePoint* point,
              double& offset,
              double& distance,
              interval_pointer& interval)
  {
    interval_pointer p = m_first;
    distance = GEODESIC_INF;
    interval = NULL;

    double x,y;
    m_edge->local_coordinates(point, x, y);

    while(p)
    {
      if(p->min()<GEODESIC_INF)
      {
        double o, d;
        p->find_closest_point(x, y, o, d);
        if(d < distance)
        {
          distance = d;
          offset = o;
          interval = p;
        }
      }
      p = p->next();
    }
  };

  unsigned number_of_intervals()
  {
    interval_pointer p = m_first;
    unsigned count = 0;
    while(p)
    {
      ++count;
      p = p->next();
    }
    return count;
  }

  interval_pointer last()
  {
    interval_pointer p = m_first;
    if(p)
    {
      while(p->next())
      {
        p = p->next();
      }
    }
    return p;
  }

  double signal(double x)
  {
    interval_pointer interval = covering_interval(x);

    return interval ? interval->signal(x) : GEODESIC_INF;
  }

  interval_pointer& first(){return m_first;};
  edge_pointer& edge(){return m_edge;};
private:
  interval_pointer m_first;     //pointer to the first member of the list
  edge_pointer m_edge;        //edge that owns this list
};

class SurfacePointWithIndex : public SurfacePoint
{
public:
  unsigned index(){return m_index;};

  void initialize(SurfacePoint& p, unsigned index)
  {
    SurfacePoint::initialize(p);
    m_index = index;
  }

  bool operator()(SurfacePointWithIndex* x, SurfacePointWithIndex* y) const //used for sorting
  {
    assert(x->type() != UNDEFINED_POINT && y->type() !=UNDEFINED_POINT);

    if(x->type() != y->type())
    {
      return x->type() < y->type();
    }
    else
    {
      return x->base_element()->id() < y->base_element()->id();
    }
  }

private:
  unsigned m_index;
}; 

class SortedSources : public std::vector<SurfacePointWithIndex>
{
private: 
  typedef std::vector<SurfacePointWithIndex*> sorted_vector_type;
public:
  typedef sorted_vector_type::iterator sorted_iterator;
  typedef std::pair<sorted_iterator, sorted_iterator> sorted_iterator_pair;

  sorted_iterator_pair sources(base_pointer mesh_element)
  {
    m_search_dummy.base_element() = mesh_element;

    return equal_range(m_sorted.begin(),
               m_sorted.end(),
               &m_search_dummy,
               m_compare_less);
  }

  void initialize(std::vector<SurfacePoint>& sources) //we initialize the sources by copie
  {
    resize(sources.size());
    m_sorted.resize(sources.size());
    for(unsigned i=0; i<sources.size(); ++i)
    {
      SurfacePointWithIndex& p = *(begin() + i);

      p.initialize(sources[i],i);
      m_sorted[i] = &p;
    }

    std::sort(m_sorted.begin(), m_sorted.end(), m_compare_less);
  };

  SurfacePointWithIndex& operator[](unsigned i)
  {
    assert(i < size());
    return *(begin() + i);
  }

private:
  sorted_vector_type m_sorted;
  SurfacePointWithIndex m_search_dummy;   //used as a search template
  SurfacePointWithIndex m_compare_less;     //used as a compare functor
};


inline void Interval::find_closest_point(double const rs, 
                     double const hs,
                     double& r,
                     double& d_out)     //find the point on the interval that is closest to the point (alpha, s)
  {
    if(m_d == GEODESIC_INF)
    {
      r = GEODESIC_INF;
      d_out = GEODESIC_INF;
      return;
    }

    double hc = -m_pseudo_y;
    double rc = m_pseudo_x;
    double end = stop();

    double local_epsilon = SMALLEST_INTERVAL_RATIO*m_edge->length();
    if(std::abs(hs+hc) < local_epsilon)
    {
      if(rs<=m_start)
      {
        r = m_start;
        d_out = signal(m_start) + std::abs(rs - m_start);
      }
      else if(rs>=end)
      {
        r = end;
        d_out = signal(end) + fabs(end - rs);
      }
      else
      {
        r = rs;
        d_out = signal(rs);
      }
    }
    else
    {
      double ri = (rs*hc + hs*rc)/(hs+hc);

      if(ri<m_start)
      {
        r = m_start;
        d_out = signal(m_start) + hypotenuse(m_start - rs, hs);
      }
      else if(ri>end)
      {
        r = end;
        d_out = signal(end) + hypotenuse(end - rs, hs);
      }
      else
      {
        r = ri;
        d_out = m_d + hypotenuse(rc - rs, hc + hs);
      }
    }
  }


inline void Interval::initialize(edge_pointer edge, 
                 SurfacePoint* source,
                 unsigned source_index)
{
  m_next = NULL;
  //m_geodesic_previous = NULL;
  m_direction = UNDEFINED_DIRECTION;
  m_edge = edge;
  m_source_index = source_index;

  m_start = 0.0;
  //m_stop = edge->length();
  if(!source)
  {
    m_d = GEODESIC_INF;
    m_min = GEODESIC_INF;
    return;
  }
  m_d = 0;

  if(source->base_element()->type() == VERTEX)
  {
    if(source->base_element()->id() == edge->v0()->id())
    {
      m_pseudo_x = 0.0;
      m_pseudo_y = 0.0;
      m_min = 0.0;
      return;
    }
    else if(source->base_element()->id() == edge->v1()->id())
    {
      m_pseudo_x = stop();
      m_pseudo_y = 0.0;
      m_min = 0.0;
      return;
    }
  }

  edge->local_coordinates(source, m_pseudo_x, m_pseudo_y);
  m_pseudo_y = -m_pseudo_y;

  compute_min_distance(stop());
} 

}   //geodesic

#endif //GEODESIC_ALGORITHM_EXACT_ELEMENTS_20071231