Scipy Delaunay triangulation is not oriented

[patricksnape]

PC details

• scipy==0.12.0
• Linux dana 3.8.0-30-generic #44-Ubuntu SMP Thu Aug 22 20:52:24 UTC 2013 x86_64 x86_64 x86_64 GNU/Linux
• Intel(R) Xeon(R) CPU E5-1650 0 @ 3.20GHz
• Linked against Intel MKL (though issue replicated on a machine linked against Atlas)

Issue

With this mesh I am seeing inconsistent results from scipy.spatial.Delaunay. Assuming that the mesh has been loaded using numpy.loadtxt into a variable called points, I duplicate the issue as follows:

from scipy.spatial import Delaunay
scipy_t = Delaunay(points).simplices

I have some code that checks the consistency of a triangulation (in terms of ordering) and the above code generates an inconsistent triangulation.

Ground Truth

In order to sanity check I compared the triangulation against the one created by pyhull. This is created using the same point set as follows:

from pyhull.delaunay import DelaunayTri
pyhull_t = np.asarray(DelaunayTri(points).vertices)

I've also generated the triangulation in Matlab using both delaunay and delaunayn and it yields results identical to pyhull_t. My triangulation consistency verifier also confirms that the pyhull_t triangulation is correct.

Options

In order to try and ensure that scipy.spatial.Delaunay is applying the same algorithm as pyhull, I dove in to their source code. There I noticed that they always set the following flags:

i Qt Qc Qbb

but running with those flags as follows:

from scipy.spatial import Delaunay
scipy_t = Delaunay(points, qhull_options='Qbb Qc').simplices

still does not yield correct results.

Example output

An example of the triangulation I am seeing using the flags Qbb Qc:

print scipy_t

[[16549 35521 21080]
 [58451 78938 91348]
 [38031 32638 33401]
 ..., 
 [47968 47967 48281]
 [47968 47656 47969]
 [47968 47655 47656]]

print pyhull_t

[[35521 16549 21080]
 [78938 58451 91348]
 [32638 38031 33401]
 ..., 
 [47968 47967 48281]
 [47656 47968 47969]
 [47655 47968 47656]]

where we can see that the same indices are being generated, but the first two are sometimes swapped.

[pv]

By "consistent", do you mean counterclockwise ordering of the neighbor triangles? Such ordering is indeed not guaranteed.

The Scipy code works essentially in the Fv Fn mode of qdelaunay, so it doesn't do the postprocessing used in the i mode. As noted in the documentation, the (only) guarantee is that the k-th neighbor is opposite to the k-th vertex.

[pv]

Based on this information, I'm marking this as an enhancement request. Producing oriented output is also possible, and probably would be the least surprising result in 2D.

[jabooth]

For 2D surface triangulation consistent orientation of the triangles would be the behaviour I would expect by default. From casual experiment it seems this is the default behaviour in Matlab, for instance.

[patricksnape]

Under what conditions would we expect the triangulation to default to oriented? 2D and 3D delaunay? I'm happy to open a pull request for this as I'm pretty sure I've identified where in the code to fix this. Maybe we would like to be able to handle the i flag?

[pv]

qdelaunay man page seems to claim it produces oriented output for N-d. Fixing the orientation in the same way is quite cheap --- the only thing that needs to be taken care of is doing the same permutation in the neighbor simplex map too.

Due to the fact that we are using qhull as a library rather than command-line program, fixing up the orientation is slightly more involved than just passing on flags. The relevant Qhull code is in qhull/src/io.c:1178

[patricksnape]

Yes the relevant code seems to be Line 1175:

case qh_PRINTincidences:
case qh_PRINToff:
case qh_PRINTtriangles:
    if (qh hull_dim == 3 && format != qh_PRINTtriangles)
      qh_printfacet3vertex(fp, facet, format);
    else if (facet->simplicial || qh hull_dim == 2 || format == qh_PRINToff)
      qh_printfacetNvertex_simplicial(fp, facet, format);
    else
      qh_printfacetNvertex_nonsimplicial(fp, facet, qh printoutvar++, format);
    break;

as you stated.

Now looking at this it alludes to being able to handle ND simplicial facets. Drilling in to that code we find:

if ((facet->toporient ^ qh_ORIENTclock) || (qh hull_dim > 2 && !facet->simplicial)) {
    FOREACHvertex_(facet->vertices)
      qh_fprintf(fp, 9130, "%d ", qh_pointid(vertex->point));
  }else {
    FOREACHvertexreverse12_(facet->vertices)
      qh_fprintf(fp, 9131, "%d ", qh_pointid(vertex->point));
  }

assuming we ignore the clockwise orientation flag then it only actually flips the vertices if the mesh is 2D and simplicial. Otherwise it just prints them in order.

Therefore, it appears it only actually corrects 2D meshes - unless I am missing somewhere else that it does reorientation?

Thus, a proposed solution would be something like (Around Line 609):

# Orient the returned facets by default
# If we are two dimensional then swap the first and second
# vertex (a facet size of 3 means 2D)
if (self._is_delaunay and not facet.toporient
    and facet_ndim == 3 and facet.simplicial):
    temp_p = facet.vertices.e[0].p
    temp_i = facet.vertices.e[0].i
    facet.vertices.e[0].p = facet.vertices.e[1].p
    facet.vertices.e[0].i = facet.vertices.e[1].i
    facet.vertices.e[1].p = temp_p
    facet.vertices.e[1].i = temp_i

[pv]

Pretty much so. However, it's best to not modify the qhull data structures directly as they may be needed later on for the incremental mode.

I'm not sure about the > 2D case, but note that it flips the first two elements for all dimensions if the per-facet toporient flag is set.

[patricksnape]

I'm not sure about the > 2D case, but note that it flips the first two elements for all dimensions if the toporient flag is set.

I could just be tired but are you sure this is true? It looks to be (from the snippet I posted) that if toporient is True then it NEVER flips?

[patricksnape]

What about adding a kwarg ordered or the like that raises an exception if the points are not 2D? Then, if True, it adds a new property to the Delaunay object ordered_simplices. This maintains the state of the Delaunay object and won't upset anyone who relies on the current implementation.

[pv]

Ah OK, you're right, it flips only 2D.

I think any code that relies on a specific order of the vertices currently, is incorrect, because the order is not guaranteed and depends on implementation details. So we are free to change the behavior also without keywords.

Always oriented in 2D is probably a good convention, and there's no performance cost.

[jabooth]

So we are free to change the behaviour also without keywords.

Is this a little dangerous? If we fix the simplices in our Delaunay object as we pull the data out of the underlying QHull object without updating the QHull state, could this have nasty ramifications for the incremental mode, or are we able to reason that this is not a problem?

[pv]

The information in simplices never goes back to Qhull, so I don't see immediate problems here.

[patricksnape]

OK. My intuition says that in order to maintain a consistent neighbours ordering we also need to flip them?

[pv]

Yes, the same orientation correction needs to be done also to that array.

[pv]

Fixed in gh-2954

[WombatBill]

I'm still not sure what to do with the "simplices" that result from 3d triangulation; how to convert them to the coordinates for individual triangular facets.

[tylerjereddy]

@WombatBill

The simplices from a 3D Delaunay triangulation are tetrahedra--so you'll have 4 triangular facets per simplex. There should be a variety of ways to iterate through the triangular facets, but that's probably better suited for a usage site like StackOverflow than a closed enhancement request on GitHub.

There's a concise post / visualization available on the topic from Joseph O'Rourke, a well-known Professor of Computational Geometry & author of related textbooks.

[WombatBill]

Which... makes me wonder what the purpose is... It's "trianugulation," not "tetraherafication," or whatever...

…

-----Original-Nachricht----- The simplices from a 3D Delaunay triangulation are tetrahedra--so you'll have 4 triangular facets per simplex. 

[tylerjereddy]

Triangulation can be defined with no reference to triangles whatsoever and you can generalize the concept to N-dimensions.

From: Devadoss, Satyan L.; O'Rourke, Joseph. Discrete and Computational Geometry (2011) Chapter 3: Triangulations:

If you're no longer in the plane, the maximal set of non-crossing edges isn't either. There are a wide variety of reasons one may want this information. You might sum together the volumes of the tetrahedra to get the volume of the convex hull, or you may be designing a network that requires the connectivity enforced by the tetrahedra, etc.

The circumspheres that include the vertices of the simplices (tetrahedra) also have unique properties.

[WombatBill]

Well, I want the coordinates of the triangles as a result of Delaunay's triangulation process, so I can analyse those statistically in order to see whether the point distribution is clustered, random or patterned. And I'm doing this in 3D.

…

-----Original-Nachricht----- Triangulation can be defined with no reference to triangles whatsoever and you can generalize the concept to N-dimensions.