Infinite initial bounding triangle in iterative Delaunay triangulators

Question

Most iterative algorithms require an initial empty triangle to get the ball rolling. It seems like a commonly used trick is just to make the super triangle very large in comparison with the point set.

But according to "Numerical recipes: the art of scientific computing":

"...if the distance is merely finite (to the boundary points) the constructed triangulation may be not-quite Delaunay. For example its outer boundary could in unusual cases be left slightly concave, with small negative angles on the order of the diameter of the "real" point set divided by the distance to the "fictitious" (boundary) points.

So what options are there for augmenting Cartesian coordinates with points at infinity, without having to convert all the input to a different coordinate system such as homogeneous coordinates? How do these points fit in with the usual geometric predicates CCW and Incircle?

Incircle (a,b,c) Infinity -> False. provided that a,b,c are finite.

But what about when one of a,b,c is a point at infinity? Does the circle become a half-plane and the test then become a CCW check? What if 2 or more points on the circumcircle are infinite? does the circle expand into a full plane causing the test to always yield true? What about CCW? how do you classify a point in relation to a line that has one or more points at infinity?

Answer 1

It is quite easy to implement a Delaunay triangulation by adding a point at infinity. Chose a convention for the infinite vertex (e.g., vertex index -1).

At the beginning, you create an initial finite tetrahedron T0 between 4 non-coplanar points taken from the input pointset, and then create four infinite tetrahedra that connect each face of T0 to the infinite vertex 0 (and do not forget to interconnect them properly along their common infinite faces).

Then you insert each point p, one by one, as usual (Boyer-Watson algo), by (1) computing the tetrahedron T that contains the point p (locate) and (2) finding the conflict zone, as follows:

(1) locate is unmodified: you walk to p from a random finite tetrahedron. During walking, if you arrive in an infinite tetrahedron, then you stop there (this means the point is outside the convex hull of the previously inserted points)

(2) Determining the conflict zone needs a small modification :

• A point p is in conflict with a finite tetrahedron T if its in its circumscribed sphere (as usual)

• For an infinite tetrahedron T = (p1, p2, p3, p4), one of p1,p2,p3,p4 is the infinite vertex (for instance p2, then T = (p1, infinite, p3, p4)). To check whether p is in conflict with T, replace the infinite vertex with p, and compute the signed volume of the tetrahedron: signed_volume(p1, p, p3, p4) in our example, if it is positive, then T is in conflict with p. If it is negative then T is not in conflict with p. If p is exactly on the supporting plane of the three finite vertices of T, then T is in conflict with p if T' is in conflict with p, where T' is the tetrahedron adjacent to T along the finite face of T (equivalently, instead of quering T', one may check whether p is in the circumscribed circle of the finite facet of T).

See implementations in:

• geogram: http://alice.loria.fr/software/geogram/doc/html/classGEO_1_1Delaunay3d.html
• CGAL: http://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation_3Delaunay

Answer 2

I think you are making life difficult for yourself by trying to define the complete mathematics of geometries that contain points at infinite distance away. This isn't necessary for the purposes of your original problem of calculating the Delaunay triangulation accurately.

I wrote a Delaunay generator in Java a while ago and it is available here: http://open.trickl.com/trickl-graph/index.html

See http://open.trickl.com/trickl-graph/apidocs/index.html and in particular:

    // Check if fourth point is within the circumcircle defined by the first three
   private boolean isWithinCircumcircle(PlanarGraph<V, E> graph, V first,
           V second,
           V third,
           V fourth) {
      // Treat the boundary as if infinitely far away
      Coordinate p = vertexToCoordinate.get(fourth);
      if (PlanarGraphs.isVertexBoundary(graph, first)) {
         return isLeftOf(third, second, p);
      } else if (PlanarGraphs.isVertexBoundary(graph, second)) {
         return isLeftOf(first, third, p);
      } else if (PlanarGraphs.isVertexBoundary(graph, third)) {
         return isLeftOf(second, first, p);
      } else if (PlanarGraphs.isVertexBoundary(graph, fourth)) {
         return false;
      }

      Coordinate a = vertexToCoordinate.get(first);
      Coordinate b = vertexToCoordinate.get(second);
      Coordinate c = vertexToCoordinate.get(third);

      boolean within = (a.x * a.x + a.y * a.y) * getDblOrientedTriangleArea(b, c, p)
              - (b.x * b.x + b.y * b.y) * getDblOrientedTriangleArea(a, c, p)
              + (c.x * c.x + c.y * c.y) * getDblOrientedTriangleArea(a, b, p)
              - (p.x * p.x + p.y * p.y) * getDblOrientedTriangleArea(a, b, c) > 0;

      return within;
   }

Here, the boundary points are explicitly checked for when checking the circumcircle condition, so they can effectively be treated as "infinitely" away. This is much easier than figuring out all the geometrical implications than you describe.