Circles.h

Go to the documentation of this file.

/* $Id$ */
#ifndef Circles_HEADER
#define Circles_HEADER

#include <boost/array.hpp>

#ifdef HAVE_LAPACK
#  include <clapack.h>
#endif

#include <math.h>

#include "Point.h"
#include "PointConvert.h"
#include "predicates.h"

namespace Geometry {
using namespace std;
using namespace boost;

template <size_t d> struct CenterRadius {
  Point<d> center;
  double radius2;
  CenterRadius(Point<d> c, double r2)
    : center(c), radius2(r2) { }
  CenterRadius()
    : center(0.0), radius2(0.0) { }

  double dist2(const CenterRadius& other) {
    return center.dist2(other.center);
  }

  bool approximatelyIntersects(const CenterRadius& other) {
    double d2 = dist2(other);
    // using branches to avoid the 3 sqrts here might be useful, depending
    // on architecture.
    double r = sqrt(radius2) + sqrt(other.radius2);
    return sqrt(d2) < r;
  }

  std::string toString() const {
    char buffer[4096];
    snprintf(buffer, sizeof(buffer), "%s // %g", center.toString().c_str(), radius2);
    return std::string(buffer);
  }
};

#ifdef HAVE_LAPACK
/***************************************************************************
 * Compute the circumcenter of the simplex.  If we live in a lower dimension
 * than the ambient, compute the center of the diametral ball.
 *
 * This is a general-dimensional algorithm, fairly accurate except on slivers,
 * but amazingly slow because of all the marshalling / unmarshalling code
 * talking to LAPACK.
 *
 * First, translate the simplex by q so that it includes the origin.
 *
 * We want a point x that is equidistant to all the points:
 *   (x-pi)^2 = (x-pj)^2  for all i, any j != i
 * Let pj be the origin; that yields:
 *   (x-pi)^2 = x^2       for all i
 * Algebra then gives the system:
 *   pi^T x = 0.5 pi^2    for all i
 * Rewrite as:
 *   P^T x = 0.5 b        where P has pi as its columns, b = {pi^2}
 *
 * If P is a lower-dimensional space, then we also want to impose that
 * x be on the affine plane of P (which goes through the origin):
 *   x = P alpha
 * Substitute in to see:
 *   P^T P alpha = 0.5 b
 * Solve for alpha, then substitute for x.
 *
 * Having found x, we translate back by q.
 ***************************************************************************/
template <unsigned ambient, unsigned topological>
CenterRadius<ambient>
circumcenter(const array<Point<ambient>, topological+1>& points) {
  typedef Point<topological> PPoint;
  typedef Point<ambient> Point;

  // Compute P.
  array<Point, topological> P;
  const Point& q = points[topological];
  for(unsigned i = 0; i < topological; i++) {
    P[i] = ( points[i] - q );
  }

  // Compute b.
  PPoint b;
  for(unsigned i = 0; i < topological; i++) {
    b[i] = P[i].norm2() / 2.0;
  }

  // Compute A, which is P^T or P^T P depending on dimension.
  // A is column-major.
  array<double, topological * topological> A;
  for(unsigned i = 0; i < topological; i++) {
    for(unsigned j = 0; j < topological; j++) {
      if (ambient == topological) {
        A[topological * j + i] = P[i][j];
      } else {
        A[topological * j + i] = P[i].dot(P[j]);
      }
    }
  }

  // Call LAPACK.  This is a bit ugly.
  {
    integer n = topological, nrhs = 1, lda = topological,
        pivots[topological], ldb = topological, info;
    dgesv_(&n, &nrhs, A.c_array(), &lda, pivots, b.c_array(), &ldb, &info);
    assert(info == 0);
  }

  // Now, A and b have been clobbered.  b is our answer.
  // We have to turn the answer into the circumcenter.
  Point c; // uninitialized
  if (ambient == topological) {
    c = pointCast<ambient,topological>(PPoint(b));
  } else {
    // b is now alpha; need to compute P alpha
    c.assign(0.0); // init to zero and sum up
    for (unsigned i = 0; i < topological; i++) {
      c += P[i] * b[i];
    }
  }

  // For basically free, we can get the radius along with the center, so
  // return it too.
  double r2 = c.norm2();

  // Don't forget to translate back by q.
  c += q;

  //printf("center %s, radius %g\n", c.toString().c_str(), sqrt(r2));
  return CenterRadius<ambient>(c, r2);
}
#endif

/***************************************************************************
 * typedef the CenterRadius and Point properly.
 * assumes we want 'CenterRadius' to mean CenterRadius<ambient>
 * Use a macro in case gcc decides to change how this is supposed to be
 * written.  I undef the macro at the end of the file.
 ***************************************************************************/
#define typedef_ambient_types(ambient) \
  typedef Geometry::Point< (ambient) > Point; \
  typedef Geometry::CenterRadius< (ambient) > CenterRadius /* no ; */

/***************************************************************************
  * Special-dimension circumcenter code: for segments.
  *
  * This is very cheap, very robust.
 ***************************************************************************/
template <unsigned ambient>
CenterRadius<ambient>
circumcenter_segment(const Point<ambient>& a, const Point<ambient>& b) {
  return CenterRadius<ambient>(
      (a + b) * 0.5, /* center is average of the points */
      (a - b).norm2() * 0.25 /* radius is |a-b|/2; we return radius^2 */);
}

/***************************************************************************
  * Special-dimension circumcenter code: for topological triangles.
  *
  * For ambient dimensions 2 and 3, I have special-dimensional code which is
  * much faster than calling the LAPACK-based circumcenter() routine above.
  * The special-dimensional code solves the same linear system using Cramer's
  * Rule.  For larger ambient dimensions, the code calls the LAPACK-based
  * routine (of course, HAVE_LAPACK must be defined then).
  *
  * I define cctr_tri_struct and specialize it because at the moment (2007),
  * C++ doesn't allow template specialization of functions.
  * Call circumcenter_triangle() and ignore the structures.
 ***************************************************************************/
template <unsigned ambient> struct cctr_tri_struct
#ifdef HAVE_LAPACK
{
  typedef_ambient_types(ambient);
  CenterRadius operator() (const Point& a, const Point& b, const Point& c)
  {
    array<Point, 3> points;
    points[0] = a;
    points[1] = b;
    points[2] = c;
    return circumcenter<ambient, 2>(points);
  }
}
#endif
;

template <> struct cctr_tri_struct<2> {
  static const unsigned ambient = 2;
  typedef_ambient_types(ambient);
  CenterRadius operator() (const Point& a, const Point& b, const Point& c)
  {
    /* Straight out of Shewchuk's robnotes */
    double denom = 2.0 * orient2d(a.data(), b.data(), c.data());
    Point r = a - c;
    Point s = b - c;
    double r2 = r.norm2();
    double s2 = s.norm2();
    double numx = r2 * s[1] - s2 * r[1];
    double numy = s2 * r[0] - r2 * s[0];
    Point q;
    q[0] = numx / denom;
    q[1] = numy / denom;
    return CenterRadius(c + q, q.norm2());
  }
};

template <> struct cctr_tri_struct<3> {
  static const unsigned ambient = 3;
  typedef_ambient_types(ambient);
  CenterRadius operator() (const Point& a, const Point& b, const Point& c)
  {
    /* Straight out of Shewchuk's robnotes */
    Point r = a - c;
    Point s = b - c;
    double r2 = r.norm2();
    double s2 = s.norm2();
    Point rs = crossProduct(r, s);
    Point r2s_s2r = (s * r2 - r * s2);
    Point numerator = crossProduct(r2s_s2r, rs);
    double denominator = 2.0 * rs.norm2();
    Point q = numerator / denominator;
    return CenterRadius(c + q, q.norm2());
  }
};

template <unsigned ambient>
CenterRadius<ambient>
circumcenter_triangle(const Point<ambient>& a, const Point<ambient>& b, const Point<ambient>& c) {
  return cctr_tri_struct<ambient>()(a,b,c);
}

/***************************************************************************
  * Special-dimension circumcenter code: for topological tetrahedra.
  *
  * The code in ambient dimension 3 is both faster, and on sliver is far more
  * numerically stable than the general-dimension code.
  * TODO: it is terrible on skinny tets -- this code is solving a linear system
  * Ax = b, where A is degenerate on slivers.
  *
  * For higher ambient dimension, HAVE_LAPACK must be defined.
  *
  * I define cctr_tet_struct and specialize it because at the moment (2007),
  * C++ doesn't allow template specialization of functions.
  * Call circumcenter_tetrahedron() and ignore the structures.
  *
 ***************************************************************************/
template <unsigned ambient> struct cctr_tet_struct
#ifdef HAVE_LAPACK
{
  typedef_ambient_types(ambient);
  CenterRadius operator() (const Point& a, const Point& b, const Point& c, const Point& d)
  {
    array<Point, 4> points;
    points[0] = a;
    points[1] = b;
    points[2] = c;
    points[3] = d;
    return circumcenter<ambient, 3>(points);
  }
}
#endif
;

template <> struct cctr_tet_struct<3> {
  static const unsigned ambient = 3;
  typedef_ambient_types(ambient);
  CenterRadius operator() (const Point& a, const Point& b, const Point& c, const Point& d)
  {
    /* Code from jrs' robnotes. */
    double denominator = 2.0 * orient3d(a.data(), b.data(), c.data(), d.data());

    /* Denominator is related to the volume; if the volume is small,
     * and we have a sliver, then all four points are almost cocircular,
     * so an arbitrary face's circumcenter is the global circumcenter.
     *
     * TODO: What if it's not a sliver?
     * TODO: What is small?  We've seen problems at 10^-15, but is 10^-10 too big?
     */
    if (fabs(denominator) < 1e-10) {
      return circumcenter_triangle<ambient>(a, b, c);
    }

    Point t = d - a; double t2 = t.norm2();
    Point u = d - b; double u2 = u.norm2();
    Point v = d - c; double v2 = v.norm2();

    Point uv = crossProduct(u,v);
    Point vt = crossProduct(v,t);
    Point tu = crossProduct(t,u);
    Point numerator = uv*t2 + vt*u2 + tu*v2;
    Point q = numerator / denominator;
    return CenterRadius(d + q, q.norm2());
  }
};

template <unsigned ambient>
CenterRadius<ambient>
circumcenter_tetrahedron(const Point<ambient>& a, const Point<ambient>& b, const Point<ambient>& c, const Point<ambient>& d) {
  return cctr_tet_struct<ambient>()(a,b,c,d);
}

#undef typedef_ambient_types

}; // namespace Geometry

#endif

---