PLCRefiner.h

Go to the documentation of this file.

#ifndef PLC_REFINER_HEADER
#define PLC_REFINER_HEADER


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

#include "MeshConstructor.h"
#include "MeshSet.h"
#include "MeshTypes.h"
#include "SplitData.h"
#include "WorkItem.h"
#include "WorkQueue.h"
#include <boost/function_output_iterator.hpp>
#include <delaunay/DelaunayHelpers.h>
#include <io/PLC.h>
#include <iterator>
#include <refine-common/RefineConstants.h>
#include <refine-common/Refiner.h>
#include <topology/SimplicialHelpers.h>
#include <utilities/OrderedHash.h>

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

namespace SVR {

/***************************************************************************
 * The refiner itself.
 *
 * Read <i>Hudson, Miller, Phillips 2006</i> to understand at a high level what
 * is going on here.
 *
 * The refiner has three main tasks:
 * -# initialization.  Create meshes for each dimension and link them up; then
 *    fill the work queue.  The creation work is mostly off-loaded to
 *    MeshConstructor, with the exception of encroachment fixing.
 * -# refinement.  While the work queue is non-empty, pop off the next element
 *    and process it.  The work queue is handled by the WorkQueue class.
 *    Much of the processing is handled by the ProtectedBall class.  There is
 *    some level of schizophrenia as to who handles what.
 * -# output.  Iterate over all the vertices ever created and output them;
 *    then, iterate over all the simplices in the top-level mesh and output
 *    them.  This is mostly handled in this very class.
 ***************************************************************************/
template <size_t ambient>
class PLCRefiner : public MeshRefiner {

  typedef Geometry::CenterRadius<ambient> CenterRadius;
  typedef Geometry::Point<ambient> Point;
  typedef ::MeshSet<ambient> MeshSet;
  typedef ::SplitData<ambient> SplitData;
  typedef ::WorkItem<ambient> WorkItem;
  typedef ::WorkQueue<ambient> WorkQueue;
  typedef typename MeshTypes<ambient>::Ball Ball;
  typedef typename MeshTypes<ambient>::Vertex Vertex;
  typedef typename MeshTypes<ambient>::GenericMesh GenericMesh;
  typedef typename Ball::BallSet BallSet;
  typedef typename Mesh<ambient, ambient>::Complex Complex;
  typedef typename Complex::Simplex Simplex;
  typedef typename Complex::OSimplex OSimplex;
  typedef typename Vertex::VertexUniqueList VertexUniqueList;

  const RefineConstants consts_;

  MeshSet meshes_; //<! A top-level mesh, and its lower-d submeshes.
  WorkQueue queue_; //<! Handle all the queues of protected balls.
  size_t nextID_; //<! Next ID to allocate for a vertex actually inserted into the mesh.
  unsigned nsimplices_; //<! Number of simplices in the top-dimensional mesh.  I keep this only because it saves me one pass over the mesh during output.
  std::vector<Vertex*> bbox_corners_;

  hudson::SmallList<Simplex> simplices_out_;

  unsigned niterations_; //<! Number of iterations of the main loop so far.
  unsigned message_niterations_; //<! Print a progress message at this iteration.

  private:
  Mesh<ambient,ambient> *getTopMesh() {
    return meshes_.getTopMesh();
  }
  const Mesh<ambient,ambient> *getTopMesh() const {
    return meshes_.getTopMesh();
  }

  public:
  size_t getID(const Vertex *v) const {
    return v->getID();
  }
  void setID(Vertex *v) {
    v->setID(nextID_++);
  }


  /***************************************************************************
   * Initialize a refiner.
   *
   * This is a static function instead of a constructor to allow munging the
   * constants a bit.  See the constructor for details.
   ***************************************************************************/
  public:
  static PLCRefiner *initialize(const SVR_IO::PLC *plc,
      RefineConstants consts) {
    consts.regularize(ambient, false);
    return new PLCRefiner(plc, consts);
  }

  /***************************************************************************
   * Initialize the set of meshes and the work queue, given the PLC.
   *
   * We use MeshConstructor to do most of the work of actual initialization.
   * We then use MC again to iterate over all the protected balls in the mesh,
   * and throw them into the queue.
   *
   * Finally, we initialize the variables we use during output.
   ***************************************************************************/
  private:
  PLCRefiner(const SVR_IO::PLC *plc, const RefineConstants& consts)
    : consts_(consts)
  {
    /* Build the meshes in each dimension, linking up their protective balls. */
    dprintf("Init from PLC\n");
    MeshConstructor<ambient>::construct(plc, consts_.rho, meshes_);

    /* Make an output iterator that calls addToQueues; then invoke it to run on
     * every ball in any dimension. */
    dprintf("Filling queues\n");
    meshes_.outputAllBalls(
        boost::make_function_output_iterator (
          std::bind1st(
            std::mem_fun( &PLCRefiner::checkBall )
            , this)
          )
        );
    Simplicial::collectVertices(
        getTopMesh()->getSimplicialComplex(),
        boost::make_function_output_iterator (
          std::bind1st(
            std::mem_fun( &PLCRefiner::checkVertex )
            , this)
          )
        );
    vector<Vertex*> verts;
    Simplicial::collectVertices(getTopMesh()->getSimplicialComplex(), verts);
    bbox_corners_.reserve(verts.size() - 1);
    BOOST_FOREACH(Vertex *v, verts) {
      if (v->getCD() == ambient) {
        bbox_corners_.push_back(v);
      }
    }

    /* Count the number of top-dimensional simplices. */
    dprintf("Counting top-dimensional balls\n");
    nsimplices_ = meshes_.template outputSomeBalls<ambient>(hudson::count_iterator<Ball*>(0));
    dprintf("  %u\n", nsimplices_);

    /* Find the next free ID to assign. */
    dprintf("Counting vertices\n");
    nextID_ = meshes_.getVertices().size() + plc->getFirstID(0);

    /* No iterations so far. */
    niterations_ = 0;
    message_niterations_ = 1;
  }

  private:
  struct EncroachmentChecker {
    WorkQueue *queue_;
    EncroachmentChecker(WorkQueue *q): queue_(q) {}
    void operator() (Ball *b) {
      BallSet uppers;
      b->gatherUppers(uppers);
      BOOST_FOREACH(Ball *up, uppers) {
        VertexUniqueList verts;
        up->gatherVertices(verts);
        BOOST_FOREACH(Vertex *v, verts) {
          if (b->encroachedBy(v)) {
            queue_->push( new typename WorkItem::Encroached(b) );
          }
        }
      }
    }
  };

  void fixInitialEncroachment() {
    WorkItem *item;
    EncroachmentChecker echeck(&queue_);
    meshes_.outputAllBalls(boost::make_function_output_iterator(echeck));
    dprintf("Start fixing initial encroachments.\n");
    while( (item = queue_.pop()) ) {
      if (item->getReason() > WorkItem::ENCROACHED) {
        /* Done dealing with initial encroachment. */
        queue_.push(item);
        break;
      } else {
        SplitData sdata(consts_);
        split(item, sdata);
        BOUNDED_FOREACH(Ball *b, sdata.begin_born(), sdata.end_born()) {
          echeck(b);
        }
      }
    }
    dprintf("Done fixing initial encroachments.\n");
  }

  friend class CheckFinal;
  struct CheckFinal {
    PLCRefiner *svr_;
    bool isGood;
    CheckFinal(PLCRefiner *svr): svr_(svr), isGood(true) { }
    void operator() (Ball *b) {
      WorkItem *w = svr_->classify(b);
      if (w) {
        fprintf(stderr, "Checking: bad simplex causes %s\n", w->toString().c_str());
        isGood = false;
      }
    }
  };

  void finalCheck() {
    CheckFinal final_mesh(this); // name chosen to make the assertion pretty
    meshes_.outputAllBalls(boost::make_function_output_iterator(final_mesh));
    assert(final_mesh.isGood);
  }

  /***************************************************************************
   * Run the SVR algorithm.
   *
   * First, we eliminate initial encroachment.  The published algorithm
   * requires that there be no initial encroachment, but we can easily fix them.
   * We lose runtime guarantees when we do this, but in practice it should be
   * fast on almost all inputs.
   *
   * Next, we iterate over the work queue until it is done.  This is the core of
   * SVR.
   *
   * At the end, we optionally remove all simplices that are outside the shape;
   * it had better be watertight if we try.
   *
   * Unless NDEBUG is on, we also iterate over all simplices and make sure we
   * didn't forget to do any work.
   ***************************************************************************/
  public:
  virtual void refine() {
    dprintf("Refining\n");

    fixInitialEncroachment();

    WorkItem *w;
    while( (w = queue_.pop()) ) {
      SplitData sdata(consts_);
      split(w, sdata);
    }

    if (consts_.cleanup) {
      cleanup();
    }

#ifndef NDEBUG
    finalCheck();
#endif

    printStats();
  }

  /***************************************************************************
   * Print some basic statistics.
   ***************************************************************************/
  public:
  void printStats() {
    printf("Output size %u simplices\n", nsimplices_);
  }


  /***************************************************************************
   * Basic SVR operation: process a work item.
   *
   * The work of this operation is divided.  The work item farms out the work
   * of finding a point to insert to the ball or vertex we must split, and
   * performs the Delaunay insertion, and updates the point location structures.
   *
   * This class does the work of finding new work to do based on the new
   * structures: we check for encroachment, skinnyness, etc.  Not published in
   * the 2006 IMR paper: we also deal with new vertices on the boundary and push
   * them up to higher dimensions immediately as the need arises.
   ***************************************************************************/
  public:
  void split(WorkItem *w, SplitData& sdata) {
    niterations_++;
    if (niterations_ >= message_niterations_) {
      printf("Iteration %u, we have %u simplices so far\n", niterations_,
          nsimplices_);
      message_niterations_ *= 2;;
    }
    dprintf("Iteration %u: %s (%s)\n", niterations_,
        w->isActive() ? "*splitting*" : "moot",
        w->toString().c_str());

    /* First, check if w is still active. */
    if (!w->isActive()) {
      delete w;
      return;
    }

    size_t topo = w->topological();
    bool didsplit = w->split(sdata);

    /* If we encroached, push the encroached items on the queue. */
    BOUNDED_FOREACH(Ball *encroached, sdata.begin_encroached(), sdata.end_encroached()) {
      queue_.push( new typename WorkItem::Encroached(encroached) );
    }

    /* If we actually did split, do all sorts of new work.  TODO factor this out. */
    if (didsplit) {

      /* If the vertex is new, set its ID and add to the list of vertices.  */
      Vertex *inserted = sdata.getVertex();
      bool isnew = inserted->getCD() == topo;
      if (isnew) {
        meshes_.add(inserted);
        setID(inserted);
        dprintf("  %s is a new vertex in dimension %u at %s\n",
            inserted->toString().c_str(), unsigned(inserted->getCD()),
            inserted->toPoint().toString().c_str());
      }

      /* If we inserted the vertex into the top-d mesh, check if it's crowded. */
      if (topo == ambient) {
        checkVertex(inserted);
      }

      /* Look at the top-d uppers, possibly push their vertices if they are
       * now crowded and weren't before (TODO: I just push if they're crowded
       * now, regardless whether they were before) */
      if (topo < ambient && isnew) {
        BallSet uppers;
        VertexUniqueList vset;
        BOUNDED_FOREACH(Ball *killed, sdata.begin_killed(), sdata.end_killed()) {
          killed->gatherUppers(uppers);
        }
        BOOST_FOREACH(Ball *upper, uppers) {
          if (upper->topological() != ambient) { continue; }
          upper->gatherVertices(vset);
        }
        BOOST_FOREACH(Vertex *upv, vset) {
          checkVertex(upv);
        }
      }

      /* Look at one of the simplices we split (any of them), call it s.
       * Look at the higher-d meshes that know about s.  For each upper
       * mesh m, if all of the vertices of s appear in m, add a new workqueue
       * item for an arbitrary s' in m telling it to split at v. */
      if (topo < ambient && isnew) {
        Ball *s = *sdata.begin_killed();
        // compute the uppers, split them up by mesh
        BallSet alluppers;
        hudson::OrderedHashSet<const GenericMesh*, hudson::HashPointer<GenericMesh> > meshes;
        hudson::FastHashMap<const GenericMesh*, BallSet, hudson::HashPointer<GenericMesh> > uppers;
        s->gatherUppers(alluppers);
        BOOST_FOREACH(Ball *upper, alluppers) {
          const GenericMesh *m = upper->getMesh();
          if (m->topological() == ambient) { continue; }
          uppers[m].insert(upper);
          meshes.insert(m);
        }
        VertexUniqueList verts;
        s->gatherVertices(verts);
        BOOST_FOREACH(const GenericMesh *m, meshes) {
          dprintf("  checking for immediate inclusion in %u-mesh[%p]...",
              m->topological(), m);
          VertexUniqueList mverts;
          BOOST_FOREACH(Ball *upper, uppers[m]) {
            upper->gatherVertices(mverts);
          }
          // check if s has all verts resolved
          bool allresolved = true;
          BOOST_FOREACH(Vertex *v, verts) {
            if (!mverts.has(v)) {
              dprintf("no\n");
              allresolved = false;
              break;
            }
          }
          if (allresolved) {
            // Find a simplex to blame.  The 'uppers' set is unrobust, but at
            // least one will be robustly good, so find it.
            Ball *topush = NULL;
            BOOST_FOREACH(Ball *upper, uppers[m]) {
              if (upper->inSphereExact(inserted)) {
                topush = upper;
                break;
              }
            }
            dprintf("split at %s\n", topush->toString().c_str());
            assert(topush);
            queue_.push( new typename WorkItem::NewBoundaryVertex(topush, inserted) );
          }
        }
      }

      /* Look at the new simplices, classify them, and possibly push them on the
         work queue skinny or sliver. */
      BOUNDED_FOREACH(Ball *newborn, sdata.begin_born(), sdata.end_born()) {
        checkBall(newborn);
      }

      /* Update the count of simplices. */
      if (topo == ambient) {
        nsimplices_ += sdata.nborn();
        nsimplices_ -= sdata.nkilled();
        dprintf("There are now %u simplices\n", unsigned(nsimplices_));
      }

      // TODO: put this somewhere useful
#if 0
      /* Animation hack: print out every triangle that's all the same feature. */
      if (topo == ambient && inserted->getCD() < ambient) {
        animation_print();
      }
#endif
    }

    /* Finally, if the work item is still active, push it back on.
       If not, delete it. */
    if (w->isActive()) {
      queue_.push(w);
    } else {
      delete w;
    }
  }


  /***************************************************************************
    * Printing, in my own special little format easy to perlify into other formats.
   ***************************************************************************/
  typedef std::pair<boost::array<unsigned, ambient>, unsigned> face_type;
  struct animator {
    hashers::hash_map<Vertex*, unsigned, hudson::hash_by_cast<Vertex*> > ids;
    hashers::hash_map<const GenericMesh*, unsigned, hudson::hash_by_cast<const GenericMesh*> > colours;

    std::vector<Vertex*> verts;
    std::vector<face_type> faces;
    const Complex& c;

    animator(const Complex& c) : c(c) { }

    static bool canEnter(const Simplex&) { return true; }

    bool choose(const Simplex& s) {
      OSimplex os = c.orientSimplex(s);
      for(unsigned i = 0; i <= ambient; ++i) { // for each of d+1 vertices of s
        processFace(c.getFace(i, os));
      }
      return false;
    }

    void processFace(const OSimplex& fi) {
      // Collect the intersection set of meshes that contain these points.
      // This uses the assumption that the meshes are stored in sorted order.
      std::vector<const GenericMesh*> meshes;
      std::copy(fi[0]->begin_uppers(ambient-1), fi[0]->end_uppers(ambient-1),
          std::back_inserter(meshes));
      for(unsigned j = 1; j < ambient; ++j) { // for each of d vertices of fi
        std::vector<const GenericMesh*> tmp;
        std::set_intersection(meshes.begin(), meshes.end(),
            fi[j]->begin_uppers(ambient - 1), fi[j]->end_uppers(ambient - 1),
            std::back_inserter(tmp));
        meshes.swap(tmp);
      }
      assert(meshes.size() == 0 || meshes.size() == 1);
      if (!meshes.empty()) {
        face_type face;
        for(unsigned j = 0; j < ambient; ++j) {
          face.first[j] = getVertexID(fi[j]);
        }
        face.second = getMeshColour(meshes[0]);
        faces.push_back(face);
      }
    }
    unsigned getVertexID(Vertex *v) {
      if(ids.find(v) == ids.end()) {
        unsigned id = ids.size();
        ids[v] = id;
        verts.push_back(v);
      }
      return ids[v];
    }

    unsigned getMeshColour(const GenericMesh *m) {
      if (colours.find(m) == colours.end()) {
        unsigned col = colours.size();
        colours[m] = col;
      }
      return colours[m];
    }
  };


  void animation_print() {
    /* Collect the faces. */
    const Complex& c = getTopMesh()->getSimplicialComplex();
    animator anim(c);
    c.dfsBySimplex(anim, c.getHandle());
    assert(ambient == 3);

    unsigned nvertices = anim.verts.size();
    unsigned nfaces = anim.faces.size();
    FILE *out;

    {
      char buffer[4096];
      snprintf(buffer, sizeof(buffer), "faces-%u.off", nvertices);
      out = fopen(buffer, "w");
    }
    assert(out);

    /* .off format */
    fprintf(out, "%u %u 0\n", nvertices, nfaces);

    BOOST_FOREACH(const Vertex *v, anim.verts) {
      Point p = v->toPoint();
      fprintf(out, "%g %g %g\n", p[0], p[1], p[2]);
    }
    BOOST_FOREACH(face_type face, anim.faces) {
      // three vertices, plus the colour.
      fprintf(out, "3 %u %u %u %u\n", face.first[0], face.first[1], face.first[2], face.second);
    }
    fclose(out);
  }

  /***************************************************************************
    * Cleanup after meshing: remove the exterior.
   ***************************************************************************/
  void cleanup() {
    bool watertight = MeshCleaner<ambient,ambient>::clean(
        getTopMesh()->getDelaunayDangerously(),
        bbox_corners_.begin(), bbox_corners_.end(),
        &simplices_out_);

    if (!watertight) {
      nsimplices_ = 0;
    } else {
      nsimplices_ = std::distance(simplices_out_.begin(), simplices_out_.end());
    }

    // In a perfect world, delete the excess vertices now.  Instead, I just
    // leak them.
  }

  /***************************************************************************
   ** Printing the top-dimensional mesh.
   ** TODO: print lower-dimensional ones (for debugging) also.
   ***************************************************************************/
  private:
  friend struct DoPrintSimplex;
  struct DoPrintSimplex {
    size_t n_;
    FILE *out_;
    const PLCRefiner *svr_;
    typedef typename Mesh<ambient, ambient>::Simplex Simplex;

    DoPrintSimplex(const PLCRefiner *svr, FILE *out): n_(0), out_(out), svr_(svr) { }
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& s) {
      boost::array<size_t, ambient + 1> arr;
      for (size_t i = 0; i <= ambient; ++i) {
        arr[i] = svr_->getID(s[i]);
      }
      Delaunay::printEle(n_++, arr, out_);
      return false;
    }
  };

  template <class v_iterator> void printVerts(FILE *node, v_iterator begin, v_iterator end) const {
    while (begin != end) {
      Vertex *v = *begin;
      ++begin;
      Delaunay::printPoint(getID(v), v->toPoint(), node);
    }
  }

  public:
  virtual void printToNodeEle(FILE *node, FILE *ele) const {
    /* print the .node file */
    fprintf(node, "%u %u 0 0\n", (unsigned)meshes_.getVertices().size(), unsigned(ambient));
    printVerts(node, meshes_.getVertices().begin(), meshes_.getVertices().end());

    /* print the .ele file */
    fprintf(ele, "%u %u 0\n", unsigned(nsimplices_), unsigned(ambient + 1));
    DoPrintSimplex dps(this, ele);
    if (consts_.cleanup) {
      BOOST_FOREACH(const Simplex& s, simplices_out_) {
        dps.choose(s);
      }
    } else {
      const Complex& c = getTopMesh()->getSimplicialComplex();
      c.dfsBySimplex(dps, c.getHandle());
    }
  }


  /***************************************************************************
    * Classify a ball's reason to be added to the work queue.
    * The reason may be NONE, in which case we won't add it.
   ***************************************************************************/
  private:
  WorkItem *classify(Ball *b) const {
    /* Assume unencroached; we check elsewhere. */

    double re2 = b->radiusEdge2();
    if (re2 > consts_.rhoprime * consts_.rhoprime) {
      return new typename WorkItem::VerySkinny(b);
    }

    if (consts_.sigma != 0.0 && b->topological() >= 3) {
      double sigma = b->computeSigma();
      if (sigma < consts_.sigma) {
        return new typename WorkItem::Sliver(b);
      }
    }

    if (re2 > consts_.rho * consts_.rho) {
      return new typename WorkItem::MediumSkinny(b);
    }

    /* Check if the simplex is supposed to be represented in the global mesh;
     * if so, and it's not, complain about that. */
    if (!b->isResolved()) {
      return new typename WorkItem::Unresolved(b);
    }

    return NULL;
  }

  /***************************************************************************
   * Check if a ball needs to be added to the work queue; if so, add it.
   ***************************************************************************/
  void checkBall (Ball *b) {
    WorkItem *w = classify(b);
    if (w) {
      queue_.push(w);
    }
  }

  void checkVertex (Vertex *v) {
    if (v->isCrowded(getTopMesh())) {
      if (v->getCD() == ambient) {
        queue_.push( new typename WorkItem::CrowdedSteiner(v, getTopMesh()) );
      } else {
        queue_.push( new typename WorkItem::CrowdedInput(v, getTopMesh()) );
      }
    }
  }

  bool isActive(const WorkItem *item) const {
    if (!item->isBall()) {
      return item->getVertex()->isCrowded(getTopMesh());
    } else {
      Ball *b = item.getBall();
      if (!b->isLive()) { return false; }
      if (item.getReason() == WorkQueue::UNRESOLVED) {
        return !isResolved(b);
      }
      return true;
    }
  }

}; // class SVR::FeatureRefiner<d>
} // namespace SVR

#undef dprintf


#endif

---