Do not correct orientation of triangles returned by Qhull

ianthomas23 · ianthomas23 · commit 178d453df1ad · 2015-05-29T10:23:29.000+01:00
diff --git a/lib/matplotlib/tests/test_triangulation.py b/lib/matplotlib/tests/test_triangulation.py
@@ -985,6 +985,29 @@ def test_trirefiner_fortran_contiguous_triangles():
     assert_array_equal(fine_triang1.triangles, fine_triang2.triangles)
 
 
+def test_qhull_triangle_orientation():
+    # github issue 4437.
+    xi = np.linspace(-2, 2, 100)
+    x, y = map(np.ravel, np.meshgrid(xi, xi))
+    w = np.logical_and(x > y - 1, np.logical_and(x < -1.95, y > -1.2))
+    x, y = x[w], y[w]
+    theta = np.radians(25)
+    x1 = x*np.cos(theta) - y*np.sin(theta)
+    y1 = x*np.sin(theta) + y*np.cos(theta)
+
+    # Calculate Delaunay triangulation using Qhull.
+    triang = mtri.Triangulation(x1, y1)
+
+    # Neighbors returned by Qhull.
+    qhull_neighbors = triang.neighbors
+
+    # Obtain neighbors using own C++ calculation.
+    triang._neighbors = None
+    own_neighbors = triang.neighbors
+
+    assert_array_equal(qhull_neighbors, own_neighbors)
+
+
 if __name__ == '__main__':
     import nose
     nose.runmodule(argv=['-s', '--with-doctest'], exit=False)
diff --git a/lib/matplotlib/tri/_tri.cpp b/lib/matplotlib/tri/_tri.cpp
@@ -222,15 +222,17 @@ Triangulation::Triangulation(const CoordinateArray& x,
                              const TriangleArray& triangles,
                              const MaskArray& mask,
                              const EdgeArray& edges,
-                             const NeighborArray& neighbors)
+                             const NeighborArray& neighbors,
+                             int correct_triangle_orientations)
     : _x(x),
       _y(y),
       _triangles(triangles),
       _mask(mask),
       _edges(edges),
       _neighbors(neighbors)
 {
-    correct_triangles();
+    if (correct_triangle_orientations)
+        correct_triangles();
 }
 
 void Triangulation::calculate_boundaries()
diff --git a/lib/matplotlib/tri/_tri.h b/lib/matplotlib/tri/_tri.h
@@ -186,13 +186,17 @@ class Triangulation
      *          once.
      *   neighbors: Optional int array of shape (ntri,3) indicating which
      *              triangles are the neighbors of which TriEdges, or -1 if
-     *              there is no such neighbor. */
+     *              there is no such neighbor.
+     *   correct_triangle_orientations: Whether or not should correct triangle
+     *                                  orientations so that vertices are
+     *                                  ordered anticlockwise. */
     Triangulation(const CoordinateArray& x,
                   const CoordinateArray& y,
                   const TriangleArray& triangles,
                   const MaskArray& mask,
                   const EdgeArray& edges,
-                  const NeighborArray& neighbors);
+                  const NeighborArray& neighbors,
+                  int correct_triangle_orientations);
 
     /* Calculate plane equation coefficients for all unmasked triangles from
      * the point (x,y) coordinates and point z-array of shape (npoints) passed
diff --git a/lib/matplotlib/tri/_tri_wrapper.cpp b/lib/matplotlib/tri/_tri_wrapper.cpp
@@ -35,15 +35,17 @@ static int PyTriangulation_init(PyTriangulation* self, PyObject* args, PyObject*
     Triangulation::MaskArray mask;
     Triangulation::EdgeArray edges;
     Triangulation::NeighborArray neighbors;
+    int correct_triangle_orientations;
 
     if (!PyArg_ParseTuple(args,
-                          "O&O&O&O&O&O&",
+                          "O&O&O&O&O&O&i",
                           &x.converter, &x,
                           &y.converter, &y,
                           &triangles.converter, &triangles,
                           &mask.converter, &mask,
                           &edges.converter, &edges,
-                          &neighbors.converter, &neighbors)) {
+                          &neighbors.converter, &neighbors,
+                          &correct_triangle_orientations)) {
         return -1;
     }
 
@@ -80,7 +82,8 @@ static int PyTriangulation_init(PyTriangulation* self, PyObject* args, PyObject*
 
     CALL_CPP_INIT("Triangulation",
                   (self->ptr = new Triangulation(x, y, triangles, mask,
-                                                 edges, neighbors)));
+                                                 edges, neighbors,
+                                                 correct_triangle_orientations)));
     return 0;
 }
 
diff --git a/lib/matplotlib/tri/triangulation.py b/lib/matplotlib/tri/triangulation.py
@@ -107,7 +107,7 @@ def get_cpp_triangulation(self):
         if self._cpp_triangulation is None:
             self._cpp_triangulation = _tri.Triangulation(
                 self.x, self.y, self.triangles, self.mask, self._edges,
-                self._neighbors)
+                self._neighbors, not self.is_delaunay)
         return self._cpp_triangulation
 
     def get_masked_triangles(self):
