Improved qhull triangulations with large x,y offset

pllim · Jul 14, 2017 · be920ad · be920ad
1 parent 8e517d8
commit be920ad
Show file tree

Hide file tree

Showing 5 changed files with 29 additions and 2 deletions.
diff --git a/doc/api/api_changes/2017-07-14-IT_qhull_large_offset.rst b/doc/api/api_changes/2017-07-14-IT_qhull_large_offset.rst
@@ -0,0 +1,6 @@
+Improved Delaunay triangulations with large offsets
+```````````````````````````````````````````````````
+
+Delaunay triangulations now deal with large x,y offsets in a better
+way. This can cause minor changes to any triangulations calculated
+using Matplotlib.
diff --git a/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_contouring.png b/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_contouring.png
diff --git a/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_gradient.png b/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_gradient.png
diff --git a/lib/matplotlib/tests/test_triangulation.py b/lib/matplotlib/tests/test_triangulation.py
@@ -1123,3 +1123,14 @@ def test_internal_cpp_api():
     with pytest.raises(ValueError) as excinfo:
         trifinder.find_many([0], [0, 1])
     excinfo.match(r'x and y must be array_like with same shape')
+
+
+def test_qhull_large_offset():
+    # github issue 8682.
+    x = np.asarray([0, 1, 0, 1, 0.5])
+    y = np.asarray([0, 0, 1, 1, 0.5])
+
+    offset = 1e10
+    triang = mtri.Triangulation(x, y)
+    triang_offset = mtri.Triangulation(x + offset, y + offset)
+    assert len(triang.triangles) == len(triang_offset.triangles)
diff --git a/src/qhull_wrap.c b/src/qhull_wrap.c
@@ -108,6 +108,8 @@ delaunay_impl(int npoints, const double* x, const double* y,
     PyArrayObject* neighbors = NULL;
     int* triangles_ptr;
     int* neighbors_ptr;
+    double x_mean = 0.0;
+    double y_mean = 0.0;
 
     QHULL_LIB_CHECK
 
@@ -119,10 +121,18 @@ delaunay_impl(int npoints, const double* x, const double* y,
         goto error_before_qhull;
     }
 
+    /* Determine mean x, y coordinates. */
+    for (i = 0; i < npoints; ++i) {
+        x_mean += x[i];
+        y_mean += y[i];
+    }
+    x_mean /= npoints;
+    y_mean /= npoints;
+
     /* Prepare points array to pass to qhull. */
     for (i = 0; i < npoints; ++i) {
-        points[2*i  ] = x[i];
-        points[2*i+1] = y[i];
+        points[2*i  ] = x[i] - x_mean;
+        points[2*i+1] = y[i] - y_mean;
     }
 
     /* qhull expects a FILE* to write errors to. */