PointRefineVertex.h

Go to the documentation of this file.

#ifndef PointRefineVertex_HEADER
#define PointRefineVertex_HEADER

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

#include <geometry/Point.h>
#include <delaunay/PointVertex.h>
#include <utilities/IntrusiveList.h>
#include <utilities/SmallList.h>
#include <refine-common/ConcentricShells.h>

// #undef MAINTAIN_NN

#if 0
#define dprintf(...) printf(__VA_ARGS__);
#define dputs(str) puts(str);
#else
#define dprintf(...)
#define dputs(str)
#endif



/***************************************************************************
 * The vertex class used by the refinement.
 *
 * I use the vertices to store the uninserted points (namely, those that lie
 * in their Voronoi cells), hence the need to specialize the Vertex class.
 *
 * Theoretically, there is an infinite set of annuli around each vertex,
 * with radii in geometric progression.  We store all the annuli that are
 * non-empty, along with their outradius and the points they contain.  This
 * means we look much (constant factors) less frequently at input points than
 * if we stored all the points in a single bucket.
 *
 * There are some facts we do not use: first, note that Steiner points have all
 * their points in O(1) annuli (since there are no points allowed close to a
 * Steiner point).  Second, annuli don't distinguish in which direction points
 * lie: a search from the left will investigate points that are very far to the
 * right.
 *
 * The uninserted points are stored in a number of lists.  Rather than
 * allocating and freeing list buckets as we reassign points, we use an
 * IntrusiveList: the list buckets are the vertices themselves.  This saves
 * a bit of memory allocation.
 ***************************************************************************/
template <unsigned d, class Simplex>
class PointRefineVertex : public PointVertex<d> {

  typedef Geometry::Point<d> Point;
  typedef PointRefineVertex Uninserted;

  public:
  PointRefineVertex(const Point& p, size_t ID = size_t(-1)) : PointVertex<d>(p) {
    visited_ = 0;
    inMesh_ = false;
    ID_ = ID;
    isSteiner_ = (ID == size_t(-1));
#ifdef MAINTAIN_NN
    NN2_ = std::numeric_limits<double>::infinity();
#endif
  }

  ~PointRefineVertex() {
    if (inMesh()) {
      get().~Simplex();
    }
  }

#ifdef NDEBUG
# define unused(sym)
#else
# define unused(sym) sym
#endif

  void *operator new(size_t unused(sz)) {
    assert(sz == sizeof(PointRefineVertex));
    return hudson::SmallMemoryPool::allocate<sizeof(PointRefineVertex)>();
  }

  void *operator new(size_t unused(sz), void *ptr) {
    assert(sz == sizeof(PointRefineVertex));
    return ptr;
  }

  void operator delete(void *v, size_t unused(sz)) {
    assert(sz == sizeof(PointRefineVertex));
    hudson::SmallMemoryPool::deallocate<sizeof(PointRefineVertex)>(v);
  }
#undef unused


  /************************************************************
   * \name Name.
   ************************************************************/
  size_t id() const { return ID_; }
  bool isSteiner() const { return isSteiner_; }
  void setID(size_t i) { ID_ = i; }

  /************************************************************
   * \name Distance.
   ************************************************************/
  double dist2(const PointRefineVertex *other) const {
    return this->toPoint().dist2(other->toPoint());
  }
  double dist2(const Point& other) const {
    return this->toPoint().dist2(other);
  }

#ifdef MAINTAIN_NN
  double getNN2() const {
    return NN2_;
  }
#endif

  /************************************************************
   * \name Queries and relocation
   * Foreach annulus, from outermost to innermost:
   *   if R(a) < |pq| - r then quit
   *   else iterate through the points.
   * r is |pq|/2 in the reassign case; it's the best point in
   * the query case.
   *
   * The newv is presumed to be a vertex which may be the new nearest
   * neighbour.
   ************************************************************/
  void reassign(PointRefineVertex *newv) {
#ifdef MAINTAIN_NN
    double pq2 = dist2(newv);
    if (pq2 < NN2_) { NN2_ = pq2; }
#endif

    annuli_.reassignTo(this->toPoint(), newv->annuli_, newv->toPoint());
  }

  pair<Uninserted*, double> findClosest(const Point& q, double r2) const {
    return annuli_.findClosest(this->toPoint(), q, r2);
  }

  pair<Uninserted*, double> findPoint(const Point& q, double r2) const {
    return annuli_.findPoint(this->toPoint(), q, r2);
  }

  void assign(Uninserted *u) {
    // If we were uninserted and just got inserted, we shall be assigned
    // ourselves.  Simply ignore the request.
    if (u == this) return;
    annuli_.add(this->toPoint(), u);
  }

  /************************************************************
   * Get a far-away point to insert.
   ************************************************************/
  bool isCrowded() const {
    return !annuli_.empty();
  }
  Uninserted *getFarPoint() const {
    return annuli_.getFarItem();
  }

  /************************************************************
   * \name Back-pointer
   * We need a handle from the vertex back to the mesh.
   ************************************************************/
  const Simplex& getHandle() const { return get(); }
  void setHandle(const Simplex& s) { get() = s; }
  void initHandle(const Simplex& s) {
    assert(!inMesh());
    new (&getUnchecked()) Simplex(s);
    inMesh_ = true;
  }

  /************************************************************
   * \name Visited flags
   *
   * We often want to visit vertices.  We only ever have one search
   * at a time.  Therefore, we simply store a flag in the vertex.
   * The flag is simply cleared every 2^32 searches (which typically
   * means we never bother cleaning the flag).
   ************************************************************/
  typedef unsigned FlagType;
  FlagType visitedFlag() const { return visited_; }
  void setVisited(FlagType i) { visited_ = i; }
  void clearVisited() { visited_ = 0; }


  private:

  /************************************************************
   * \name Management of the simplex handle
   *
   * Since we don't have a no-arg ctor for Simplex, and since
   * we don't get a simplex at Vertex ctor time, we have to
   * hack around it.  Paying for a scoped_ptr is too expensive.
   * Implications:
   * 1. we use placement new to initialize the simplex, in initData
   * 2. we have to explicitely destroy the simplex in ~Vertex
   * 3. use get() to access the simplex.
   ************************************************************/
  bool inMesh() const { return inMesh_; }
  Simplex& getUnchecked() { return *(Simplex*)&handle_[0]; }
  Simplex& get() {
    assert(inMesh());
    return getUnchecked();
  }
  const Simplex& get() const {
    assert(inMesh());
    return *(Simplex*)&handle_[0];
  }


  /************************************************************
   * Data members.  Descriptions are above.
   ************************************************************/
#ifdef MAINTAIN_NN
  double NN2_; // nearest-neighbour
#endif

  char handle_[sizeof(Simplex)];
  // bool inMesh_; // packed with a later field

  ConcentricShells<Uninserted, Point> annuli_;

  size_t ID_:(sizeof(size_t)*8 - 2); // 30 or 62 bits
  bool isSteiner_:1;
  bool inMesh_:1;
  FlagType visited_;
};
#undef dprintf
#undef dputs
#undef MAINTAIN_NN

#endif

---