delaunay() - why does it not work for coplanar points in 3-D?

Why exactly is it a condition that the points not be co-planar?

I am trying on a set and getting the following error:

Error computing the Delaunay triangulation. The points may be coplanar or collinear.

I am pretty sure I am inputting in the correct format as specified by the documentation.

Browsing the wikipedia, there seems to be no mention of this requirement for a triangulation to exist:

https://en.wikipedia.org/wiki/Delaunay_triangulation

I am not sure if this is just an arbitrary restriction imposed by MATLAB or if I am missing something...

Answer 1

[John D'Errico]

Why should it work?

For example, we can try this:

XYZ1 = dec2bin(0:7) - '0'

XYZ1 = 8×3

0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 1

So we have the 8 vertices of a unit cube. There are subsets of these points that are coplanar. Of course, delaunayTriangulation will have no problem with the set.

T1 = delaunayTriangulation(XYZ1)

T1 = 
  delaunayTriangulation with properties:

              Points: [8×3 double]
    ConnectivityList: [6×4 double]
         Constraints: []

Howver, suppose we try a simple set of points that lie entirely in a plane?

XYZ2 = randn(8,2)*randn(2,3)

XYZ2 = 8×3

1.4814 -2.7935 3.5939 0.0841 -3.3857 1.4239 -0.3402 0.9955 -0.9591 -0.7697 3.5052 -2.6436 0.0657 0.3543 -0.0214 -0.4894 0.3221 -0.9601 0.3676 -0.4206 0.7887 0.7759 -3.1797 2.5313

These points all lie in a single plane. We can see that is true, by looking at the results of the singular value decomposition. It will have one zero element in S.

format long g
svd(XYZ2)

ans = 3×1

8.61903296375816 1.86223136243366 5.13026536773492e-17

plot3(XYZ2(:,1),XYZ2(:,2),XYZ2(:,3),'o')
box on
grid on

I should rotate it around to show the points lie exactly in a plane, but I'm too lazy now.

T2 = delaunayTriangulation(XYZ2)

T2 = 
  delaunayTriangulation with properties:

              Points: [8×3 double]
    ConnectivityList: [11×4 double]
         Constraints: []

See that delaunayTriangulation survives, although off those simplexes will have zero volume.

However, some of the older tools will fail.

delaunayn(XYZ2)

Error using matlab.internal.math.qhull
qhull precision warning:
The initial hull is narrow (cosine of min. angle is 1.0000000000000000).
A coplanar point may lead to a wide facet. Options 'QbB' (scale to unit box)
or 'Qbb' (scale last coordinate) may remove this warning. Use 'Pp' to skip
this warning. See 'Limitations' in qh-impre.htm.
QH6114 qhull precision error: initial simplex is not convex. Distance=-6.7e-16

While executing: | qhull d Qt Qbb Qc
Options selected for Qhull 2010.1 2010/01/14:
run-id 958828864 delaunay Qtriangulate Qbbound-last Qcoplanar-keep
_pre-merge _zero-centrum Qinterior-keep Pgood _max-width 6.9
Error-roundoff 1.4e-14 _one-merge 1.3e-13 Visible-distance 8.3e-14
U-coplanar-distance 8.3e-14 Width-outside 1.7e-13 _wide-facet 5e-13
_narrow-hull 0

precision problems (corrected unless 'Q0' or an error)
1 flipped facets
1 zero divisors during gaussian elimination

The input to qhull appears to be less than 4 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet. The maximum round off error for
computing distances is 1.4e-14. The center point, facets and distances
to the center point are as follows:


facet p7 p1 p0 p3 distance= -1.3e-15
facet p2 p1 p0 p3 distance= -2.2e-16
facet p2 p7 p0 p3 distance= -7e-18
facet p2 p7 p1 p3 distance= -7e-17
facet p2 p7 p1 p0 distance= -2.3e-16

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates. Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
0: -0.7697 1.481 difference= 2.251
1: -3.386 3.505 difference= 6.891
2: -2.644 3.594 difference= 6.238
3: -6.939e-18 6.891 difference= 6.891

If the input should be full dimensional, you have several options that
may determine an initial simplex:
- use 'QJ' to joggle the input and make it full dimensional
- use 'QbB' to scale the points to the unit cube
- use 'QR0' to randomly rotate the input for different maximum points
- use 'Qs' to search all points for the initial simplex
- use 'En' to specify a maximum roundoff error less than 1.4e-14.
- trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
- use 'QJ' to joggle the input and make it full dimensional
- use 'Qbk:0Bk:0' to delete coordinate k from the input. You should
pick the coordinate with the least range. The hull will have the
correct topology.
- determine the flat containing the points, rotate the points
into a coordinate plane, and delete the other coordinates.
- add one or more points to make the input full dimensional.


This is a Delaunay triangulation and the input is co-circular or co-spherical:
- use 'Qz' to add a point "at infinity" (i.e., above the paraboloid)
- or use 'QJ' to joggle the input and avoid co-circular data


Error in delaunayn (line 101)
t = matlab.internal.math.qhull(x', 'd ', opt);

Effectively, the problem probably arises from the older codes, which were dependent on qhull libraries. It seems that tools like delaunayTriangulation survive. Note that a common solution in those older codes was to use a joggle, so a tiny additive random perturbation, that would result in the points no longer being coplanar.