/Users/barrand/private/dev/softinex/old/inexlib-1.2/inlib/inlib/dct

Go to the documentation of this file.

// Copyright (C) 2010, Guy Barrand. All rights reserved.
// See the file inlib.license for terms.

#ifndef inlib_dct
#define inlib_dct

// A Delaunay triangulation algorithm.
// From Geoff Leach dct C code. See credit at end.
// C++ version done by G.Barrand in June 2010.

#include <vector>
#include <string>
#include <ostream>

namespace inlib {
namespace dct {

typedef float real;

class edge;

class point {
public:
  real x,y;
  edge* entry_pt;
};

class edge {
public:
  point* org;
  point* dest;
  edge* onext;
  edge* oprev;
  edge* dnext;
  edge* dprev;
};

class triangle {
public:
  triangle(unsigned int a_1,unsigned int a_2,unsigned int a_3)
  :m_1(a_1),m_2(a_2),m_3(a_3)
  {}
public:
  triangle(const triangle& a_from)
  :m_1(a_from.m_1),m_2(a_from.m_2),m_3(a_from.m_3)
  {}
  triangle& operator=(const triangle& a_from) {
    m_1 = a_from.m_1;
    m_2 = a_from.m_2;
    m_3 = a_from.m_3;
    return *this;
  }
  unsigned int index_1() const { return m_1;}
  unsigned int index_2() const { return m_2;}
  unsigned int index_3() const { return m_3;}
private:
  unsigned int m_1;
  unsigned int m_2;
  unsigned int m_3;
};

class alg {
public:
  typedef enum {right, left} side;
  typedef unsigned int index_t; /* G.Barrand : have _t. */

public:
  alg(std::ostream& a_out)
  :m_out(a_out)
  ,p_number(0),p_array(0)
  ,e_array(0)
  ,free_list_e(0)
  ,n_free_e(0)
  {}
  virtual ~alg(){
    delete [] p_array;
    delete [] e_array;  
    delete [] free_list_e;  
  }
private:
  alg(const alg& a_from):m_out(a_from.m_out){}
  alg& operator=(const alg&){return *this;}
public:
  bool process(const std::vector< std::pair<real,real> >& a_data){
    delete [] p_array;
    delete [] e_array;  
    delete [] free_list_e;  

    p_number = a_data.size();
    //alloc :
    p_array = new point[p_number];
    // edges :
    n_free_e = 3 * p_number;   /* Eulers relation */
    e_array = new edge[n_free_e];
    if(!e_array) {
      panic("Not enough memory\n");
      return false;
    }
    typedef edge* ept;
    free_list_e = new ept[n_free_e];
    if(!free_list_e) {
      panic("Not enough memory\n");
      return false;
    }
   {edge* e = e_array;
    for(unsigned int i=0;i<n_free_e;i++,e++) free_list_e[i] = e;}
  
    point* pos = p_array;
    std::vector< std::pair<real,real> >::const_iterator it;  
    for (it=a_data.begin();it!=a_data.end();++it,++pos) {
      (*pos).x = (*it).first;
      (*pos).y = (*it).second;
      (*pos).entry_pt = 0;
    }


    // sort :
    typedef point* ppt;
    ppt* p_sorted = new ppt[p_number];
    if(!p_sorted) {
      panic("Not enough memory\n");
      return false;
    }

    ppt* p_temp = new ppt[p_number];
    if(!p_temp) {
      panic("Not enough memory\n");
      return false;
    }
   {for (unsigned int i = 0; i < p_number; i++) p_sorted[i] = p_array + i;}
    merge_sort(p_sorted, p_temp, 0, p_number-1);  
    delete [] p_temp;


    // triangulate :
    edge* l_cw;
    edge* r_ccw;
    divide(p_sorted, 0, p_number-1, &l_cw, &r_ccw);
  
    delete [] p_sorted;

    return true;
  }

  std::vector<triangle> get_triangles() {
    std::vector<triangle> ts;
    //  Print the ring of triangles about each vertex.
    edge *e_start, *e, *next;
    point *u, *v, *w;
    index_t i;
    point *t;
  
    for (i = 0; i < p_number; i++) {
      u = &p_array[i];
      e_start = e = u->entry_pt;
      do
      {
        v = Other_point(e, u);
        if (u < v) {
        next = Next(e, u);
        w = Other_point(next, u);
        if (u < w)
          if (Identical_refs(Next(next, w), Prev(e, v))) {  
            /* Triangle. */
            if (v > w) { t = v; v = w; w = t; }

            ts.push_back(triangle((unsigned int)(u - p_array),
                                  (unsigned int)(v - p_array),
                                  (unsigned int)(w - p_array)));
          }
        }
  
        /* Next edge around u. */
        e = Next(e, u);
      } while (!Identical_refs(e, e_start));
    }

    return ts;
  }

  void print_triangles() {
    //  Print the ring of triangles about each vertex.
    edge *e_start, *e, *next;
    point *u, *v, *w;
    index_t i;
    point *t;
  
    for (i = 0; i < p_number; i++) {
      u = &p_array[i];
      e_start = e = u->entry_pt;
      do
      {
        v = Other_point(e, u);
        if (u < v) {
        next = Next(e, u);
        w = Other_point(next, u);
        if (u < w)
          if (Identical_refs(Next(next, w), Prev(e, v))) {  
            /* Triangle. */
            if (v > w) { t = v; v = w; w = t; }
            m_out << (unsigned long)(u - p_array)
                  << " " << (unsigned long)(v - p_array)
                  << " " << (unsigned long)(w - p_array)
                  << std::endl;
          }
        }
  
        /* Next edge around u. */
        e = Next(e, u);
      } while (!Identical_refs(e, e_start));
    }
  }

  void divide(point *p_sorted[], index_t l, index_t r, edge **l_ccw, edge **r_cw)
  {
    unsigned int n;
    index_t split;
    edge *l_ccw_l, *r_cw_l, *l_ccw_r, *r_cw_r, *l_tangent;
    edge *a, *b, *c;
    real c_p;
    
    n = r - l + 1;
    if (n == 2) {
      /* Bottom of the recursion. Make an edge */
      *l_ccw = *r_cw = make_edge(p_sorted[l], p_sorted[r]);
    } else if (n == 3) {
      /* Bottom of the recursion. Make a triangle or two edges */
      a = make_edge(p_sorted[l], p_sorted[l+1]);
      b = make_edge(p_sorted[l+1], p_sorted[r]);
      splice(a, b, p_sorted[l+1]);
      c_p = Cross_product_3p(p_sorted[l], p_sorted[l+1], p_sorted[r]);
      
      if (c_p > 0.0)
      {
        /* Make a triangle */
        c = join(a, p_sorted[l], b, p_sorted[r], right);
        *l_ccw = a;
        *r_cw = b;
      } else if (c_p < 0.0) {
        /* Make a triangle */
        c = join(a, p_sorted[l], b, p_sorted[r], left);
        *l_ccw = c;
        *r_cw = c;
      } else {
        /* Points are collinear,  no triangle */ 
        *l_ccw = a;
        *r_cw = b;
      }
    } else if (n  > 3) {
      /* Continue to divide */
  
      /* Calculate the split point */
      split = (l + r) / 2;
    
      /* Divide */
      divide(p_sorted, l, split, &l_ccw_l, &r_cw_l);
      divide(p_sorted, split+1, r, &l_ccw_r, &r_cw_r);
  
      /* Merge */
      merge(r_cw_l, p_sorted[split], l_ccw_r, p_sorted[split+1], &l_tangent);
  
      /* The lower tangent added by merge may have invalidated 
         l_ccw_l or r_cw_r. Update them if necessary. */
      if (Org(l_tangent) == p_sorted[l])
        l_ccw_l = l_tangent;
      if (Dest(l_tangent) == p_sorted[r])
        r_cw_r = l_tangent;
  
      /* Update edge refs to be passed back */ 
      *l_ccw = l_ccw_l;
      *r_cw = r_cw_r;
    }
  }
  
private:
  /*
   *  Determines the lower tangent of two triangulations. 
   */
  void lower_tangent(edge *r_cw_l, point *s, edge *l_ccw_r, point *u,
                            edge **l_lower, point **org_l_lower,
                            edge **r_lower, point **org_r_lower)
  {
    edge *l, *r;
    point *o_l, *o_r, *d_l, *d_r;
    bool finished;
  
    l = r_cw_l;
    r = l_ccw_r;
    o_l = s;
    d_l = Other_point(l, s);
    o_r = u;
    d_r = Other_point(r, u);
    finished = false;
  
    while (!finished)
      if (Cross_product_3p(o_l, d_l, o_r) > 0.0) {
        l = Prev(l, d_l);
        o_l = d_l;
        d_l = Other_point(l, o_l);
      } else if (Cross_product_3p(o_r, d_r, o_l) < 0.0) {
        r = Next(r, d_r);
        o_r = d_r;
        d_r = Other_point(r, o_r);
      } else
        finished = true;
  
    *l_lower = l;
    *r_lower = r;
    *org_l_lower = o_l;
    *org_r_lower = o_r;
  }
  
  /* 
   *  The merge function is where most of the work actually gets done.  It is
   *  written as one (longish) function for speed.
   */ 
  void merge(edge *r_cw_l, point *s, edge *l_ccw_r, point *u, edge **l_tangent)
  {
    edge *base, *l_cand, *r_cand;
    point *org_base, *dest_base;
    real u_l_c_o_b, v_l_c_o_b, u_l_c_d_b, v_l_c_d_b;
    real u_r_c_o_b, v_r_c_o_b, u_r_c_d_b, v_r_c_d_b;
    real c_p_l_cand, c_p_r_cand;
    real d_p_l_cand, d_p_r_cand;
    bool above_l_cand, above_r_cand, above_next, above_prev;
    point *dest_l_cand, *dest_r_cand;
    real cot_l_cand, cot_r_cand;
    edge *l_lower, *r_lower;
    point *org_r_lower, *org_l_lower;
  
    /* Create first cross edge by joining lower common tangent */
    lower_tangent(r_cw_l, s, l_ccw_r, u, &l_lower, &org_l_lower, &r_lower, &org_r_lower);
    base = join(l_lower, org_l_lower, r_lower, org_r_lower, right);
    org_base = org_l_lower;
    dest_base = org_r_lower;
    
    /* Need to return lower tangent. */
    *l_tangent = base;
  
    /* Main merge loop */
    do 
    {
      /* Initialise l_cand and r_cand */
      l_cand = Next(base, org_base);
      r_cand = Prev(base, dest_base);
      dest_l_cand = Other_point(l_cand, org_base);
      dest_r_cand = Other_point(r_cand, dest_base);
  
      /* Vectors for above and "in_circle" tests. */
      Vector(dest_l_cand, org_base, u_l_c_o_b, v_l_c_o_b);
      Vector(dest_l_cand, dest_base, u_l_c_d_b, v_l_c_d_b);
      Vector(dest_r_cand, org_base, u_r_c_o_b, v_r_c_o_b);
      Vector(dest_r_cand, dest_base, u_r_c_d_b, v_r_c_d_b);
  
      /* Above tests. */
      c_p_l_cand = Cross_product_2v(u_l_c_o_b, v_l_c_o_b, u_l_c_d_b, v_l_c_d_b);
      c_p_r_cand = Cross_product_2v(u_r_c_o_b, v_r_c_o_b, u_r_c_d_b, v_r_c_d_b);
      above_l_cand = c_p_l_cand > 0.0;
      above_r_cand = c_p_r_cand > 0.0;
      if (!above_l_cand && !above_r_cand)
        break;        /* Finished. */
  
      /* Advance l_cand ccw,  deleting the old l_cand edge,  until the 
         "in_circle" test fails. */
      if (above_l_cand)
      {
        real u_n_o_b, v_n_o_b, u_n_d_b, v_n_d_b;
        real c_p_next, d_p_next, cot_next;
        edge *next;
        point *dest_next;
  
        d_p_l_cand = Dot_product_2v(u_l_c_o_b, v_l_c_o_b, u_l_c_d_b, v_l_c_d_b);
        cot_l_cand = d_p_l_cand / c_p_l_cand;
  
        do 
        {
          next = Next(l_cand, org_base);
          dest_next = Other_point(next, org_base);
          Vector(dest_next, org_base, u_n_o_b, v_n_o_b);
          Vector(dest_next, dest_base, u_n_d_b, v_n_d_b);
          c_p_next = Cross_product_2v(u_n_o_b, v_n_o_b, u_n_d_b, v_n_d_b);
          above_next = c_p_next > 0.0;
  
          if (!above_next) 
            break;    /* Finished. */
  
          d_p_next = Dot_product_2v(u_n_o_b, v_n_o_b, u_n_d_b, v_n_d_b);
          cot_next = d_p_next / c_p_next;
  
          if (cot_next > cot_l_cand)
            break;    /* Finished. */
  
          delete_edge(l_cand);
          l_cand = next;
          cot_l_cand = cot_next;
    
        } while (true);
      }
  
      /* Now do the symmetrical for r_cand */
      if (above_r_cand)
      {
        real u_p_o_b, v_p_o_b, u_p_d_b, v_p_d_b;
        real c_p_prev, d_p_prev, cot_prev;
        edge *prev;
        point *dest_prev;
  
        d_p_r_cand = Dot_product_2v(u_r_c_o_b, v_r_c_o_b, u_r_c_d_b, v_r_c_d_b);
        cot_r_cand = d_p_r_cand / c_p_r_cand;
  
        do
        {
          prev = Prev(r_cand, dest_base);
          dest_prev = Other_point(prev, dest_base);
          Vector(dest_prev, org_base, u_p_o_b, v_p_o_b);
          Vector(dest_prev, dest_base, u_p_d_b, v_p_d_b);
          c_p_prev = Cross_product_2v(u_p_o_b, v_p_o_b, u_p_d_b, v_p_d_b);
          above_prev = c_p_prev > 0.0;
  
          if (!above_prev) 
            break;    /* Finished. */
  
          d_p_prev = Dot_product_2v(u_p_o_b, v_p_o_b, u_p_d_b, v_p_d_b);
          cot_prev = d_p_prev / c_p_prev;
  
          if (cot_prev > cot_r_cand)
            break;    /* Finished. */
  
          delete_edge(r_cand);
          r_cand = prev;
          cot_r_cand = cot_prev;
  
        } while (true);
      }
  
      /*
       *  Now add a cross edge from base to either l_cand or r_cand. 
       *  If both are valid choose on the basis of the in_circle test . 
       *  Advance base and  whichever candidate was chosen.
       */
      dest_l_cand = Other_point(l_cand, org_base);
      dest_r_cand = Other_point(r_cand, dest_base);
      if (!above_l_cand || (above_l_cand && above_r_cand && cot_r_cand < cot_l_cand))
      {
        /* Connect to the right */
        base = join(base, org_base, r_cand, dest_r_cand, right);
        dest_base = dest_r_cand;
      } else {
        /* Connect to the left */
        base = join(l_cand, dest_l_cand, base, dest_base, right);
        org_base = dest_l_cand;
      }
  
    } while (true);
  }

private:
  edge *make_edge(point *u, point *v)
  {
    edge *e;
  
    e = get_edge();
  
    e->onext = e->oprev = e->dnext = e->dprev = e;
    e->org = u;
    e->dest = v;
    if (u->entry_pt == NULL)
      u->entry_pt = e;
    if (v->entry_pt == NULL)
      v->entry_pt = e;
    
  
    return e;
  }
  
  /* 
   *  Add an edge to a ring of edges. 
   */
  void splice(edge *a, edge *b, point *v)
  {
    edge *next;
  
    
    /* b must be the unnattached edge and a must be the previous 
       ccw edge to b. */
  
    if (Org(a) == v) { 
      next = Onext(a);
      Onext(a) = b;
    } else {
      next = Dnext(a);
      Dnext(a) = b;
    }
  
    if (Org(next) == v)
      Oprev(next) = b;
    else
      Dprev(next) = b;
  
    if (Org(b) == v) {
      Onext(b) = next;
      Oprev(b) = a;
    } else {
      Dnext(b) = next;
      Dprev(b) = a;
    }
  }
  
  /* 
   *  Creates a new edge and adds it to two rings of edges.
   */
  edge *join(edge *a, point *u, edge *b, point *v, side s)
  {
    edge *e;
  
    /* u and v are the two vertices which are being joined.
       a and b are the two edges associated with u and v res.  */
  
    e = make_edge(u, v);
    
    if (s == left) {
      if (Org(a) == u)
        splice(Oprev(a), e, u);
      else
        splice(Dprev(a), e, u);
      splice(b, e, v);
    } else {
      splice(a, e, u);
      if (Org(b) == v)
        splice(Oprev(b), e, v);
      else
        splice(Dprev(b), e, v);
    }
  
    return e;
  }
  
  /* 
   *  Remove an edge.
   */
  void delete_edge(edge *e)
  {
    point *u, *v;
  
    /* Cache origin and destination. */
    u = Org(e);
    v = Dest(e);
  
    /* Adjust entry points. */
    if (u->entry_pt == e)
      u->entry_pt = e->onext;
    if (v->entry_pt == e)
      v->entry_pt = e->dnext;
  
     /* Four edge links to change */
    if (Org(e->onext) == u)
      e->onext->oprev = e->oprev;
    else
      e->onext->dprev = e->oprev;
  
    if (Org(e->oprev) == u)
      e->oprev->onext = e->onext;
    else
      e->oprev->dnext = e->onext;
  
    if (Org(e->dnext) == v)
      e->dnext->oprev = e->dprev;
    else
      e->dnext->dprev = e->dprev;
  
    if (Org(e->dprev) == v)
      e->dprev->onext = e->dnext;
    else
      e->dprev->dnext = e->dnext;
  
    free_edge(e);
  }

  edge *get_edge()
  {
    if (n_free_e == 0)
      panic("Out of memory for edges\n");
  
     return (free_list_e[--n_free_e]);
  }
  
  void free_edge(edge *e)
  {
     free_list_e[n_free_e++] = e;
  }

  static void merge_sort(point *p[], point *p_temp[], index_t l, index_t r)
  {
    index_t i, j, k, m;
  
    if (r - l > 0)
    {
      m = (r + l) / 2;
      merge_sort(p, p_temp, l, m);
      merge_sort(p, p_temp, m+1, r);
  
      for (i = m+1; i > l; i--)
        p_temp[i-1] = p[i-1];
      for (j = m; j < r; j++)
        p_temp[r+m-j] = p[j+1];
      for (k = l; k <= r; k++)
        if (p_temp[i]->x < p_temp[j]->x) {
          p[k] = p_temp[i];
          i = i + 1;
        } else if (p_temp[i]->x == p_temp[j]->x && p_temp[i]->y < p_temp[j]->y) {
          p[k] = p_temp[i];
          i = i + 1;
        } else {
          p[k] = p_temp[j];
          j = j - 1;
        }
    }
  }

  void panic(const char *m){m_out << std::string(m);}

private:
  static point*& Org(edge* e) {return e->org;}
  static point*& Dest(edge* e) {return e->dest;}
  static edge*& Onext(edge* e) {return e->onext;}
  static edge*& Oprev(edge* e) {return e->oprev;}
  static edge*& Dnext(edge* e) {return e->dnext;}
  static edge*& Dprev(edge* e) {return e->dprev;}

  static point* Other_point(edge* e,point* p) {
    return (e->org == p ? e->dest : e->org);
  }
  static edge* Next(edge* e,point* p) {
    return (e->org == p ? e->onext : e->dnext);
  }
  static edge* Prev(edge* e,point* p) {
    return (e->org == p ? e->oprev : e->dprev);
  }

  static void Vector(point* p1,point* p2,real& u,real& v) {
    u = p2->x - p1->x;
    v = p2->y - p1->y;
  }
  static real Cross_product_2v(real u1,real v1,real u2,real v2) {
    return u1 * v2 - v1 * u2;
  }
  static real Cross_product_3p(point* p1,point* p2,point* p3) {
    return (p2->x - p1->x)*(p3->y - p1->y) - (p2->y - p1->y)*(p3->x - p1->x);
  }

  static real Dot_product_2v(real u1,real v1,real u2,real v2) {
    return u1 * u2 + v1 * v2;
  }

  static bool Identical_refs(edge* e1,edge* e2) {
    return (e1==e2 ? true:false);
  }

private:
  std::ostream& m_out;
  unsigned int p_number;
  point* p_array;
  edge* e_array;
  edge** free_list_e;
  unsigned int n_free_e;
};

}}

#endif

/*   Author: Geoff Leach, Department of Computer Science, RMIT.
 *   email: gl@cs.rmit.edu.au
 *
 *   Date: 6/10/93
 *
 *   Version 1.0
 *   
 *   Copyright (c) RMIT 1993. All rights reserved.
 *
 *   License to copy and use this software purposes is granted provided 
 *   that appropriate credit is given to both RMIT and the author.
 *
 *   License is also granted to make and use derivative works provided
 *   that appropriate credit is given to both RMIT and the author.
 *
 *   RMIT makes no representations concerning either the merchantability 
 *   of this software or the suitability of this software for any particular 
 *   purpose.  It is provided "as is" without express or implied warranty 
 *   of any kind.
 *
 *   These notices must be retained in any copies of any part of this software.
 */