PointRefiner.h

Go to the documentation of this file.

#ifndef PointRefiner_HEADER
#define PointRefiner_HEADER

/***************************************************************************
 * Benoit Hudson (C) 2007.
 ***************************************************************************/

#include <geometry/Point.h>
#include <delaunay/IncrementalDelaunay.h>
#include <delaunay/DelaunayHelpers.h>
#include <delaunay/PointVertex.h>
#include <utilities/foreach.h>
#include <utilities/BucketPQ.h>
#include <list>
#include <math.h> // sqrt
#include <stack>
#include <stdio.h>

#include <io/PLC.h>
#include <refine-common/RefineConstants.h>
#include <refine-common/Refiner.h>

#include "PointRefineVertex.h"

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

#define COLLECT_STATS

namespace SVR {

/***************************************************************************
 * The refiner itself.
 *
 * Read <i>Hudson, Miller, Phillips 2006</i> to understand what is going on
 * here.
 ***************************************************************************/
template <unsigned ambient>
struct PointRefiner : public MeshRefiner {
  typedef Geometry::Point<ambient> Point;

  // Forward-declare so we can define IncrementalDelaunay.
  struct Vertex;

  typedef ::IncrementalDelaunay<ambient,ambient,Vertex,void, typename PointVertex<ambient>::Printer > IncrementalDelaunay;
  typedef typename IncrementalDelaunay::Simplex Simplex;
  typedef typename IncrementalDelaunay::CenterRadius CenterRadius;
  typedef typename IncrementalDelaunay::Cavity Cavity;
  typedef typename IncrementalDelaunay::Star Star;

  // really just a typedef, but there's a circular dependency with IncrementalDelaunay, and we
  // can't forward-declare typedefs.
  struct Vertex : public PointRefineVertex<ambient, Simplex> {
    public:
    Vertex(const Point& p) : PointRefineVertex<ambient, Simplex>(p) { }
    Vertex(const Point& p, size_t inputid) : PointRefineVertex<ambient, Simplex>(p, inputid) { }
  };

  typedef typename Vertex::FlagType FlagType;

  typedef PointRefineVertex<ambient, Simplex> Uninserted;
  static Vertex *promoteToVertex(Uninserted *u) { return static_cast<Vertex*>(u); }

  /***************************************************************************
   * Initialization.
   *
   * Given points and constants, create the data structures we need to support
   * a call to refine().
   *
   * After this call, the storage for the points in the range can be deleted:
   * I've made a copy.
   ***************************************************************************/
  template <class point_iterator> /* operator* must return a Point<ambient> */
  static PointRefiner *
  initialize(point_iterator begin, const point_iterator& end,
      RefineConstants consts) {
    consts.regularize(ambient, true);
    return new PointRefiner(begin, end, consts);
  }

  static PointRefiner *
  initialize(SVR_IO::PLC *plc, const RefineConstants& consts) {
    return initialize(
        plc->begin_points<ambient>(),
        plc->end_points<ambient>(),
        consts);
  }

  /***************************************************************************
    Refinement.
   ***************************************************************************/
  virtual void refine() {
    doRefine();
  }

  /***************************************************************************
    Printing.
   ***************************************************************************/
  virtual void printToNodeEle(FILE *node, FILE *ele) const {
    doPrint(node, ele);
  }

  /***************************************************************************
    Direct access to the underlying mesh.
   ***************************************************************************/
  IncrementalDelaunay *getDelaunay() { return del_; }
  static const unsigned dimension = ambient;
  virtual unsigned dim() const { return ambient; }



  /***************************************************************************
   ***************************************************************************
   ***************************************************************************
   * From now on, the entire implementation is private.
   ***************************************************************************
   ***************************************************************************
   ***************************************************************************/

  /***************************************************************************
   * \name Trace handling.
   *
   * Useful in generating pretty animations of SVR.
   * Off most of the time.
   ***************************************************************************/
  enum TraceItemType { QUERY, INSERT };
  void printItem(const Simplex& s, const Point& p, TraceItemType t) {
    if (!traceFile_) return;
    fprintf(traceFile_, "* : %c : %s", t==QUERY ? 'T' : 'F', p.toString().c_str());
    BOOST_FOREACH(Vertex *v, s) {
      fprintf(traceFile_, ": %s", v->toPoint().toString().c_str());
    }
    fprintf(traceFile_, ": %g : %s\n", sqrt(center(s).radius2), center(s).center.toString().c_str());
  }
  void printQuery(const Simplex& s, const Point& p) {
    printItem(s, p, QUERY);
  }
  void printInsert(const Cavity& cavity, const Vertex *v) {
    printItem(*cavity.begin(), v->toPoint(), INSERT);
  }


  /***************************************************************************
    Data members.
   ***************************************************************************/
  private:
  IncrementalDelaunay *del_;

  // Work queues, in order.  First, skinny simplices, to get
  // always-good-quality.  Next, crowded Voronoi cells.  Because crowded
  // Steiner cells are so much easier to deal with, we do those first.
  // After that, slivers.
  // The simplex queues are bucketed by size (circumradius^2), and FIFO
  // within a bucket.  This means we quite strongly bias towards smallest-
  // first, but in O(1) per operation. TODO: do the same for vertex queues.
  //typedef stack<Vertex*> VertexStack;
  typedef hudson::FrontStack<Vertex*> VertexStack;
  typedef hudson::BucketDoublePQ<Simplex> SimplexStack;

  SimplexStack skinny_;  
  VertexStack crowdedSteiner_; 
  VertexStack crowdedInput_; 
  SimplexStack sliver_;  
  SimplexStack moderateSkinny_; 

  std::vector<Vertex*> meshvertices_; 
  size_t ninputs_;               
  size_t nuninserted_;         
  size_t nsimplices_;          

  FlagType searchID_; 

  const double rho2_;     // intermediate radius/edge ratio
  const double k2_;       // how far to look when warping
  const bool fastSearch_; // if true, stop the k*r search at the first match
  const double userRho2_; // user-demanded radius/edge ratio
  const bool useOffcenters_;
  const double sigma_; // radius/radius ratio
  const double sliverGrowth_;
  const double perturbFactor_;
  static const size_t SLIVER_RETRIES = 10;
  const RefineConstants::NodeOrder outputNodeOrder_;
  FILE * const traceFile_;


  // stats, largely for debugging
#ifndef COLLECT_STATS
# define statop(op)
# define statopconst(op)
# define statmax(stat, val)
# define statmaxforce(stat, val)
#else
# define statop(op) stats.op
# define statopconst(op) const_cast<struct stats&>(stats).op
# define statmax(stat, val) stats.stat = max(stats.stat, (val))
# define statmaxconst(stat, val) const_cast<struct stats&>(stats).stat = max(stats.stat, (val))
  struct stats {
    unsigned long nrounds_;
    unsigned long npoints_;
    unsigned long ninputs_;
    unsigned long nretries_;
    unsigned long numtri_;
    unsigned long maxtri_;
    unsigned long maxcavity_;
    unsigned long sumcavity_;

    // operations
    struct {
      unsigned long nskinny_;
      unsigned long ncrowdSteiner_;
      unsigned long ncrowdInput_;
      unsigned long nsliver_;
      unsigned long nmedium_;
    } ops;

    // point location
    unsigned long nlocate_;
    unsigned long nsupercav_;
    unsigned long sumsupercav_;
    unsigned long maxsupercav_;
    unsigned long nwarp_;
    double sumwarp_; // distance
    double maxwarp_; // distance

    // intermediate element quality
    unsigned long nslivers_;
    unsigned long nsliverretry_;
    double maxRe_;
    double maxRr_;

    stats(const IncrementalDelaunay *del) {
      memset(this, 0, sizeof(*this));
      npoints_ = Simplicial::countVertices(del->getComplex());
      numtri_ = Simplicial::countSimplices(del->getComplex());
    }

    void updateTri(size_t cavitysize, size_t starsize) {
      npoints_++;
      nrounds_++;
      numtri_ -= cavitysize;
      numtri_ += starsize;
      maxtri_ = max(numtri_, maxtri_);
      maxcavity_ = (cavitysize > maxcavity_) ? cavitysize : maxcavity_;
      sumcavity_ += cavitysize;
    }

    void supercav(size_t scavsize, double mindist2, bool isSteiner) {
      nsupercav_ ++;
      sumsupercav_ += scavsize;
      maxsupercav_ = (maxsupercav_ > scavsize) ? maxsupercav_ : scavsize;

      if (!isSteiner) {
        nwarp_ ++;
        double mindist = sqrt(mindist2);
        sumwarp_ += mindist;
        maxwarp_ = (maxwarp_ > mindist) ? maxwarp_ : mindist;
      }
    }

    void print(bool ignoreSlivers) const {
      printf("---Stats---\n");
      printf("Output: %lu points, %lu simplices\n", npoints_, numtri_);
      printf("Max triangulation size %lu\n", maxtri_);
      printf("Max cavity: %lu; average %g over %lu rounds\n", maxcavity_,
          (double)sumcavity_ / (double)nrounds_, nrounds_);
      printf("Worst radius/edge: %g\n", maxRe_);
      if (!ignoreSlivers) {
        printf("%lu slivers created, %lu retries\n", nslivers_, nsliverretry_);
        printf("Worst sliver: %g\n", maxRr_);
      }

      /* The number of times we ducked the supercavity is the number of times we located
       * but didn't supercavity.
       *
       * The number of times we did any point location at all in split(.) is
       * the number of times we did point location when we didn't retry for slivers. */
      unsigned long no_supercavity = nlocate_ - nsupercav_;
      unsigned long no_supercavity_pct = 100 * no_supercavity / nlocate_;
      printf("Point location: %lu queries (%lu [%lu%%] short-cut)\n",
          nlocate_,
          no_supercavity, no_supercavity_pct);

      printf("Point location: visited %lu simplices (%g per slow search, max %lu)\n",
          sumsupercav_, (double)sumsupercav_ / (double)nsupercav_, maxsupercav_);
      printf("Blindly inserted input %lu times (%lu%%).\n",
          (ninputs_ - nwarp_), (100 * (ninputs_ - nwarp_) / ninputs_));
      printf("Warped %lu times (%lu%%); max warp distance ratio %g, average %g\n",
          nwarp_, (100 * nwarp_ / ninputs_),
          maxwarp_, (double)maxwarp_ / (double)nwarp_);
      printf("Inserted %lu Steiner points\n", (npoints_ - ninputs_));
      printf("Op counts: %lu skinny, %lu crowd(Steiner), %lu crowd(input), %lu sliver, %lu rho\n",
          ops.nskinny_, ops.ncrowdSteiner_, ops.ncrowdInput_, ops.nsliver_, ops.nmedium_);
      printf("-----------\n");
    }
  } stats;
#endif

  /***************************************************************************
   * A mockup of an STL-like set, in order to allow easily switching
   * implementations.
   *
   * We use the visited flag of the Vertex class to help support the search.
   * The assumption is that only one search proceeds at a time.
   * When we run out of IDs, we iterate over all the vertices and clear their
   * flags.  This is exceedingly rare, and is thus stupidly cheap.
   ***************************************************************************/
  struct VertexVisitedSet;
  friend struct VertexVisitedSet;
  struct VertexVisitedSet {
    VertexVisitedSet(PointRefiner *pr): searchID_(++pr->searchID_) {
      /* If we've wrapped the id count, then clear everyone's flag. */
      if (searchID_ == 0) {
        BOOST_FOREACH(Vertex *v, pr->meshvertices_) {
          v->clearVisited();
        }
        searchID_ = ++pr->searchID_;
      }
    }
    ~VertexVisitedSet() {
    }
    struct notpair { bool second; notpair(bool b): second(b) { } };
    notpair insert(Vertex *v) {
      if (v->visitedFlag() == searchID_) { return false; }
      v->setVisited(searchID_);
      return true;
    }
    private:
    FlagType searchID_;
  };


  /***************************************************************************
   * Initialize the data structures.
   ***************************************************************************/
  template <class point_iterator>
  PointRefiner(const point_iterator& begin, const point_iterator& end, const RefineConstants& consts)
  : del_(IncrementalDelaunay::createMeshBox(begin, end, consts.rho)),
    rho2_(consts.rhoprime * consts.rhoprime),
    k2_(consts.k * consts.k),
    fastSearch_(true), //consts.shortcutSearch),
    userRho2_(consts.rho * consts.rho), useOffcenters_ (consts.useOffcenters),
    sigma_(consts.sigma), sliverGrowth_(consts.sliverGrowth), perturbFactor_(consts.perturbFactor),
    outputNodeOrder_(consts.nodeOrder),
    traceFile_(NULL) //consts.traceFile)
#ifdef COLLECT_STATS
    , stats(del_)
#endif
  {
    searchID_ = 0;
    nuninserted_ = (end-begin);
    dprintf("Refining on %lu points\n", (unsigned long)nuninserted_);

#ifndef NDEBUG
    {
      printf("Radius/edge %g (%g intermediate)", sqrt(userRho2_), sqrt(rho2_));
      if (ambient == 2) {
        printf("; off-centers %s\n", useOffcenters_ ? "enabled" : "disabled");
      } else {
        if (sigma_ == 0.0) {
          printf("; slivers ignored\n");
        } else {
          printf("; sliver quality %g\n", sigma_);
        }
      }
    }
#endif

    initialize(begin, end);
  }

  struct CompareVertices {
    bool operator()(const Vertex *a, const Vertex *b) const {
      return a->id() < b->id();
    }
  };

  void doRefine() {
    // This is the entire SVR algorithm, except for the details.
    dprintf("Running SVR\n");
    while(!queueEmpty()) {
      performQueueOperation();
    }
    statop(print(sigma_ == 0.0));
    assert(isClean());
    assert(nuninserted_ == 0);

    switch(outputNodeOrder_) {
      case RefineConstants::ID_ORDER:
        sort(meshvertices_.begin(), meshvertices_.end(), CompareVertices());
        break;
      case RefineConstants::INSERTION_ORDER_RENUMBER:
        for(size_t i = 0; i < meshvertices_.size(); ++i) {
          meshvertices_[i]->setID(i);
        }
        break;
      case RefineConstants::INSERTION_ORDER_NORENUMBER:
        break;
    }
  }

  /***************************************************************************
   * Print the mesh.
   ***************************************************************************/
  void doPrint(FILE *node, FILE *ele) const {
    /* print the .node file */
    fprintf(node, "%lu %u 0 0\n", (unsigned long)meshvertices_.size(), ambient);
    BOOST_FOREACH(Vertex *v, meshvertices_) {
      Delaunay::printPoint(v->id(), v->toPoint(), node);
    }

    /* print the .ele file */
    fprintf(ele, "%lu %u 0\n", (unsigned long)nsimplices_, ambient+1);
    DoPrintSimplex dps(ele);
    del_->getComplex().dfsBySimplex(dps, del_->getHandle());
  }
  struct DoPrintSimplex {
    size_t n_;
    FILE *out_;
    DoPrintSimplex(FILE *out): n_(0), out_(out) { }
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& s) {
      boost::array<size_t, ambient+1> arr;
      for (unsigned i = 0; i <= ambient; ++i) {
        arr[i] = s[i]->id();
      }
      Delaunay::printEle(n_++, arr, out_);
      return false;
    }
  };

  /***************************************************************************
   * Initialize the point-location tables and the work queues.
   ***************************************************************************/
  template <class point_iterator>
  void initialize(point_iterator begin, point_iterator end) {
    vector<Simplex> simps;
    Simplicial::collectSimplices(del_->getComplex(), simps);
    nsimplices_ = simps.size();

    dprintf("Starting with %u simplices\n", nsimplices_);

    // Find the mesh vertices, and create handles for each.
    VertexVisitedSet visited(this);
    BOOST_FOREACH(const Simplex& s, simps) {
      BOOST_FOREACH(Vertex *v, s) {
        if (!visited.insert(v).second) { continue; }
        initVertex(v, s);
      }
    }
    dprintf("Starting with %u vertices\n", (unsigned)meshvertices_.size());

    // Create the uninserted vertices and assign them to the mesh
    // Voronoi cells (i.e. the nearest mesh vertex).
    size_t id = 0;
    BOUNDED_FOREACH(const Point& p, begin, end) {
      Uninserted *u = new Uninserted(p, id++);
      Vertex *best = bestPoint<Vertex*>(meshvertices_.begin(), meshvertices_.end(), p).first;
      assert(best);
      best->assign(u);
    }
    statop(ninputs_ = id);

    // Queue up the work: skinny tets, crowded vertices.
    dprintf("Filling queues\n");
    addToQueues(simps.begin(), simps.end());
    BOOST_FOREACH(Vertex *v, meshvertices_) {
      addToQueues(v);
    }
  }

  /***************************************************************************
   * Maintain the points table.
   * Now just a call to the vertices.
   ***************************************************************************/
  void redistributePoints(Vertex *oldv, Vertex *newv) {
    oldv->reassign(newv);
  }

  /***************************************************************************
   * \name Vertex->simplex mapping.
   ***************************************************************************/
  void initVertex(Vertex *v, const Simplex& s) {
    v->initHandle(s);
    if (v->isSteiner()) {
      // Steiner point 0 is numbered after all the input vertices.
      // First, compute which Steiner point this is.  Then add ninputs_.
      size_t ninserted = ninputs_ - nuninserted_;
      size_t nsteiner = meshvertices_.size() - ninserted;
      size_t id = nsteiner + ninputs_;
      v->setID(id);
    }
    meshvertices_.push_back(v);
  }
  void setHandle(Vertex *v, const Simplex& s) {
    v->setHandle(s);
  }
  const Simplex& getHandle(Vertex *v) const {
    return v->getHandle();
  }

  /***************************************************************************
   * \name Splitting.
   *
   * Split a simplex.  Either use the candidate point, or try to use
   * the circumcenter (perhaps perturbed to avoid slivers), offcenter, or a
   * nearby uninserted input point.
   *
   * Warning: the candidate point, if a Steiner, is now owned by this routine.
   * In particular, if it causes a small sliver, I delete the candidate and
   * create a new one.
   ***************************************************************************/
  void split(const Simplex& s, Uninserted *candidate = NULL) {
    dprintf("splitting %s\n", s.toString().c_str());
    if (!isMember(s)) {
      // already removed.
      dprintf("  moot\n");
      return;
    }

    // For slivers, we want to try several times.  This requires:
    // (a) stopping after some number of tries.
    // (b) checking if new slivers are large, and thus checking our radius.
    const double r2_for_sliver = center(s).radius2;
    size_t ntries = 0;

    Cavity cavity;
    Star star;
    Vertex *v;

makecavity:

    /*
     * Pick a vertex to insert.  This will either be the one passed in
     * (if we know that it's safe to insert), or one found by splitPoint(.)
     * Both of those give us an Uninserted point; promote it to a Vertex.
     */
    if (candidate && ntries == 0) {
      v = promoteToVertex(candidate);
    } else {
      dprintf("  attempt %u\n", (unsigned)ntries);
      v = promoteToVertex(splitPoint(s, ntries > 0));
    }

    del_->computeCavity(s, v->toPoint(), cavity);
    del_->computeStar(cavity, v, star);
    dprintf("  cavity %u, star %u, point %s\n", (unsigned)cavity.size(),
        (unsigned)star.size(), v->toString().c_str());

    ntries++;

    // If we care about slivers, check for new small slivers
    // and run that again if we found one.
    if (nuninserted_ == 0
        && v->isSteiner()
        && ntries < SLIVER_RETRIES
        && hasSmallSliver(r2_for_sliver, star)) {
      cavity.clear();
      star.clear();
      delete v; // v is definitely a Steiner: we just checked.
      goto makecavity;
    }

    statop(nsliverretry_ += (ntries-1));

    // Finally, commit.
    commitStar(cavity, star, v);

    // update the queues.
    addToQueues(star.begin(), star.end());
    addToQueues(v);
  }

  Uninserted * splitPoint(const Simplex& s, bool perturb) {
    const CenterRadius& cr = center(s);
    const Point& c = cr.center;
    double r2 = cr.radius2;
    dprintf("  %s, radius %g\n", c.toString().c_str(), sqrt(r2));
    statopconst(nlocate_++);

    printQuery(s, c);

    // Try to avoid computing the supercavity, which is very expensive.
    if (perturb) {
      Point p = c + Point::randPoint(perturbFactor_ * sqrt(r2));
      dprintf("  using Steiner perturbed to %s\n", p.toString().c_str());
      return new Uninserted(p);
    }
    if (nuninserted_ == 0) {
      Point offc = offcenter(s);
      dprintf("  no uninserted; using off-center\n");
      return new Uninserted(offc);
    }

    // Compute the supercavity
    pair<Uninserted*, double> w = findWarpPoint(s, cr);
    Uninserted *best = w.first;

    if (best) {
      dprintf("  warping to %s (dist %g, ratio %g)\n",
          best->toString().c_str(), sqrt(w.second), sqrt(w.second / r2));
      return best;
    } else {
      Point offc = offcenter(s);
      dprintf("  using off-center %s\n", offc.toString().c_str());
      return new Uninserted(offc);
    }
  }

  static double dist2(const Point& p, const Point& q) {
    return p.dist2(q);
  }
  template <class V>
  static double dist2(const Point& p, const V *v) {
    return v->dist2(p);
  }
  template <class V>
  static double dist2(const V *v, const Point& p) {
    return v->dist2(p);
  }
  template <class V, class U>
  static double dist2(const V *v, const U *u) {
    return u->dist2(v);
  }

  template <class iteratorV, class V>
  static void bestPoint(iteratorV begin, iteratorV end, const Point& c, V& best, double& mindist2) {
    while (begin != end) {
      double ambient = dist2(c, *begin);
      if (ambient < mindist2) {
        best = *begin;
        mindist2 = ambient;
      }
      ++begin;
    }
  }
  template <class V, class iteratorV>
  static pair<V, double> bestPoint(iteratorV begin, iteratorV end, const Point& c) {
    assert(begin != end);
    pair<V, double> p(*begin, dist2(c, *begin));
    ++begin;
    bestPoint(begin, end, c, p.first, p.second);
    return p;
  }

  /***************************************************************************
   * Search the supercavity for nearby points.  Note that any vertex whose
   * Voronoi cell intersects the ball has at least one circumcenter that
   * intersects the ball.
   *
   * We start with a search radius of k * r.
   * When we find a point, we shrink the radius: anything we'd see further is
   * no longer useful.
   * When we find a very close point, we cancel the search.
   * Of course, we only search each mesh vertex once.
   ***************************************************************************/
  struct SCavStruct {
    PointRefiner& refine_;
    size_t nchecked_;
    Point c_;
    double threshold2_; // what's good enough? (this is static; may be infinite)
    double searchR2_;   // how far to search? (this shrinks)
    Uninserted *best_;
    VertexVisitedSet visited_; // set of vertices we've seen

    SCavStruct(PointRefiner& r, const CenterRadius& cr)
      : refine_(r), visited_(&r) {
        c_ = cr.center;
        searchR2_ = refine_.k2_ * cr.radius2; // this is the distance to beat
        threshold2_ = refine_.fastSearch_ ? std::numeric_limits<double>::max() : (searchR2_ / 4.0); // good enough in kr/2
        best_ = NULL;
        nchecked_ = 0;
      }

    // We can enter if the ball of s intersects the ball defined by (c,r)
    // The balls intersect if |pq| < r + r'
    // Equivalently, check |pq|^2 < r^2 + r'^2 + 2.0 sqrt(r2r'2)
    // It might be worth doing a branch to avoid the sqrt, but probably not.
    bool canEnter(const Simplex& s) const {
      const CenterRadius& cr = refine_.center(s);
      double pq2 = c_.dist2(cr.center);
      double r2 = cr.radius2;
      double rprime2 = searchR2_;
      return pq2 < r2 + rprime2 + 2.0 * sqrt(r2 * rprime2);
    }

    bool choose(const Simplex& s) {
      nchecked_++;
      BOOST_FOREACH(Vertex *v, s) {
        if (!visited_.insert(v).second) { continue; } // already checked

        pair<Uninserted*, double> p;
        if (threshold2_ == std::numeric_limits<double>::max()) {
          p = v->findPoint(c_, searchR2_);
        } else {
          p = v->findClosest(c_, searchR2_);
        }
        if (p.first) {
          best_ = p.first;
          searchR2_ = p.second;
          if (searchR2_ < threshold2_) {
            return true;
          }
        }
      }
      return false;
    }
  };

  pair<Uninserted*, double> findWarpPoint(const Simplex& s, const CenterRadius& cr) {
    SCavStruct finder(*this, cr);
    del_->getComplex().bfsBySimplex(finder, s);
    statopconst(supercav(finder.nchecked_, finder.searchR2_ / cr.radius2,
          finder.best_ == NULL));
    return make_pair(finder.best_, finder.searchR2_);
  }

  /***************************************************************************
   * Commit the star.
   * Update all the necessary data structures.
   ***************************************************************************/
  void commitStar(const Cavity& cavity, Star& star, Vertex *newv) {
    printInsert(cavity, newv);

    // Update the mesh.
    del_->insertVertex(cavity, star);
    initVertex(newv, *star.begin());
    nsimplices_ += (star.size() - cavity.size());

    if (newv->isSteiner()) {
      // The mesh now owns this vertex, not us.
      del_->giveVertexOwnership(newv);
    } else {
      nuninserted_--;
      if (nuninserted_ == 0) {
#ifndef NDEBUG
        printf("All input vertices in the mesh; removing slivers, improving quality.\n");
#endif
        //crowded_.clear();
        while (!crowdedSteiner_.empty()) crowdedSteiner_.pop();
        while (!crowdedInput_.empty()) crowdedInput_.pop();
      }
    }

    // Redistribute the points.
    // Note that we need not redistribute points if
    // there are none to be distributed.  That keeps
    // some handles alive longer than necessary.
    // Presumably this is no big deal (presumably cheaper
    // than star.size() * (d+1) iterations every time).
    if (nuninserted_ > 0) {
      VertexVisitedSet vset(this);
      /* Ideally we'd iterate over the cavity, but we need
         to set new handles.  So we iterate over the star and
         avoid touching the new vertex. */
      BOOST_FOREACH(const Simplex& s, star) {
        BOOST_FOREACH(Vertex *oldv, s) {
          if (oldv == newv) { continue; }
          if (!vset.insert(oldv).second) { continue; }
          redistributePoints(oldv, newv);
          setHandle(oldv, s);
        }
      }
    }

    statop(updateTri(cavity.size(), star.size()));
  }

  /***************************************************************************
   * \name Geometric computations and predicates
   ***************************************************************************/
  const CenterRadius& center(const Simplex& s) const {
    return del_->circumcenter(s);
  }
  Point offcenter(const Simplex& s) const {
    if (useOffcenters_) {
      return del_->offCenter(s, userRho2_);
    } else {
      return center(s).center;
    }
  }

  bool inSphere(const Simplex& s, const Point& p) const {
    return del_->inSphereAffine(s, p);
  }
  bool inSimplex(const Simplex& s, const Point& p) const {
    return del_->inSimplexAffine(s, p);
  }
  bool inSphere(const Simplex& s, Vertex *v) const {
    return del_->inSphereAffine(s, v->toPoint());
  }
  bool inSimplex(const Simplex& s, Vertex *v) const {
    return del_->inSimplexAffine(s, v->toPoint());
  }

  bool isMember(const Simplex& s) const {
    return del_->getComplex().isMember(s);
  }


  enum Skinnyness {
    FAT = 0,
    MEDIUM,
    SKINNY
  };

  Skinnyness computeSkinniness(const Simplex& s) const {
    double re2 = del_->radiusEdge2(s);
    if (re2 > rho2_) {
      dprintf("  R/e = %g ==> skinny: %s\n", sqrt(re2), s.toString().c_str());
      return SKINNY;
    } else if (re2 > userRho2_) {
      dprintf("  R/e = %g ==> moderate skinny: %s\n", sqrt(re2), s.toString().c_str());
      return MEDIUM;
    } else {
      return FAT;
    }
  }

  bool isCrowded(const Vertex *v) const {
    if (nuninserted_ == 0) {
      return false;
    } else {
      return v->isCrowded();
    }
  }

  bool isSliver(const Simplex& s) const {
    if (sigma_ == 0.0) return false;
    if (computeSkinniness(s) != FAT) {
      return false;
    }
    double Rr = del_->radiusRadius(s);
    if (Rr > sigma_) {
      dprintf("  R/r = %g ==> sliver: %s\n", Rr, s.toString().c_str());
    }
    return Rr > sigma_; // TODO!!
  }

  bool hasSmallSliver(double r2, const Star& star) const {
    if (sigma_ == 0.0) return false;
    BOOST_FOREACH(const Simplex& s, star) {
      if (center(s).radius2 < sliverGrowth_ * r2) {
        if (isSliver(s)) {
          return true;
        }
      }
    }
    return false;
  }

  // for debugging at the end, to make sure we didn't miss something
  struct FindBad {
    const PointRefiner& refine_;
    FindBad( const PointRefiner& r  ): refine_(r) { }
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& s) {
      Skinnyness sness = refine_.computeSkinniness(s);
      bool sliver = refine_.isSliver(s);
      if (sness != FAT) {
        printf("skinny: %s\n", s.toString().c_str());
      }
      if (sliver) {
        printf("sliver: %s\n", s.toString().c_str());
      }
      bool crowd = false;
      BOOST_FOREACH(Vertex *v, s) {
        if (refine_.isCrowded(v)) {
          crowd = true;
          printf("crowded: %s\n", v->toString().c_str());
        }
      }
      return (sness!=FAT) || sliver || crowd;
    }
  };

  bool isClean() const {
    FindBad fb(*this);
    return !del_->getComplex().dfsBySimplex(fb, del_->getHandle());
  }

  /***************************************************************************
   * \name Queue handling
   ***************************************************************************/
  bool queueEmpty() const {
    return skinny_.empty()
      && crowdedSteiner_.empty()
      && crowdedInput_.empty()
      && sliver_.empty()
      && moderateSkinny_.empty();
  }

  struct FindCrowdedSimplex {
    Vertex *v_;
    Uninserted *u_;
    IncrementalDelaunay *del_;
    FindCrowdedSimplex (Vertex *v, Uninserted *u, IncrementalDelaunay *del): v_(v), u_(u), del_(del) { }
    bool canEnter(const Simplex& s) const { return s.has(v_); }
    bool choose(const Simplex& s) const { return del_->inSphereAffine(s, u_->toPoint()); }
  };

  void performQueueOperation() {
    if (!skinny_.empty()) {
      Simplex s = skinny_.top();
      skinny_.pop();
      dprintf("popped skinny: %s\n", s.toString().c_str());
      assert(isMember(s));
      statop(ops.nskinny_++);
      split(s);
    }
    else if (!crowdedSteiner_.empty()) {
      Vertex *v = crowdedSteiner_.top();
      // do NOT pop -- we're likely to still be crowded!
      dprintf("popped crowded: %s\n", v->toString().c_str());
      assert(isCrowded(v));
      assert(v->isSteiner());
      statop(ops.ncrowdSteiner_++);

      /* Pick the farthest point from v */
      Uninserted *u = v->getFarPoint();

      /* Find a simplex around v that has u in its circumball. */
      FindCrowdedSimplex fcs(v, u, del_);
      Simplex s = *del_->getComplex().findSimplex(fcs, getHandle(v));

      /* We know already that u is safe to insert -- else we'd have warped to u */
      split(s, u);
    }
    else if (!crowdedInput_.empty()) {
      Vertex *v = crowdedInput_.top();
      // do NOT pop -- we're likely to still be crowded!
      dprintf("popped crowded: %s\n", v->toString().c_str());
      assert(isCrowded(v));
      assert(!v->isSteiner());
      statop(ops.ncrowdInput_++);

      /* Pick the farthest point from v */
      Uninserted *u = v->getFarPoint();

      /* Find a simplex around v that has u in its circumball. */
      FindCrowdedSimplex fcs(v, u, del_);
      Simplex s = *del_->getComplex().findSimplex(fcs, getHandle(v));

      /* u may not be safe to insert, so find an appropriate warp point. */
      split(s);
    }
    else if (!sliver_.empty()) {
      Simplex s = sliver_.top();
      sliver_.pop();
      dprintf("popped sliver: %s\n", s.toString().c_str());
      assert(isMember(s));
      statop(ops.nsliver_++);
      split(s);
    }
    else if (!moderateSkinny_.empty()) {
      Simplex s = moderateSkinny_.top();
      moderateSkinny_.pop();
      dprintf("popped moderate skinny: %s\n", s.toString().c_str());
      assert(isMember(s));
      statop(ops.nmedium_++);

      /* We know there are no uninserted points left, because moderate-skinnies are
         processed after all crowding.  Therefore, just blindly insert the offcenter
         (blindly except for slivers). */
      assert(nuninserted_ == 0);

      Uninserted *offc = new Uninserted(offcenter(s));
      split(s, offc);
    }
    else {
      assert(0 && "queue is empty");
    }
  }

  template <class iteratorS>
  void addToQueues(iteratorS begin, iteratorS end) {
    // Clean up the queues.  This allows the complex to reap simplices,
    // and it allows our queues to do so also.
    while(!skinny_.empty() && !isMember(skinny_.top())) { skinny_.pop(); }
    while(!sliver_.empty() && !isMember(sliver_.top())) { sliver_.pop(); }
    while(!moderateSkinny_.empty() && !isMember(moderateSkinny_.top())) { moderateSkinny_.pop(); }

    BOUNDED_FOREACH(const Simplex& s, begin, end) {
      Skinnyness sness = computeSkinniness(s);
      if (sness == SKINNY) {
        dprintf("push skinny: %s\n", s.toString().c_str());
        skinny_.push(center(s).radius2, s);
        statmax(maxRe_, sqrt(del_->radiusEdge2(s)));
      }
      // else is not strictly necessary, but avoids
      // duplicate insertions.
      else if (isSliver(s)) {
        dprintf("push sliver: %s\n", s.toString().c_str());
        sliver_.push(center(s).radius2, s);
        statop(nslivers_++);
        statmax(maxRr_, del_->radiusRadius(s));
      }
      else if (sness == MEDIUM) {
        dprintf("push moderate skinny: %s\n", s.toString().c_str());
        moderateSkinny_.push(center(s).radius2, s);
        statmax(maxRe_, sqrt(del_->radiusEdge2(s)));
      }
      else assert(sness == FAT);
    }
  }

  void addToQueues(Vertex *v) {
    while(!crowdedSteiner_.empty() && !isCrowded(crowdedSteiner_.top())) { crowdedSteiner_.pop(); }
    while(!crowdedInput_.empty() && !isCrowded(crowdedInput_.top())) { crowdedInput_.pop(); }
    if (isCrowded(v)) {
      if (v->isSteiner()) {
        crowdedSteiner_.push(v);
      } else {
        crowdedInput_.push(v);
      }
    }
  }

  /***************************************************************************
   * FIN.
   ***************************************************************************/
#undef dprintf
}; // end class
} // end namespace

#endif

---