Delaunay1.h

Go to the documentation of this file.

#ifndef Delaunay1_HEADER
#define Delaunay1_HEADER

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif

#include <topology/Complex1.h>
#include "IncrementalDelaunay.h"

#undef DEBUG_DELAUNAY
#define DEBUG_DELAUNAY 1

template <unsigned ambient,
  class Vertex_, class SimplexData,
  class VertexPrinter>
class IncrementalDelaunay<ambient, 1, Vertex_, SimplexData, VertexPrinter> {
public:
static const unsigned topological = 1;
typedef Vertex_ Vertex;
typedef Geometry::Point<ambient> Point; 
typedef Geometry::Point<topological> PPoint; 
typedef Geometry::ProjectionPlane<ambient, topological> ProjectionPlane;
typedef Geometry::BoundingBox<topological> BoundingBox;
typedef Geometry::CenterRadius<ambient> CenterRadius;
typedef SimplicialComplex<topological, Vertex_, SimplexData, VertexPrinter> Complex;
typedef typename Complex::Simplex Simplex;
typedef typename Complex::Cavity Cavity;
typedef typename Complex::Star Star;

/***************************************************************************
 ***************************************************************************/
public:
IncrementalDelaunay(const ProjectionPlane& plane, Vertex *a, Vertex *b)
  : plane_(plane)
  , complex_(isOriented(a,b) ? a : b, isOriented(a,b) ? b : a)
{
  assert(isOriented());
}

static IncrementalDelaunay *fromSimplex(const ProjectionPlane& plane, Vertex *a, Vertex *b) {
  return new IncrementalDelaunay(plane, a, b);
}

static IncrementalDelaunay *fromSimplex(Vertex *a, Vertex *b) {
  boost::array<Point, 2> p;
  p[0] = a->toPoint();
  p[1] = b->toPoint();
  return new IncrementalDelaunay(ProjectionPlane(p.begin(), p.end()), a, b);
}

static IncrementalDelaunay *fromSimplex(const vector<Vertex*>& verts) {
  assert(verts.size() == 2);
  return fromSimplex(verts[0], verts[1]);
}

template<class vertex_iterator>
static IncrementalDelaunay *computeDelaunaySlowlyFromVertices (vertex_iterator begin, vertex_iterator end)
{
  assert(begin != end);
  Vertex *a = *begin;
  ++begin;
  assert(begin != end);
  Vertex *b = *begin;
  assert(++begin == end);

  return fromSimplex(a, b);
}

static IncrementalDelaunay *makeBoundingSimplex(
    const ProjectionPlane& plane,
    BoundingBox box, // copy
    double buffer) {
  if (buffer != 0.0) { box = box.scale(buffer); }
  box.scale(buffer);
  Vertex *a = new Vertex(plane.unproject(box.mincorner));
  Vertex *b = new Vertex(plane.unproject(box.mincorner + PPoint(box.extent)));
  IncrementalDelaunay *id = new IncrementalDelaunay(plane, a, b);
  id->myVertices_.push_back(a);
  id->myVertices_.push_back(b);
  return id;
}

template <class point_iterator>
static IncrementalDelaunay *makeBoundingSimplex(point_iterator begin, point_iterator end, double buffer) {
  ProjectionPlane plane(begin, end);
  BoundingBox box(begin, end, plane);
  return makeBoundingSimplex(plane, box, buffer);
}

/***************************************************************************
 ***************************************************************************/
const Complex& getComplex() const {
  return complex_;
}
Complex& getComplexRW() {
  return complex_;
}
const Simplex& getHandle() const {
  return complex_.getHandle();
}

/***************************************************************************
 ***************************************************************************/
#ifdef NDEBUG
#define unused(p)
#else
#define unused(p) p
#endif
void computeCavity(const Simplex& s, const Point& unused(p), Cavity& cavity, bool append = false) const {
  // The cavity in 1-d is just the single simplex.
  if (!append) { assert(cavity.empty()); }
  assert(inSimplexAffine(s, p));
  cavity.add(s);
}

void computeCavity(const Simplex& s, const Point& unused(p), Cavity& cavity, vector<Simplex>& neighbours, bool append = false) const {
#undef unused
  assert(isMember(s));
  if (!append) {
    assert(cavity.empty());
    assert(neighbours.empty());
  }
  assert(inSimplexAffine(s, p));
  cavity.add(s);

  // Get both (or just one, or none) the neighbours of s.
  Simplex n = complex_.getNeighbour(0, s);
  if (complex_.isMember(n)) { neighbours.push_back(n); }
  n = complex_.getNeighbour(1, s);
  if (complex_.isMember(n)) { neighbours.push_back(n); }
}

void computeStar(const Cavity& cavity, Vertex *v, Star& star) {
  assert(cavity.size() == 1);
  complex_.computeStar(*cavity.begin(), v, star);
}

void insertVertex(const Cavity& cavity, Star& star) {
  complex_.replaceCavity(cavity, star);
}

void insertVertex(const Simplex& s, Vertex *v) {
  complex_.replaceCavity(s, v);
}
void insertVertex(const Cavity& cavity, Vertex *v) {
  complex_.replaceCavity(cavity, v);
}

void expensiveInsertVertex(Vertex *v) {
  insertVertex(expensiveLocate(v), v);
}

Simplex expensiveLocate(const Vertex *v) {
  return expensiveLocate(v->toPoint());
}

Simplex expensiveLocate(const Point& p) {
  Simplex s = getHandle();
  double q = x(p);
  double a = x(s[0]);
  double b = x(s[1]);
  assert(a < b);

  if (q > a) {
    // go right until we find our simplex (could be right here!)
    if (q < b) {
      return s;
    } else {
      while (q > b) {
        s = complex_.getNeighbour(0, s);
        assert(x(s[0]) == b);
        // a = x(s[0]);
        b = x(s[1]);
      }
      return s;
    }
  } else {
    // go left until we find our simplex
    while (q < a) {
      s = complex_.getNeighbour(1, s);
      assert(x(s[1]) == a);
      a = x(s[0]);
      // b = x(s[1]);
    }
    return s;
  }
}


/***************************************************************************
 ***************************************************************************/
private:
double x(const Vertex *v) const {
  return plane_.project(v->toPoint())[0];
}
double x(const Point& p) const {
  return plane_.project(p)[0];
}

public:
const ProjectionPlane& getPlane() const {
  return plane_;
}

CenterRadius circumcenter(const Simplex& s) const {
  assert(isMember(s));
  return Geometry::template circumcenter_segment<ambient>(s[0]->toPoint(), s[1]->toPoint());
}

static double radiusEdge2(const Simplex&) {
  return 1.0;
}

static double computeSigma(const Simplex&) {
  return 1.0;
}

static double radiusRadius(const Simplex&) {
  return 1.0;
}

/***************************************************************************
 ***************************************************************************/
bool isMember(const Simplex& s) const {
  return getComplex().isMember(s);
}

bool inSphere(const Simplex& s, const Point& p) const {
  // TODO: make robust
  CenterRadius cr = circumcenter(s);
  return p.dist2(cr.center) < cr.radius2;
}

bool inSphereAffine(const Simplex& s, const Point& p) const {
  double a = x(s[0]);
  double b = x(s[1]);
  double q = x(p);
  assert (a < b);
  return a < q && q < b;
}

bool inSimplexAffine(const Simplex& s, const Point& p) const {
  return inSphereAffine(s, p);
}

bool isOriented(const Simplex& s) const {
  return isOriented(s[0], s[1]);
}
bool isOriented(const Vertex *a, const Vertex *b) const {
  bool isO = x(a) < x(b);
#if DEBUG_DELAUNAY
  if (!isO) {
    printf("Unoriented segment: Vertices %s and %s (at %s and %s)\n",
        a->toString().c_str(), b->toString().c_str(),
        a->toPoint().toString().c_str(), b->toPoint().toString().c_str());
    printf("Project to %g and %g\n", x(a), x(b));
  }
#endif
  return isO;
}

private:
struct CheckOriented {
  const IncrementalDelaunay *del;
  CheckOriented(const IncrementalDelaunay *del): del(del) { }
  static bool canEnter(const Simplex&) { return true; }
  bool choose(const Simplex& s) { return !del->isOriented(s); }
};
public:
bool isOriented() const {
  CheckOriented co(this);
  return !complex_.dfsBySimplex(co, getHandle());
}
static bool isDelaunay() { return true; }

/***************************************************************************
 ***************************************************************************/
public:

/***************************************************************************
 * \name Data access
 * We have to be careful with what we return, in order to handle void data.
 * We forward the request on to the complex, and then we pull out the
 * user's part of it.
 ***************************************************************************/
typedef typename hudson::Bucket<SimplexData>::reference data_ref;
typedef typename hudson::Bucket<SimplexData>::const_reference data_const_ref;
data_ref getDataRW(const Simplex& s) {
  return complex_.getDataRW(s);
}

data_const_ref getData(const Simplex& s) const {
  return complex_.getData(s);
}



/***************************************************************************
 ***************************************************************************/

private:
ProjectionPlane plane_; // MUST be before complex_
Complex complex_;
std::vector<Vertex*> myVertices_;

};

#endif

---