Mesh.h

Go to the documentation of this file.

#ifndef Mesh_HEADER
#define Mesh_HEADER

#include "GenericMesh.h"
#include "ProtectedBall.h"
#include "RefineVertex.h"
#include "SimplexData.h"
#include "SliverChecker.h"
#include <algorithm>
#include <boost/foreach.hpp>
#include <boost/function_output_iterator.hpp>
#include <delaunay/IncrementalDelaunay.h>
#include <functional>
#include <stdlib.h>
#include <utilities/hash.h>
#include <vector>

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

/* signature MESH =
   sig
     type vertex
     type simplex
     type point
     val center : simplex -> (point * real)
     val split : mesh * simplex -> { v: vertex, kill: simplex list, born: simplex list }
     val encroached : mesh * simplex * point -> bool
     val notify : simplex * { kill: ball list, born: ball list, v: vertex } -> unit
     val getHandle : mesh -> simplex
   end
   */

template <size_t ambient, size_t topo>
struct Mesh : public GenericMesh<ambient> {
  typedef ::SimplexData<ambient, topo> SimplexData;
  typedef IncrementalDelaunay<ambient, topo, RefineVertex<ambient>, SimplexData> Delaunay;
  typedef typename Delaunay::Complex Complex;
  typedef typename Delaunay::Cavity Cavity;
  typedef typename Delaunay::Star Star;
  typedef typename Delaunay::Simplex Simplex;
  typedef typename Delaunay::Vertex Vertex;
  typedef typename MeshTypes<ambient>::Ball Ball;
  typedef typename MeshTypes<ambient>::BallSet BallSet;
  typedef typename MeshTypes<ambient>::GenericMesh GenericMesh;
  typedef Geometry::CenterRadius<ambient> CenterRadius;
  typedef Geometry::Point<ambient> Point;
  typedef Geometry::Point<topo> PPoint;
  typedef ::SplitData<ambient> SplitData;

  static const size_t ambient_dimension = ambient;
  static const size_t topological_dimension = topo;

  private:
  struct InitVertex {
    Mesh *mesh;
    InitVertex(Mesh *mesh): mesh(mesh) { }
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& s) {
      BOOST_FOREACH(Vertex *v, s) {
        if (topo == ambient) {
          v->replaceHandle(mesh->toBall(s));
          // the 'upper' mesh is NULL, to save a word per vertex.
        } else {
          v->addUpperMesh(mesh);
        }
      }
      return false;
    }
  };


  public:
  Mesh(Delaunay *del) : del_(del) {
    dprintf("%u-mesh %p: Projection plane %s\n", unsigned(topo), this,
        del_->getPlane().toString().c_str());
    assert(del_->isOriented());
    InitVertex sh(this);
    dfs(sh);
  }

  virtual unsigned topological() const { return topo; }

  const Delaunay *getDelaunay() const {
    return del_;
  }

  const Complex& getSimplicialComplex() const {
    return del_->getComplex();
  }

  Delaunay *getDelaunayDangerously() {
    return del_;
  }

  const Simplex& getHandle() const {
    return del_->getHandle();
  }

  template <class DFSData>
  void dfs(DFSData& data) const {
    del_->getComplex().dfsBySimplex(data, getHandle());
  }


  /***************************************************************************
   ***************************************************************************/
  private:
  template <class output_iterator>
  struct ConvertToBalls {
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& s) {
      *out_ = mesh_->toBall(s);
      ++out_;
      return false;
    }
    const Mesh *mesh_;
    output_iterator out_;
    ConvertToBalls(const Mesh *m, output_iterator o): mesh_(m), out_(o) { }
  };

  public:
  template <class output_iterator>
  output_iterator collectAllBalls(output_iterator out) const {
    ConvertToBalls<output_iterator> ctb(this, out);
    dfs(ctb);
    return ctb.out_;
  }

  /***************************************************************************
   * Convert a simplex into its protected ball.
   * Each simplex has a unique ball, so this is const-like.  If the ball has
   * not yet been created, create it here (violating constness).
   ***************************************************************************/
  Ball *toBall(const Simplex& s) const {
    Ball *b = getData(s).toBall();
    if (b == NULL) {
      // Create a new ball.  The ball must know the mesh as non-const,
      // and it calls setBall below upon creation.  Thus I violate
      // constness.
      b = makeProtectedBall(const_cast<Mesh*>(this), s);
    }
    return b;
  }

  void setBall(const Simplex& s, Ball *b) {
    assert(getData(s).toBall() == NULL);
    getDataRW(s).setBall(b);
  }

  /***************************************************************************
   * \name Geometric and combinatorial information about a simplex.
   ***************************************************************************/

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

  CenterRadius circumcenter(const Simplex& s) const {
    return del_->circumcenter(s);
  }

  double radiusEdge2(const Simplex& s) const {
    return del_->radiusEdge2(s);
  }
  double computeSigma(const Simplex& s) const {
    return del_->computeSigma(s);
  }

  bool encroachedBy(const Simplex& s, const Vertex *v) const {
    // Check if v is a vertex of s.
    if (s.has(v)) {
      dprintf("  vertex %p does not encroach %s\n", v, toBall(s)->toString().c_str());
      return false;
    }

    // Otherwise, precision isn't necessary here, so use the cheap computation.
    CenterRadius cr = circumcenter(s);
    bool b = v->toPoint().dist2(cr.center) < cr.radius2;
    dprintf("  vertex %p (%s) %s <%s,%g>  (%s)\n",
        v, v->toPoint().toString().c_str(),
        b ? "encroaches" : "does not encroach",
        cr.center.toString().c_str(), sqrt(cr.radius2),
        toBall(s)->toString().c_str());
    return b;
  }

  bool inSphereExact(const Simplex& s, const Vertex *v) const {
    assert(v->hasUpperMesh(this));
    if (s.has(v)) { return false; }
    return del_->inSphereAffine(s, v->toPoint());
  }

  bool isResolved(const Simplex& s) const {
    if (topological_dimension == ambient) {
      // obviously, since this ball is *in* the top-d mesh.
      return true;
    }

    // Look at each top-d simplex that we intersect.
    BOUNDED_FOREACH(Ball *b, getData(s).begin_uppers(), getData(s).end_uppers()) {
      if (b->topological() != ambient) { continue; }

      // Look for all the vertices.  If we find them, we're resolved by this simplex.
      bool hasAll = true;
      BOOST_FOREACH(Vertex *v, s) {
        if (!b->has(v)) { hasAll = false; break; }
      }
      if (hasAll) { return true; }
    }

    // if we made it here, the ball is unresolved
    return false;
  }


  void gatherUppers(const Simplex& s, BallSet& set) const {
    return getData(s).gatherUppers(set);
  }
  void gatherLowers(const Simplex& s, BallSet& set) const {
    return getData(s).gatherLowers(set);
  }

  bool approximatelyIntersects(const Simplex& s, const Ball *b) const {
    return b->approximatelyIntersects(circumcenter(s));
  }

  /***************************************************************************
    \name Modification routines
   ***************************************************************************/

  private:
  bool isDegenerate(const Simplex& s) const {
    if (topo == 1) {
      // we never split vertices, so this can't be a problem.
      return false;
    }
    vector<const GenericMesh*> meshes;
    std::copy(s[0]->begin_uppers(topo - 1), s[0]->end_uppers(topo - 1),
        std::back_inserter(meshes));
    // std::sort(meshes.begin(), meshes.end());

    for(size_t i = 1 ; i <= topo; ++i) {
      // vector<const GenericMesh*> mis;
      // std::copy(s[i]->begin_uppers(topo - 1), s[i]->end_uppers(topo - 1),
      //     std::back_inserter(mis));
      // std::sort(mis.begin(), mis.end());

      // Compute the set intersection.  Sadly, not in place, but at least we
      // use the swap to avoid a redundent copy.  By using the assumption that
      // the uppers of s[i] are sorted, we at least eliminate one copy.
      vector<const GenericMesh*> tmp;
      std::set_intersection(meshes.begin(), meshes.end(),
          s[i]->begin_uppers(topo - 1), s[i]->end_uppers(topo - 1),
          std::back_inserter(tmp));
      meshes.swap(tmp);
    }

    // Degeneracy only occurs if all the meshes lie on a common (topo-1)-feature.
    // We therefore return whether the list of such is non-empty.
    assert(meshes.size() == 0 || meshes.size() == 1);
    return !meshes.empty();
  }

  private:
  struct WarpClosure {
    Vertex *best_;
    double radius2_;
    const Point& center_;
    Mesh *mesh_;

    WarpClosure(double radius2, const Point& center, Mesh *mesh)
      : radius2_(radius2), center_(center), mesh_(mesh) {
        best_ = NULL;
      }

    bool canEnter(const Simplex& s) {
      // do we need this precise test?  Seems unlikely...
      return mesh_->getDelaunay()->inSphereAffine(s, center_);
    }

    bool choose(const Simplex& s) {
      BOOST_FOREACH(Vertex *v, s) {
        pair<Vertex*,double> p = v->findClosest(mesh_, center_, radius2_);
        if (p.first) {
          best_ = p.first;
          radius2_ = p.second;
        }
      }
      return false;
    }
  };

  Vertex *findWarpVertex(const Simplex& s, double k2) {
    WarpClosure wc(circumcenter(s).radius2 * k2, circumcenter(s).center, this);
    del_->getComplex().bfsBySimplex(wc, s);
    return wc.best_;
  }

  public:
  bool split(const Simplex& s, SplitData& data, Vertex *inserthint = NULL, bool force = false) {
    dprintf("Splitting in %u-mesh [%p] at %s\n", unsigned(topo), this, s.toString().c_str());
    if (inserthint) { dprintf("  Using insert-hint of %s\n", inserthint->toString().c_str()); }

    /* Find a vertex to warp to. */
    Vertex *toinsert = inserthint ? inserthint : findWarpVertex(s, data.k * data.k);
    bool isSteiner = toinsert == NULL;

    // Create the new vertex if needed.  Remember to delete it if we don't
    // insert it.
    if (isSteiner) {
      /*
       * If we're in ambient dimension 3 or higher, check if we're about to
       * create small slivers in higher dimension.  If so, perturb the point
       * until we don't.
       */
      Point newpoint = circumcenter(s).center;
      if (ambient >= 3 && data.sigma != 0.0) {
        SliverChecker<ambient> checker(getData(s).begin_uppers(),
            getData(s).end_uppers(), data, circumcenter(s).radius2);
        enum { MAX_SLIVER_ITERATIONS = 20 };
        unsigned ntries = 0;
        while(++ntries < MAX_SLIVER_ITERATIONS && checker.formsSliver(newpoint)) {
          /* Choose a better newpoint.  This uses the fact that projection
           * preserves distances. */
          PPoint c = del_->getPlane().project(circumcenter(s).center);
          c += PPoint::randPoint(data.perturbFactor);
          newpoint = del_->getPlane().unproject(c);
        }
        if (ntries == MAX_SLIVER_ITERATIONS) {
          printf("warning: can't avoid sliver after %u tries\n",
              MAX_SLIVER_ITERATIONS);
        }
      }
      toinsert = new Vertex(newpoint);
      toinsert->setContainingDimension(topo);
      if (topo < ambient) { toinsert->addUpperMesh(this); }
    }
    dprintf("  split at %s (%s) %s\n", toinsert->toString().c_str(),
        toinsert->toPoint().toString().c_str(),
        isSteiner ? "(Steiner)" : "(input)");

    /* Compute the cavity.  If we end up inserting, we will need the
     * neighbours: the union of the circumballs of cavity and neighbours covers
     * the star.  If we don't, the cost to have collected them is miniscule. */
    Cavity cavity;
    vector<Simplex> neighbours;
    del_->computeCavity(s, toinsert->toPoint(), cavity, neighbours);

    /* Check for encroachment.  If anything is encroached, it must be in the cavity.
     * TODO -- this adds a log(L/s) factor.  We should allow
     * single-encroachment for lower-containing-dimension vertices.
     */
    {
      hudson::OrderedHashSet<Ball *, hudson::HashPointer<Ball> > checked;
      BOOST_FOREACH(const Simplex& s, cavity) {
        dprintf("  checking for encroachment in %s\n", toBall(s)->toString().c_str());
        BOUNDED_FOREACH(Ball *b, getData(s).begin_lowers(), getData(s).end_lowers()) {
          if (checked.insert(b).second) { // then we've not seen it before
            if (b->encroachedBy(toinsert)) {
              data.insertEncroached(b);
            }
          }
        }
      }
    }

    // If we encroached, back out (unless forced).  The caller is responsible
    // for putting the current job back on the queue.
    if (data.encroached()) {
      if (force) {
        dprintf("  encroached; force anyway.\n");
      } else {
        dprintf("  encroached; back out and try again\n");
        if (isSteiner) { delete toinsert; }
        return false;
      }
    }

    dprintf("  starting ForceInsert\n");

    /*
     * Compute the star and actually insert the vertex.
     * In the top dimension, we have a bounding box.  In lower dimensions,
     * however, the star may include degenerate simplices where the point toinsert
     * is exactly on the boundary (because it is the result of splitting a
     * lower-dimensional protected ball).  We remove those from the star.
     */
    Star star;
    del_->computeStar(cavity, toinsert, star);
    if (topo != ambient && toinsert->getCD() < topo) {
      typename Star::iterator it = star.begin();
      while (it != star.end() /* recompute: we're modifying star */) {
        const Simplex& s = *it;
        if (isDegenerate(s)) {
          dprintf("   reject degenerate %s\n", s.toString().c_str());
          it = star.erase(it);
        } else {
          dprintf("   accept non-degenerate %s\n", s.toString().c_str());
          ++it;
        }
      }
    }
#ifndef NDEBUG
    BOOST_FOREACH(const Simplex& s, star) {
      dprintf("   verifying %s\n", s.toString().c_str());
      // This may occur in d=3, and isn't a problem: we do have
      // slivers of 4 points "on" a plane due to numerical error.
      // assert(!isDegenerate(s));

      assert(del_->isOriented(s));
    }
#endif

    dprintf("  star computed\n");
    del_->insertVertex(cavity, star);

    dprintf("  starting to reassign uppers/lowers\n");

    /* Link in the new star to its uppers and lowers. */
    BallSet lowers;
    BallSet uppers;
    typedef typename Vertex::VertexUniqueList VertexSet;
    typedef hudson::FastHashMap<const GenericMesh*, VertexSet, hudson::HashPointer<GenericMesh> >
      MeshVertexMap;
    MeshVertexMap upvertices;
    std::vector<const GenericMesh*> upmeshes;

    BOOST_FOREACH(const Simplex& s, cavity) {
      const SimplexData& sdata = getData(s);
      lowers.insert(sdata.begin_lowers(), sdata.end_lowers());
      uppers.insert(sdata.begin_uppers(), sdata.end_uppers());
    }

    BOOST_FOREACH(const Simplex& s, neighbours) {
      const SimplexData& sdata = getData(s);
      lowers.insert(sdata.begin_lowers(), sdata.end_lowers());
      uppers.insert(sdata.begin_uppers(), sdata.end_uppers());
    }

    // Gather the upper vertices, ordered by mesh.
    BOOST_FOREACH(const Ball *b, uppers) {
      const GenericMesh *m = b->getMesh();
      // find or insert the mesh into upvertices, then gather this
      // ball's vertices into the set.
      std::pair<typename MeshVertexMap::iterator, bool> p =
        upvertices.insert(m);
      if (p.second) { upmeshes.push_back(m); }
      VertexSet& vset = ((*p.first).second);
      b->gatherVertices(vset);
    }

    // Set uppers and lowers for the new star, one by one.
    BOOST_FOREACH(const Simplex& s, star) {
      BOOST_FOREACH(Ball *b, lowers) {
        checkAddLower(s, b);
      }
      Ball *sball = toBall(s);
      BOOST_FOREACH(Ball *b, uppers) {
        b->checkAddLower(sball);
      }
    }

    // Add the new vertex to each higher-dimensional mesh.
    BOOST_FOREACH(const GenericMesh *m, upmeshes) {
      const VertexSet& vset = upvertices.findOrDie(m);
      Vertex *closest = *std::min_element(vset.begin(), vset.end(),
          typename Vertex::IsCloser(toinsert));
      closest->addUninserted(m, toinsert);
    }

    // Reassign uninserted points in the current mesh to the new vertex.
    {
      hudson::FastHashSet<Vertex*, hudson::HashPointer<Vertex> > verts;
      BOOST_FOREACH(const Simplex& s, cavity) {
        BOOST_FOREACH(Vertex *v, s) {
          if (verts.insert(v).second) {
            v->reassignTo(this, toinsert);
          }
        }
      }
    }

    // Unlink the cavity.
    BOOST_FOREACH(const Simplex& s, cavity) {
      unlink(s);
    }

    // Set new handles for each vertex, if we're in top-dimension.  This may
    // set the handle for a vertex several times.  Doesn't cost anything.
    if (topo == ambient) {
      BOOST_FOREACH(const Simplex& s, star) {
        BOOST_FOREACH(Vertex *v, s) {
          v->replaceHandle(toBall(s));
        }
      }
    }

    // Set the return data.
    data.insertKill(this, cavity.begin(), cavity.end());
    data.insertBorn(this, star.begin(), star.end());
    data.setVertex(toinsert);

    // We inserted something, so say so.
    return true;
  }

  private:
  const Simplex& toSimplex(const Ball *b) const {
#ifndef NDEBUG
    // This cast should never fail, so only test it at run-time in debug mode.
    assert(b);
    assert(dynamic_cast<const PBall<Mesh>*>(b));
#endif
    return static_cast<const PBall<Mesh>*>(b)->toSimplex();
  }

  public:
  bool split(Vertex *v, SplitData& data) {
    BOOST_STATIC_ASSERT(topo == ambient);
    assert(v->isCrowded(this));

    // Collect the orbit around the vertex.
    vector<Simplex> orbit;
    Simplicial::collectStar(getSimplicialComplex(), v,
        toSimplex(v->getHandle()), std::back_inserter(orbit));

    if (v->getCD() == topo) {
      // this is a Steiner point; there can't be any nearby uninserted points,
      // or we'd have either warped to it or encroached on something.
      // So just find any point and insert it.  To do that, we need to find a
      // simplex encroached by the uninserted point.
      Vertex *uninserted = v->getFarVertex(this);
      BOOST_FOREACH(const Simplex& s, orbit) {
        if (getDelaunay()->inSphereAffine(s, uninserted->toPoint())) {
          return split(s, data, uninserted);
        }
      }
      // *someone* must include that vertex in its ball
      abort();
    } else {
      // This is an "input" point (it lies on lower-dimensional input).
      // Insert the farthest Voronoi node -- i.e. split the largest-radius simplex.
      size_t best = 0;
      double best_r2 = circumcenter(orbit[best]).radius2;
      size_t n = orbit.size();
      for(size_t i = 0; i < n; ++i) {
        double r2 = circumcenter(orbit[i]).radius2;
        if (r2 > best_r2) {
          best = i;
          best_r2 = r2;
        }
      }
      return split(orbit[best], data);
    }
  }


  /***************************************************************************
   * Check if s and b intersect; if so, link s and b.
   ***************************************************************************/
  bool checkAddLower(const Simplex& s, Ball *b) {
    if (approximatelyIntersects(s, b)) {
      dprintf("  intersection: balls %s and %s\n",
          toBall(s)->toString().c_str(), b->toString().c_str());
      getDataRW(s).blindlyAddLower(b);
      b->blindlyAddUpper(toBall(s));
      return true;
    }
    dprintf("  non-intersection: balls %s and %s\n",
        toBall(s)->toString().c_str(), b->toString().c_str());
    return false;
  }

  /***************************************************************************
   * Add an up pointer.  Does no checking to see whether this is appropriate.
   ***************************************************************************/
  public:
  void blindlyAddUpper(const Simplex& s, Ball *up) {
    assert(up->topological() > topo);
    getDataRW(s).blindlyAddUpper(up);
  }

  /***************************************************************************
   * Add a down pointer.  Does no checking to see whether this is appropriate.
   ***************************************************************************/
  public:
  void blindlyAddLower(const Simplex& s, Ball *lower) {
    dprintf("lower d=%u %s\nto d=%u %s\n",
        (unsigned)lower->topological(), lower->toString().c_str(),
        (unsigned)topo, s.toString().c_str());
    assert(lower->topological() < topo);
    getDataRW(s).blindlyAddLower(lower);
  }

  /***************************************************************************
   * Remove an up pointer.  Does no checking to see whether this is appropriate.
   ***************************************************************************/
  public:
  void removeUpper(const Simplex& s, Ball *up) {
    assert(up->topological() > topo);
    getDataRW(s).removeUpper(up);
  }

  /***************************************************************************
   * Remove a down pointer.  Does no checking to see whether this is appropriate.
   ***************************************************************************/
  public:
  void removeLower(const Simplex& s, Ball *lower) {
    assert(lower->topological() < topo);
    getDataRW(s).removeLower(lower);
  }


  /***************************************************************************
   * Remove all pointers to this simplex.
   ***************************************************************************/
  private:
  void unlink(const Simplex& s) {
    getDataRW(s).unlink();
  }

  /***************************************************************************
   * \name Initialization functions.
   *
   * Iterate over the entire mesh. Used in initialization.
   ***************************************************************************/
  private:
  struct BallNotifier {
    BallNotifier(Mesh *mesh, Ball *b)
      : mesh_(mesh), b_(b) { }
    Mesh *mesh_;
    Ball *b_;
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& s) const {
      mesh_->checkAddLower(s, b_);
      return false;
    }
  };

  public:
  void checkAddLowerToAll(Ball *b) {
    assert(b->topological() < topo);
    BallNotifier notifier(this, b);
    dfs(notifier);
  }

  private:
  const SimplexData& getData(const Simplex& s) const { return del_->getData(s); }
  SimplexData& getDataRW(const Simplex& s) { return del_->getDataRW(s); }

  Delaunay *del_;
};

#undef dprintf

#endif

---