IncrementalDelaunay.h

Go to the documentation of this file.

00001 #ifndef IncrementalDelaunay_HEADER
00002 #define IncrementalDelaunay_HEADER
00003 
00004 /***************************************************************************
00005  * Benoit Hudson and Todd Phillips (C) 2007.
00006  ***************************************************************************/
00007 
00008 #include <boost/array.hpp>
00009 #include <geometry/BoundingBox.h>
00010 #include <geometry/Circles.h>
00011 #include <geometry/Point.h>
00012 #include <geometry/PointConvert.h>
00013 #include <geometry/ProjectionPlane.h>
00014 #include <geometry/predicates.h>
00015 #include <topology/SimplicialComplex.h>
00016 #include <topology/SimplicialHelpers.h>
00017 #include <geometry/VertexToPointIterator.h>
00018 #include "SimplexCenter.h"
00019 
00021 #define DEBUG_DELAUNAY 0
00022 #if 0
00023 #define dprintf(...) printf(__VA_ARGS__);
00024 #else
00025 #define dprintf(...)
00026 #endif
00027 
00028 /***************************************************************************
00029  * A class supporting incremental insertions into a Delaunay triangulation.
00030  *
00031  * Because of the usage in SVR, it makes sense to hold the center and radius
00032  * within each simplex.  We also can attach any additional user data may be
00033  * useful per simplex.
00034  *
00035  * Internally, we use the SimplicialComplex from <topology/SimplicialComplex.h>
00036  * Simplices are oriented in CCW or right-hand-rule order (see isOriented).
00037  *
00038  * Class Vertex must match:
00039  * - Vertex(Point) constructor
00040  * - Point toPoint() const
00041  *
00042  * The 1-d code is specialized in Delaunay1.h , included at the end of this file.
00043  ***************************************************************************/
00044 template <unsigned ambient, unsigned topological,
00045   class Vertex_, class SimplexData = void,
00046   class VertexPrinter = typename Vertex_::Printer>
00047 class IncrementalDelaunay {
00048 public:
00049 typedef Geometry::Point<ambient> Point; 
00050 typedef Geometry::Point<topological> PPoint; 
00051 typedef Geometry::ProjectionPlane<ambient, topological> ProjectionPlane;
00052 typedef Geometry::BoundingBox<topological> BoundingBox;
00053 typedef Geometry::CenterRadius<ambient> CenterRadius;
00054 
00055 typedef Vertex_ Vertex;
00056 private:
00057 typedef hudson::PayloadBucket<CenterRadius, SimplexData> SimplexStorage;
00058 
00059 public:
00060 typedef SimplicialComplex<topological, Vertex, SimplexStorage, VertexPrinter> Complex;
00061 typedef typename Complex::Simplex Simplex;
00062 typedef typename Complex::OSimplex OSimplex;
00063 typedef typename Complex::Cavity Cavity;
00064 typedef typename Complex::Star Star;
00065 
00066 private:
00067 Complex complex_;
00068 ProjectionPlane projectionPlane_;
00069 vector<Vertex*> myVertices_;
00070 
00071 public:
00074 ~IncrementalDelaunay() {
00075   // Delete all the vertices we own.
00076   BOOST_FOREACH(Vertex *v, myVertices_) {
00077     delete v;
00078   }
00079   // The complex goes away automatically.
00080 }
00081 
00082 /***************************************************************************
00083  * \name Constructors.
00084  * - fromSimplex(): takes d+1 points and produces a simplex.
00085  * - makeBoundingSimplex(): takes in any number of vertices and produces
00086  *      a simplex that contains the vertices.
00087  * - createSquare(): takes in any number of vertices and produces a 2^d-point
00088  *      box that contains the vertices.
00089  * - create4box(): takes in any number of vertices and produces a box
00090  *      useful for meshing: it centers the min/max box in a box 4L on
00091  *      a side, split by a point every L (thus, four points per side of the
00092  *      square, 16 points per face of the cube, ...)
00093  * TODO: properly document them all.
00094  ***************************************************************************/
00096 public:
00100 IncrementalDelaunay(const vector<Vertex*>& verts,
00101                     const vector<vector<size_t> >& indices,
00102                     const ProjectionPlane& pplane)
00103   : complex_(verts, indices), projectionPlane_(pplane)
00104 {
00105   if (!isOriented(getHandle())) {
00106     complex_.flipOrientation();
00107   }
00108   assert(isOriented());
00109 }
00110 
00111 public:
00112 static IncrementalDelaunay *fromSimplex(const vector<Vertex*>& verts, const ProjectionPlane& plane) {
00113    assert(verts.size() == topological + 1);
00114    vector<size_t> thesimplex;
00115    thesimplex.resize(topological + 1);
00116    for(unsigned i = 0; i < topological + 1; i++) {
00117      thesimplex[i] = i;
00118    }
00119    vector<vector<size_t> > elements;
00120    elements.push_back(thesimplex);
00121    return new IncrementalDelaunay(verts, elements, plane);
00122  }
00123 
00124 static IncrementalDelaunay *fromSimplex(const vector<Vertex*>& verts) {
00125    ProjectionPlane plane (verts.begin(), verts.end(), true);
00126    return fromSimplex(verts, plane);
00127  }
00129 
00130 
00131 /***************************************************************************
00132  * \name Data access
00133  * We have to be careful with what we return, in order to handle void data.
00134  * We forward the request on to the complex, and then we pull out the
00135  * user's part of it.
00136  ***************************************************************************/
00137 typedef typename hudson::Bucket<SimplexData>::reference data_ref;
00138 typedef typename hudson::Bucket<SimplexData>::const_reference data_const_ref;
00139 data_ref getDataRW(const Simplex& s) {
00140   return complex_.getDataRW(s).payload_ref();
00141 }
00142 
00144 data_const_ref getData(const Simplex& s) const {
00145   return complex_.getData(s).payload_cref();
00146 }
00147 
00148 /***************************************************************************
00149  * \name Dimension
00150  * There are two notions of dimension.  We require callers to specify
00151  * which one they want.
00152  * - The ambient dimension is the dimension in which the points live.
00153  * - The topological dimension is that the triangulation lives on a sub-space
00154  * of the ambient dimension.
00155  *
00156  * Everything is of course much easier when the two notions are equal.
00157  ***************************************************************************/
00158 public:
00159 static unsigned ambientDimension() {
00160   return ambient;
00161 }
00162 static unsigned topologicalDimension() {
00163   return topological;
00164 }
00165 
00166 const ProjectionPlane& getPlane() const {
00167   return projectionPlane_;
00168 }
00169 
00170 /***************************************************************************
00171  * \name Geometric computations
00172  ***************************************************************************/
00174 private:
00176 CenterRadius
00177 do_circumcenter(const Simplex& simplex) const {
00178   assert(isOriented(simplex));
00179   return SimplexCenter<Simplex, ambient, topological>()(simplex);
00180 }
00181 
00182 public:
00189 const CenterRadius&
00190 circumcenter(const Simplex& s) const {
00191   // violate constness for the sake of the cache
00192   CenterRadius& cr = const_cast<Complex&>(complex_).getDataRW(s).data;
00193   if (topological != 0 && cr.radius2 == 0.0) {
00194     cr = do_circumcenter(s);
00195   }
00196   return cr;
00197 }
00198 
00199 /***************************************************************************
00200  * Compute Ungor's offcenter for the simplex: the point that makes the
00201  * smallest face of s achieve quality rho.
00202  ***************************************************************************/
00203 public:
00204 Point offCenter (const Simplex& s, double rho2) const {
00205   /* TODO!!! */
00206   if (topological != 2) {
00207     return circumcenter(s).center;
00208   }
00209 
00210   /* We want to walk x units along centerVect starting at midpoint,
00211    * until the point forms a well-shaped triangle with e.
00212    *
00213    * We first should walk y units to the circumcenter of the new
00214    * triangle, with radius r = rho * |e|.  y^2 = r^2 - (e/2)^2 = (rho^2 - 1/4) e^2.
00215    * Then we walk another r units to find the apex rather than the
00216    * center.   x = y + r = (sqrt(rho^2-1/4) + rho) e.
00217    *
00218    * If the offcenter goes past the circumcenter, we use the circumcenter instead.
00219    *
00220    * We use rho - epsilon as the target quality, so we don't have rounding error problems.
00221    */
00222   // hack: avoid rounding bugs
00223   rho2 -= 1e-8;
00224 
00225   pair<unsigned, unsigned> e = shortestEdge(s);
00226   const CenterRadius& cr = circumcenter(s);
00227   Point a = s[e.first]->toPoint();
00228   Point b = s[e.second]->toPoint();
00229   Point midpoint = (a+b) / 2.0;
00230   double elen2 = (a-b).norm2();
00231   Point centerVect = cr.center - midpoint;
00232   double norm2 = centerVect.norm2();
00233 
00234   double x2 = 2.0 * elen2 * (rho2 - 0.125 + sqrt(rho2 * (rho2 - 0.25)));
00235 
00236   if (norm2 < x2) {
00237     dprintf("off-center rejected: %g < %g\n", sqrt(norm2), sqrt(x2));
00238     return cr.center;
00239   } else {
00240     Point offcenter = centerVect * sqrt(x2 / norm2);
00241     dprintf("off-center accepted: %s rather than %s\n",
00242         Point(offcenter + midpoint).toString().c_str(),
00243         cr.center.toString().c_str());
00244     return offcenter + midpoint;
00245   }
00246 }
00247 
00248 
00249 /***************************************************************************
00250  * Find the shortest edge (by index).
00251  ***************************************************************************/
00252 public:
00253 static pair<unsigned, unsigned> shortestEdge (const Simplex& s) {
00254   pair<unsigned, unsigned> best = make_pair(0,1);
00255   double bestlen2 = s[0]->toPoint().dist2(s[1]->toPoint());
00256 
00257   for(unsigned i = 0; i < topological; ++i) {
00258     Point p = s[i]->toPoint();
00259     for(unsigned j = i+1; j <= topological; ++j) {
00260       Point q = s[j]->toPoint();
00261       double len2 = p.dist2(q);
00262       if (len2 < bestlen2) {
00263         bestlen2 = len2;
00264         best = make_pair(i, j);
00265       }
00266     }
00267   }
00268   return best;
00269 }
00270 
00272 static double edgeLen2(const Simplex& s, const pair<unsigned, unsigned>& e) {
00273   return s[e.first]->toPoint().dist2(s[e.second]->toPoint());
00274 }
00275 static double edgeLen2(const Simplex& s) {
00276   // check all edges
00277   double e2 = (s[0]->toPoint() - s[1]->toPoint()).norm2();
00278   for(unsigned i = 0; i < topological; ++i) {
00279     Point p = s[i]->toPoint();
00280     for(unsigned j = i+1; j <= topological; ++j) {
00281       Point q = s[j]->toPoint();
00282       e2 = std::min(e2, (p-q).norm2());
00283     }
00284   }
00285   return e2;
00286 }
00287 
00288 /***************************************************************************
00289  * Compute various quality ratios.
00290  ***************************************************************************/
00291 public:
00292 // Ratio of circumradius to shortest edge.  Small is good.
00293 double radiusEdge2 (const Simplex& s) const {
00294   double r2 = circumcenter(s).radius2;
00295   double e2 = edgeLen2(s);
00296   assert(e2 > 0.0);
00297   return r2/e2;
00298 }
00299 
00300 double volume(const Simplex& simplex) const {
00301   if (topological == ambient) {
00302     switch(ambient) {
00303       case 1: return fabs(simplex[1]->toPoint()[0] - simplex[0]->toPoint()[0]);
00304       case 2: return fabs(orient2d(simplex[0]->x(), simplex[1]->x(), simplex[2]->x())) / 2.0;
00305       case 3: return fabs(orient3d(simplex[0]->x(), simplex[1]->x(), simplex[2]->x(), simplex[3]->x())) / 6.0;
00306       default: assert(0);
00307     }
00308   } else {
00309     boost::array<PPoint, topological+1> s;
00310     projectSimplex(s, simplex);
00311     switch(topological) {
00312       case 1: return fabs(s[1][0] - s[0][0]);
00313       case 2: return fabs(orient2d(s[0].data(), s[1].data(), s[2].data())) / 2.0;
00314       case 3: return fabs(orient3d(s[0].data(), s[1].data(), s[2].data(), s[3].data())) / 6.0;
00315       default: assert(0);
00316     }
00317   }
00318 }
00319 
00320 // Sigma, as defined by Li-Teng.  Small is bad.
00321 double computeSigma (const Simplex& s) const {
00322   assert(ambient == 3);
00323   double V = volume(s);
00324   double e2 = edgeLen2(s);
00325   assert(e2 > 0.0);
00326   double e = sqrt(e2);
00327   double e3 = e*e*e;
00328   return V / e3;
00329 }
00330 
00331 // Ratio of outradius to inradius.  Small is good.
00332 double radiusRadius (const Simplex& simplex) const {
00333   double r2 = circumcenter(simplex).radius2;
00334 
00335   switch(topological) {
00336     case 1:
00337       return 1.0;
00338     case 2:
00339       {
00340         /* From jrs' robnotes: the inscribed circle has radius 2 A / sum(l_i) */
00341 
00342         /* Area is sqrt(s*(s-a)*(s-b)*(s-c))
00343            where s=(a+b+c)/2 is the semi-perimeter. */
00344         double a2 = (simplex[0]->toPoint() - simplex[1]->toPoint()).norm2();
00345         double b2 = (simplex[0]->toPoint() - simplex[2]->toPoint()).norm2();
00346         double c2 = (simplex[1]->toPoint() - simplex[2]->toPoint()).norm2();
00347         double a = sqrt(a2);
00348         double b = sqrt(b2);
00349         double c = sqrt(c2);
00350 
00351         double s = (a + b + c) / 2.0;
00352         double A = sqrt(s*(s-a)*(s-b)*(s-c));
00353 
00354         /* from jrs: 2A / sum l_i */
00355         double rin = 2.0 * A / (a+b+c);
00356         return sqrt(r2) / rin;
00357       }
00358     case 3:
00359       {
00360         /* From jrs' robnotes: the inscribed circle has radius 3 V / sum(A_i) */
00361         double V6 = fabs(orient3d(simplex[0]->x(), simplex[1]->x(), simplex[2]->x(), simplex[3]->x()));
00362         double numerator = V6 / 2.0;
00363         double denominator = 0.0;
00364         OSimplex os = getComplex().orientSimplex(simplex);
00365         for (unsigned i = 0; i <= topological; ++i) {
00366           denominator += triangleArea(getComplex().getFace(i, os));
00367         }
00368         double rin = numerator / denominator;
00369         return sqrt(r2) / rin;
00370       }
00371     default:
00372       assert(0);
00373       return 0.0;
00374   }
00375 }
00376 
00381 double minDihedral (const Simplex& simplex) const {
00382   assert(topological == 3);
00383   assert(ambient == 3);
00384   const Point& a = simplex[0]->toPoint();
00385   const Point& b = simplex[1]->toPoint();
00386   const Point& c = simplex[2]->toPoint();
00387   const Point& d = simplex[3]->toPoint();
00388 
00389   Point ad = a - d;
00390   Point ab = a - b;
00391   Point ac = a - c;
00392   Point bd = b - d;
00393   Point bc = b - c;
00394   Point cd = c - d;
00395 
00396   Point na = crossProduct(cd,bd);
00397   Point nb = crossProduct(ad,cd);
00398   Point nc = crossProduct(bd,ad);
00399   Point nd = crossProduct(ac,bc);
00400 
00401   double V6 = orient3d(a.data(), b.data(), c.data(), d.data());
00402 
00403   double parts[][2] = {
00404     { sqrt(ab.norm2()) , nc.dot(nd) },
00405     { sqrt(ac.norm2()) , nb.dot(nd) },
00406     { sqrt(ad.norm2()) , nb.dot(nc) },
00407     { sqrt(bc.norm2()) , na.dot(nd) },
00408     { sqrt(bd.norm2()) , na.dot(nc) },
00409     { sqrt(cd.norm2()) , na.dot(nb) }
00410   };
00411   const unsigned n = sizeof(parts)/sizeof(*parts);
00412   boost::array<double, n> theta;
00413   double worst = 1e9;
00414   for (unsigned i = 0; i < n; ++i) {
00415     if (parts[i][1] == 0) {
00416       // theta is either -90 or +90, i.e. pi/2
00417       theta[i] = M_PI / 2.0;
00418     } else if (parts[i][1] < 0) {
00419       // theta is in (-90,90)
00420       theta[i] = fabs(atan(V6 * parts[i][0] / parts[i][1]));
00421     } else /*parts[i][1] > 0*/ {
00422       // theta is in (90,180) or (-90,-180)
00423       theta[i] = fabs(atan(V6 * parts[i][0] / parts[i][1])) + M_PI;
00424     }
00425     if (theta[i] < worst) worst = theta[i];
00426   }
00427   assert(worst > 0.0);
00428   return worst;
00429 }
00430 
00434 double triangleArea (const OSimplex& simplex) const {
00435   assert(topological >= 2);
00436   /* sqrt(s*(s-a)*(s-b)*(s-c))
00437      where s=(a+b+c)/2 is the semi-perimeter. */
00438   double a = sqrt((simplex[0]->toPoint() - simplex[1]->toPoint()).norm2());
00439   double b = sqrt((simplex[0]->toPoint() - simplex[2]->toPoint()).norm2());
00440   double c = sqrt((simplex[1]->toPoint() - simplex[2]->toPoint()).norm2());
00441 
00442   double s = (a + b + c) / 2.0;
00443   double A = sqrt(s*(s-a)*(s-b)*(s-c));
00444   return A;
00445 }
00446 
00447 
00448 /***************************************************************************
00449  * Project from the simplex down to the axes-parallel plane.
00450  ***************************************************************************/
00451 void projectSimplex(boost::array<PPoint, topological+1>& s, const Simplex& simplex) const {
00452   for(unsigned i = 0; i <= topological; i++) {
00453     s[i] = projectionPlane_.project(simplex[i]->toPoint());
00454   }
00455 }
00457 
00458 
00459 
00460 /***************************************************************************
00461  * \name Geometric predicates
00462  ***************************************************************************/
00464 
00465 /***************************************************************************
00466  * Is the point in the open circumball (or diametral ball) of the simplex?
00467  ***************************************************************************/
00468 public:
00469 bool
00470 inSphere(const Simplex& s, const Point& p) const {
00471   assert(isOriented(s));
00472   if(ambient == topological) {
00473     BOOST_STATIC_ASSERT(topological > 1);
00474     /* collect the points and send to Shewchuk's code */
00475     switch(ambient) {
00476       case 2:
00477         // incircle is positive if d lies inside
00478         return incircle(s[0]->x(), s[1]->x(), s[2]->x(), p.data()) > 0.0;
00479 
00480       case 3:
00481         // insphere is *negative* if e lies inside.
00482         return insphere(s[0]->x(), s[1]->x(), s[2]->x(), s[3]->x(), p.data()) < 0.0;
00483 
00484       default:
00485         ;
00486         // Fall through to default processing.
00487     }
00488   }
00489 
00490   // TODO: eliminate this fall-through and have only robust predicates.
00491   CenterRadius cr = circumcenter(s);
00492   double cp2 = (cr.center - p).norm2();
00493   return cp2 < cr.radius2;
00494 }
00495 
00496 
00497 /***************************************************************************
00498  * Assuming p is on the plane of this topology, do an in-sphere test.
00499  *
00500  * This is cheaper than the above because we project down first,
00501  * giving us a lower-dimensional in-circle test to run.
00502  * It's also more robust, since we can use precise predicates instead
00503  * of calculating the circumcenter and hoping the radius means something.
00504  ***************************************************************************/
00505 public:
00506 bool
00507 inSphereAffine(const Simplex& simplex, const Point& p) const {
00508   if (topological == ambient || topological > 3) {
00509     // then we don't gain anything by being fancy here: either
00510     // we don't need to project, or we have no predicates.
00511     return inSphere(simplex, p);
00512   }
00513   boost::array<PPoint, topological+1> s;
00514   projectSimplex(s, simplex);
00515   PPoint pp = projectionPlane_.project(p);
00516   BOOST_STATIC_ASSERT(topological == 2 || topological == 3);
00517   switch(topological) {
00518     case 2:
00519       return incircle(s[0].data(), s[1].data(), s[2].data(), pp.data()) > 0.0;
00520 
00521     case 3:
00522       return insphere(s[0].data(), s[1].data(), s[2].data(), s[3].data(), pp.data()) < 0.0;
00523   }
00524   // unreachable code
00525   assert(0);
00526 }
00527 
00528 /***************************************************************************
00529  * Assuming p is on the affine plane of the complex, it is well defined
00530  * to test whether it lies within or outside of the simplex.
00531  *
00532  * We consider the simplex to be closed: the sides are 'in.'
00533  ***************************************************************************/
00534 bool
00535 inSimplexAffine(const Simplex& simplex, const Point& p) const {
00536   if (ambient == topological) {
00537     return inSimplex(simplex, p);
00538   }
00539 
00540   PPoint pp = projectionPlane_.project(p);
00541   boost::array<PPoint, topological+1> s;
00542   projectSimplex(s, simplex);
00543   switch(topological) {
00544     case 1:
00545       // not oriented
00546       return (pp[0] >= s[0][0] && pp[0] <= s[1][0])
00547         || (pp[0] >= s[1][0] && pp[0] <= s[0][0]);
00548     case 2:
00549       // orientation is 012, so 01, 12, and 20 are oriented
00550       // orient2d returns positive on CCW order
00551       return (orient2d(s[0].data(), s[1].data(), pp.data()) >= 0.0)
00552         && (orient2d(s[1].data(), s[2].data(), pp.data()) >= 0.0)
00553         && (orient2d(s[2].data(), s[0].data(), pp.data()) >= 0.0);
00554     case 3:
00555       // orientation is 0123, so 012, 031, 023, 321
00556       // orient3d returns negative on right-hand-rule order
00557       return (orient3d(s[0].data(), s[1].data(), s[2].data(), pp.data()) <= 0.0)
00558         && (orient3d(s[0].data(), s[3].data(), s[1].data(), pp.data()) <= 0.0)
00559         && (orient3d(s[0].data(), s[2].data(), s[3].data(), pp.data()) <= 0.0)
00560         && (orient3d(s[3].data(), s[2].data(), s[1].data(), pp.data()) <= 0.0);
00561     default:
00562       assert(0);
00563       return false;
00564   }
00565 }
00566 
00574 private:
00575 bool
00576 inSimplex(const Simplex& s, const Point& p) const {
00577   assert(ambient == topological);
00578   switch(topological) {
00579     case 1:
00580       // not oriented
00581       return (p[0] >= s[0]->toPoint()[0] && p[0] <= s[1]->toPoint()[0])
00582         || (p[0] >= s[1]->toPoint()[0] && p[0] <= s[0]->toPoint()[0]);
00583     case 2:
00584       // orientation is 012, so 01, 12, and 20 are oriented
00585       // orient2d returns positive on CCW order
00586       return (orient2d(s[0]->x(), s[1]->x(), p.data()) >= 0.0)
00587         && (orient2d(s[1]->x(), s[2]->x(), p.data()) >= 0.0)
00588         && (orient2d(s[2]->x(), s[0]->x(), p.data()) >= 0.0);
00589     case 3:
00590       // orientation is 0123, so 012, 031, 023, 321
00591       // orient3d returns negative on right-hand-rule order
00592       return (orient3d(s[0]->x(), s[1]->x(), s[2]->x(), p.data()) <= 0.0)
00593         && (orient3d(s[0]->x(), s[3]->x(), s[1]->x(), p.data()) <= 0.0)
00594         && (orient3d(s[0]->x(), s[2]->x(), s[3]->x(), p.data()) <= 0.0)
00595         && (orient3d(s[3]->x(), s[2]->x(), s[1]->x(), p.data()) <= 0.0);
00596     default:
00597       assert(0);
00598       return false;
00599   }
00600 }
00601 
00602 /***************************************************************************
00603   * For debugging: Check whether the simplex is properly-oriented.
00604   *
00605   * Also used in the constructor.
00606  ***************************************************************************/
00607 private:
00608 bool
00609 isOriented_(const Simplex& simplex) const {
00610   if (topological == 1) return true;
00611   boost::array<PPoint, topological+1> s;
00612   projectSimplex(s, simplex);
00613   switch (topological) {
00614     case 1:
00615       return s[0][0] < s[1][0];
00616     case 2:
00617       // ccw is positive
00618       return orient2d(s[0].data(), s[1].data(), s[2].data()) > 0.0;
00619     case 3:
00620       // right-hand rule is negative
00621       return orient3d(s[0].data(), s[1].data(), s[2].data(), s[3].data()) < 0.0;
00622     default:
00623       assert(0); // TODO
00624       return false;
00625   }
00626 }
00627 public:
00628 bool
00629 isOriented(const Simplex& simplex) const {
00630   bool b = isOriented_(simplex);
00631 #if DEBUG_DELAUNAY
00632   if (!b) {
00633     boost::array<PPoint, topological+1> s;
00634     projectSimplex(s, simplex);
00635     printf("Unoriented simplex: %s\n", simplex.toString().c_str());
00636     printf("Projects to: ");
00637     BOOST_FOREACH(const PPoint& p, s) {
00638       printf(" %s ", p.toString().c_str());
00639     }
00640     printf("\n");
00641   }
00642 #endif
00643   return b;
00644 }
00645 
00647 
00648 
00649 /***************************************************************************
00650  * \name Point insertion
00651  *
00652  * The algorithm is from Bowyer and Watson.  To insert a point p:
00653  * -# find a simplex that contains p in its circumball
00654  * -# compute p's cavity using computeCavity()
00655  * -# compute the star using computeStar()
00656  * -# commit the star using insertVertex()
00657  *
00658  * You may combine the last three steps using a different variant of
00659  * insertVertex().
00660  ***************************************************************************/
00661 private:
00664 struct computeCavityStruct {
00665   const IncrementalDelaunay& del;
00666 
00667   Cavity& cavity;
00668   const Point& p;
00669   computeCavityStruct(const IncrementalDelaunay& del, Cavity& cavity, const Point& p)
00670     : del(del), cavity(cavity), p(p) { }
00671   bool canEnter(const Simplex& s) const {
00672     bool b = del.inSphereAffine(s, p);
00673     dprintf("Considering %s: %s\n", s.toString().c_str(), b ? "insphere" : "out");
00674     return b;
00675   }
00676   bool choose(const Simplex& s) { cavity.add(s); return false; }
00677 };
00678 
00683 struct ExtendedCavityStruct {
00684   const IncrementalDelaunay& del;
00685   Cavity& cavity;
00686   vector<Simplex>& neighbours;
00687   const Point& p;
00688   ExtendedCavityStruct(const IncrementalDelaunay& dd, Cavity& cav, vector<Simplex>& neigh, const Point& pp)
00689     : del(dd), cavity(cav), neighbours(neigh), p(pp) { }
00690   bool canEnter(const Simplex& s) {
00691     bool b = del.inSphereAffine(s, p);
00692     if (b) {
00693       cavity.add(s);
00694     } else {
00695       neighbours.push_back(s);
00696     }
00697     dprintf("Considering %s: %s\n", s.toString().c_str(), b ? "insphere" : "out");
00698     return b;
00699   }
00700   static bool choose(const Simplex&) { return false; }
00701 };
00702 public:
00703 
00711 void computeCavity(const Simplex& s, const Point& p, Cavity& cavity, bool append = false) const {
00712 #ifndef NDEBUG
00713   if (!inSphereAffine(s, p)) {
00714     fprintf(stderr,
00715         "error: %s is not in the circumsphere of %s\n",
00716         p.toString().c_str(), s.toString().c_str());
00717     if (circumcenter(s).radius2 < 1e-8) {
00718       fprintf(stderr, "error: maybe there's too much refinement here?\n");
00719     }
00720     assert(inSphereAffine(s, p));
00721   }
00722 #endif
00723   if (!append) { assert(cavity.empty()); }
00724   computeCavityStruct dat(*this, cavity, p);
00725   complex_.dfsBySimplex(dat, s);
00726 }
00727 
00739 void computeCavity(const Simplex& s, const Point& p, Cavity& cavity, vector<Simplex>& neighbours,
00740     bool append = false) const {
00741   assert(inSphereAffine(s, p));
00742   if (!append) {
00743     assert(cavity.empty());
00744     assert(neighbours.empty());
00745   }
00746   ExtendedCavityStruct dat(*this, cavity, neighbours, p);
00747   complex_.dfsBySimplex(dat, s);
00748 }
00749 
00750 /***************************************************************************
00751  * Compute the star of a cavity.
00752  *
00753  * This gives us all the simplices that would come about from actually
00754  * inserting the point.  We can use this to check for slivers.
00755  ***************************************************************************/
00756 public:
00757 void computeStar(const Cavity& cavity, Vertex *v, Star& star) {
00758   complex_.computeStar(cavity, v, star);
00759 }
00760 
00761 /***************************************************************************
00762  * Insert a vertex.
00763  *
00764  * Assumes the cavity and star have been computed.
00765  ***************************************************************************/
00766 public:
00767 void insertVertex(const Cavity& cavity, Star& star) {
00768   //dumpEPS();
00769   complex_.replaceCavity(cavity, star);
00770   //dumpEPS();
00771 #ifdef EXPENSIVE_ASSERTIONS
00772   assert(isOriented());
00773   assert(isDelaunay());
00774 #endif
00775 }
00776 
00777 /***************************************************************************
00778  * Insert a vertex.
00779  *
00780  * Wrapper around the above functions.
00781  ***************************************************************************/
00782 public:
00783 void insertVertex(const Simplex& s, Vertex *v) {
00784   Cavity cavity;
00785   computeCavity(s, v->toPoint(), cavity);
00786   Star star;
00787   computeStar(cavity, v, star);
00788   insertVertex(cavity, star);
00789 }
00790 
00791 /***************************************************************************
00792  * Insert a vertex.
00793  *
00794  * Wrapper including expensive point location to initially locate v
00795  ***************************************************************************/
00796  public:
00797 void expensiveInsertVertex(Vertex *v) {
00798   Simplex s = expensiveLocate(v);
00799   insertVertex(s,v);
00800 }
00801 
00802 /***************************************************************************
00803  * \name Point Location
00804  * Find an Simplex that v encroaches, to kickstart insertion
00805 
00806  * Runs in <b>linear time</b>, hence the 'expensive' prefix.
00807  * TODO: improve to the expected n^(1/(d+1)) ?
00808  ***************************************************************************/
00810 private:
00811 struct expensiveLocateStruct {
00812   const IncrementalDelaunay& id_;
00813   const Point& p_;
00814   expensiveLocateStruct(const IncrementalDelaunay& id, const Point& p)
00815     : id_(id), p_(p) { }
00816   static bool canEnter(const Simplex&) { return true; }
00817   bool choose(const Simplex& s) {
00818     return id_.inSphereAffine(s, p_);
00819   }
00820 };
00821 
00822 public:
00823 Simplex expensiveLocate(const Point& p) const {
00824   expensiveLocateStruct loc(*this, p);
00825   return * complex_.findSimplex(loc, getHandle());
00826 }
00827 Simplex expensiveLocate(const Vertex *v) const {
00828   return expensiveLocate(v->toPoint());
00829 }
00831 
00832 /***************************************************************************
00833  * \name Complex access
00834  * Get at the underlying complex.  See the SimplicialComplex class for details.
00835  ***************************************************************************/
00837 public:
00838 const Complex& getComplex() const { return complex_; }
00839 Complex& getComplexRW() { return complex_; }
00840 const Simplex& getHandle() const { return complex_.getHandle(); }
00842 
00843 /***************************************************************************
00844  * \name Debugging routines.
00845  ***************************************************************************/
00847 private:
00848 struct verifyOrientedStruct {
00849   const IncrementalDelaunay& id_;
00850   verifyOrientedStruct (const IncrementalDelaunay& id): id_(id) { }
00851   static bool canEnter(const Simplex&) { return true; }
00852   bool choose (const Simplex& s) {
00853     if (!id_.isOriented(s)) {
00854       dprintf("unoriented simplex: %s\n", s.toString().c_str());
00855     }
00856     return ! id_.isOriented(s);
00857   }
00858 };
00859 struct verifyDelaunayStruct {
00860   const IncrementalDelaunay& id_;
00861   verifyDelaunayStruct (const IncrementalDelaunay& id): id_(id) { }
00862   static bool canEnter(const Simplex&) { return true; }
00863   bool choose (const Simplex& s) {
00864     OSimplex oriented = id_.getComplex().orientSimplex(s);
00865     for (unsigned i = 0; i <= topological; ++i) {
00866       OSimplex neighbour = id_.getComplex().getNeighbour(topological, oriented);
00867       if (!neighbour.isInternal()) continue;
00868       Vertex *v = neighbour[topological];
00869       if (id_.inSphereAffine(s, v->toPoint())) {
00870         dprintf("non-delaunay simplex: %s vs %s\n",
00871             s.toString().c_str(), v->toString().c_str());
00872         return true;
00873       }
00874     }
00875     return false;
00876   }
00877 };
00878 public:
00879 
00881 bool isOriented() const
00882 {
00883   if (topological == 1) return true;
00884 
00885   verifyOrientedStruct vOS(*this);
00886   return !getComplex().dfsBySimplex(vOS, getHandle());
00887 }
00888 
00890 bool isDelaunay() const
00891 {
00892   if (topological == 1) return true;
00893 
00894   verifyDelaunayStruct vDS(*this);
00895   return !getComplex().dfsBySimplex(vDS, getHandle());
00896 }
00898 
00899 
00900 /***************************************************************************
00901  * The constructor routines create vertices; these need to be cleaned
00902  * up eventually, else we leak.  The destructor goes through myVertices
00903  * and deletes them all.  This is how we add them to that array.
00904  * Other callers are free to do this also.
00905  ***************************************************************************/
00906 public:
00907 void giveVertexOwnership(Vertex *v) { myVertices_.push_back(v); }
00908 
00909 
00910 
00911 /***************************************************************************
00912  * \name Fields remaining to be commented
00913  *
00914  * These functions implement the more complicated constructors.
00915  * \see Constructors.
00916  * TODO: document.
00917  ***************************************************************************/
00919 public:
00926 static IncrementalDelaunay *
00927 computeBoundingSimplex(const BoundingBox& box, const ProjectionPlane& plane, double buffer, vector<Vertex*>& corners)
00928 {
00929   /* We want a simplex that fits the big box. */
00930   BoundingBox bigbox = box.scale(buffer);
00931 
00932   double temp_buffer = 5.0 * ambient;
00933   double simplexsidelength = (ambient + 2) * temp_buffer * bigbox.extent;
00934 
00935   PPoint origin = bigbox.mincorner - (PPoint(temp_buffer * bigbox.extent));
00936 
00937   for (unsigned i = 0; i < topological; i++) {
00938     PPoint corner = origin;
00939     corner[i] += simplexsidelength;
00940     corners.push_back(new Vertex(plane.unproject(corner)));
00941   }
00942   corners.push_back(new Vertex(plane.unproject(origin)));
00943   assert(corners.size() == topological + 1);
00944 
00945   return fromSimplex(corners, plane);
00946 }
00947 
00948 template <class p_iterator>
00949 static IncrementalDelaunay *
00950 computeBoundingSimplex(p_iterator begin, p_iterator end, double buffer, vector<Vertex*>& corners)
00951 {
00952   ProjectionPlane plane (begin, end);
00953   BoundingBox box (begin, end, plane);
00954 
00955   return computeBoundingSimplex(box, plane, buffer, corners);
00956 }
00957 
00958 static vector<Vertex*>
00959 unprojectBoundingSimplex(
00960     const vector<PPoint>& points,
00961     const ProjectionPlane& plane)
00962 {
00963   vector<Vertex*> verts;
00964   BOOST_FOREACH(const PPoint& pp, points) {
00965     verts.push_back(new Vertex(plane.unproject(pp)));
00966   }
00967   return verts;
00968 }
00969 
00970 static IncrementalDelaunay *makeBoundingSimplex(
00971     const ProjectionPlane& plane,
00972     const BoundingBox& box,
00973     double buffer) {
00974   vector<Vertex*> verts;
00975   IncrementalDelaunay *del = computeBoundingSimplex(box, plane, buffer, verts);
00976   del->myVertices_ = verts;
00977   return del;
00978 }
00979 
00980 template <class point_iterator>
00981 static IncrementalDelaunay *makeBoundingSimplex(point_iterator begin, point_iterator end, double buffer) {
00982   ProjectionPlane plane(begin, end);
00983   BoundingBox box(begin, end, plane);
00984   return makeBoundingSimplex(plane, box, buffer);
00985 }
00986 
00987 /***************************************************************************
00988  * Building a bounding box.
00989  *
00990  * Step 1: compute the plane.
00991  * Step 2: compute the min corner and the extent.
00992  * Step 3: build a very large bounding simplex.
00993  * Step 4: insert the bbox points into the simplex.
00994  ***************************************************************************/
00995 public:
00996 /* ppd is the desired points per dimension */
00997 /* buffer is how much whitespace around the points as a multiple of their width */
00998 template <class p_iterator> static IncrementalDelaunay *createOrthoBox(
00999     const p_iterator& begin,
01000     const p_iterator& end,
01001     unsigned ppd, double buffer) {
01002   assert(ppd >= 2);
01003   assert(buffer >= 0.0);
01004 
01005   dprintf("Compute plane");
01006   ProjectionPlane pplane (begin, end);
01007 
01008   dprintf("Compute extent");
01009   BoundingBox box = BoundingBox(begin, end, pplane);
01010   assert (box.extent > 0.0);
01011 
01012   BoundingBox bigbox = box.scale(buffer);
01013   assert (bigbox.extent > 0.0);
01014   printf("Meshing domain: extent %g (buffer %g), %u points per side \n",
01015       bigbox.extent, buffer, ppd);
01016 
01017   /*********************************************************
01018    * Construct a Large Right Simplex to Bound the Hypercube
01019    *********************************************************/
01020   // Multiples of 'extent' are used as a buffer between the cube and the right
01021   // bounding simplex.  Make sure this is big enough to guarantee that the cube sides are Delaunay.
01022   // THIS ONLY EXISTS INTERMEDIATELY
01023   // We must rememeber to delete the corners.
01024   vector<Vertex*> vcorners;
01025   IncrementalDelaunay *boundingcomplex = computeBoundingSimplex(box, pplane, buffer, vcorners);
01026 
01027   /*********************************************************
01028    * Insert the corners of the hypercube
01029    *********************************************************/
01030 
01031   // Really this is a Point in Z^topo
01032   dprintf("Insert cube");
01033   PPoint coordinates(0.0);
01034   double surface_grid_unit = bigbox.extent / (ppd - 1.0);
01035   size_t surface_size = sizeOfSurface(ppd);
01036 
01037   for (size_t i = 0; i < surface_size; i++) {
01038     PPoint surface_point = (bigbox.mincorner + (coordinates * surface_grid_unit));
01039     Vertex *surface_v = new Vertex(pplane.unproject(surface_point));
01040     dprintf("Inserting %s", surface_v->toString().c_str());
01041     if (topological != ambient) {
01042       dprintf(" (projects to %s)\n", surface_point.toString().c_str());
01043     } else { dprintf(""); }
01044     // expensive insert since everything is constant size
01045     boundingcomplex->giveVertexOwnership(surface_v);
01046     boundingcomplex->expensiveInsertVertex(surface_v);
01047     coordinates = incrementSurface(coordinates,ppd);
01048   }
01049   dprintf("got hypercube");
01050 
01051   /*********************************************************
01052    * Remove the Large Bounding Simplex
01053    *********************************************************/
01054   boundingcomplex->removeCorners(vcorners);
01055   BOOST_FOREACH(Vertex *v, vcorners) { delete v; }
01056 
01057   return boundingcomplex;
01058 }
01059 
01060 
01066 template <class p_iterator> static IncrementalDelaunay *createMeshBox(
01067     const p_iterator& begin,
01068     const p_iterator& end,
01069     double rho) {
01070 
01071   // This routine just computes the following two numbers:
01072   unsigned k;  // number of points per side.
01073   double D;    // buffer size
01074 
01075   // This term keeps coming up.
01076   const double gamma =
01077     0.5 * sqrt(topological - 1) * (1.0 + sqrt(2.0) * rho / (rho - 1.0));
01078 
01079   // See the proof in doc/proofs.  We want to choose k, an integer, that
01080   // is larger than 1 + 2 gamma.  But, the closer it gets to be exactly 1 + 2
01081   // gamma, the larger D gets.  Therefore, I set k to be the ceiling of
01082   // 1 + 2 gamma + epsilon.
01083   const double epsilon = 0.5;
01084 
01085   k = 1 + unsigned (ceil(2.0 * gamma + epsilon));
01086   D = gamma / (double(k) - 1.0 - 2.0 * gamma) + 1e-8;
01087 
01088   return IncrementalDelaunay::createOrthoBox(
01089         begin, end, k, D);
01090 }
01091 
01092 // helper for createOrtho to loop through d-dimensional points on hypercube surface
01093 //this just iterates over {ppd}^topo
01094 private:
01095 static PPoint incrementGrid(PPoint coord, unsigned ppd) {
01096   unsigned carry_bit = 0;
01097   coord[0]++; //dunno if shorter works
01098   while (coord[carry_bit] >= ((double)ppd-0.5))
01099   {
01100     coord[carry_bit] = 0;
01101     carry_bit++;
01102     if (carry_bit >= topological) return PPoint(0.0); // overflow wrap
01103     coord[carry_bit]++;
01104   }
01105   return coord;
01106 }
01107 
01108 // is a point in {ppd}^topo on the surface
01109 private:
01110 static bool onSurface(PPoint coord, unsigned ppd) {
01111   for (unsigned j = 0; j < topological; j++) {
01112     if (coord[j] == 0) return true;
01113     if (coord[j] == ppd - 1) return true;
01114   }
01115   return false;
01116 }
01117 
01118 // iterate over {ppd}^topo, stopping only on surface
01119 private:
01120 static PPoint incrementSurface(PPoint coord, unsigned ppd) {
01121   coord = incrementGrid(coord,ppd);
01122   while (!onSurface(coord,ppd))
01123     coord = incrementGrid(coord,ppd);
01124   return coord;
01125 }
01126 
01127 private:
01128 static unsigned ipow(unsigned base, unsigned exp) {
01129   unsigned x = 1;
01130   for(unsigned i = 0; i < exp; i++) {
01131     x *= base;
01132   }
01133   return x;
01134 }
01135 
01136 private:
01137 static unsigned sizeOfSurface(unsigned ppd) {
01138   return ipow(ppd, topological) - ipow(ppd-2, topological);
01139 }
01140 
01141 // helpers for step 4
01142 
01143 private:
01144 struct CollectCorners {
01145   vector<Simplex> found_;
01146   const vector<Vertex*>& corners_;
01147   CollectCorners(const vector<Vertex*>& corners) : corners_(corners) {};
01148   static bool canEnter(const Simplex&) { return true; }
01149   bool choose(const Simplex& s) {
01150     BOOST_FOREACH(const Vertex *v, corners_) {
01151       if (s.has(v)) { found_.push_back(s); return false; }
01152     }
01153     return false;
01154   }
01155 };
01156 struct FindNotCorner {
01157   const vector<Vertex*>& corners_;
01158   FindNotCorner(const vector<Vertex*>& corners) : corners_(corners) {};
01159   static bool canEnter(const Simplex&) { return true; }
01160   bool choose(const Simplex& s) const {
01161     BOOST_FOREACH(const Vertex *v, corners_) {
01162       if (s.has(v)) { return false; }
01163     }
01164     return true;
01165   }
01166 };
01167 
01168 
01169 private:
01170 void removeCorners(const vector<Vertex*>& corners) {
01171   CollectCorners kf(corners);
01172   getComplex().dfsBySimplex(kf, getHandle());
01173 
01174   FindNotCorner nc(corners);
01175   Simplex safeHandle = *getComplex().findSimplex(nc, getHandle());
01176 
01177   // iterate over the corner simplices we found and remove them
01178   dprintf("found %u to kill\n", (unsigned)kf.found_.size());
01179   dprintf("new handle %s\n", safeHandle.toString().c_str());
01180   complex_.setHandle(safeHandle);
01181   BOOST_FOREACH(const Simplex& s, kf.found_) {
01182     complex_.checkedRemoveSimplex(s);
01183   }
01184 }
01185 
01186 
01187 public:
01188 template <class p_iterator> static IncrementalDelaunay *createSquare(
01189     const p_iterator begin, const p_iterator& end,
01190     double buffer = 1.0) {
01191   return createOrthoBox(begin, end, 2, buffer);
01192 }
01193 template <class p_iterator> static IncrementalDelaunay *create4box(
01194     const p_iterator begin, const p_iterator& end) {
01195   return createOrthoBox(begin, end, 4, 1.500001);
01196 }
01198 
01199 /***************************************************************************
01200  * Computes the Delaunay triangulation -- slowly, for some inputs.
0