Repository: scipy/scipy
Pull request: ENH, MAINT: voronoi_plot_2d nicer inf lines
* Fixes #19653

* we currently set the plotted length of infinitely-long Voronoi edges as a function of the largest x-or-y Voronoi generator range; in cases where the range of one axis is substantially larger than the other, it seems that we can produce a more appealing representation of the infinitely-long line by accounting for the aspect ratio; otherwise, there's a danger the (infinitely-long) Voronoi edge may not be visible with our default `voronoi_plot_2d` settings

* note that we only have very crude testing for `voronoi_plot_2d`, so there is no guarantee this doesn't introduce new problems for other scenarios (so not keen on backporting, even if folks thought this was a genuine bug), and I'm also not super keen on extracting line data from matplotlib like this, but I do also think this probably is a step in the right testing direction if we're going to respond to user feedback in a sane way; the function already takes 8 parameters, and to some extent the user should be expected to make reasonable adjustments, but since the Voronoi edge was invisible in this default case I thought we might try to fix it

* while in the neighborhood, remove some local test shims for `matplotlib` `1.x`; `matplotlib` is already at `3.x` so I think we can drop those

[skip cirrus]

cc @bagrow and @dpapadop for feedback, since they were vocal on the original ticket

Before/after for the generators in question (it may be good to spot check a few more cases and perhaps even suggest another regression case or two):

![image](https://github.com/scipy/scipy/assets/7903078/57f199dc-6ef9-4974-aa0a-279f837da411)

![image](https://github.com/scipy/scipy/assets/7903078/016c53d4-6b5a-4267-b2a8-5d99cc8927f5)


Problem statement:
ENH: Voronoi diagram gives unexpected results from scipy.spatial
### Describe your issue.

When I try to plot the Voronoi diagram using the scipy.spatial, the results show (please see in the reproducing code example section), insstead of getting 5 areas, I am getting 4, even though I provide 5 points in the 2D space. Moreover the space separation seems to be a bit off and they are all rectangles

### Reproducing Code Example

```python
import pandas as pd
pd.DataFrame({'cl': {0: 'A', 1: 'C', 2: 'H', 3: 'M', 4: 'S'},
 'd': {0: 245.059986986012,
  1: 320.49044143557785,
  2: 239.79023081978914,
  3: 263.38325791238833,
  4: 219.53334398353175},
 'p': {0: 10.971011721360075,
  1: 10.970258360366753,
  2: 13.108487516946218,
  3: 12.93241352743668,
  4: 13.346107628161008}})

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d

centers2 = np.array(
    centers_dt[['d', 'p']]
)

scatter_x = np.array(centers_dt['d'])
scatter_y = np.array(centers_dt['p'])
group = np.array(centers_dt['cl'])
cdict = {'C': 'red', 'A': 'blue', 'H': 'green', 'M': 'yellow', 'S': 'black'}

fig, ax = plt.subplots()
for g in np.unique(group):
    ix = np.where(group == g)
    ax.scatter(scatter_x[ix], scatter_y[ix], c = cdict[g], label = g, s = 100)
ax.legend()

vor = Voronoi(centers2)
fig = voronoi_plot_2d(vor,plt.gca())

plt.show()
plt.close()
```


### Error message

```shell
There is no error message, it is just that the results are unexpected.
```


### SciPy/NumPy/Python version and system information

```shell
1.7.0 1.21.6 sys.version_info(major=3, minor=7, micro=10, releaselevel='final', serial=0)
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
```

Hints from issue discussion:
Two potential problems I'm seeing here:

1. When I try to run the reproducing code I get:

```
Traceback (most recent call last):
  File "/Users/treddy/rough_work/test.py", line 19, in <module>
    centers_dt[['d', 'p']]
    ^^^^^^^^^^
NameError: name 'centers_dt' is not defined
```

2. You're using SciPy `1.7.0`, which is over two years old. I suppose it isn't ancient, but could you try a more recent version as well?
I can reproduce this bug with the following simplified code (scipy v1.11.4):

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d

# Renaming the array to 'points'
points = np.array([
    [245.059986986012, 10.971011721360075],
    [320.49044143557785, 10.970258360366753],
    [239.79023081978914, 13.108487516946218],
    [263.38325791238833, 12.93241352743668],
    [219.53334398353175, 13.346107628161008]
])

# Create a Voronoi tessellation for these points
vor = Voronoi(points)

# Plot the Voronoi tessellation
fig, ax = plt.subplots()
voronoi_plot_2d(vor, ax)

# Overlay the original points
ax.plot(points[:, 0], points[:, 1], 'ko')

ax.set_title('Voronoi Tessellation of Points')
plt.show()
```

The resulting plot appears to have four cells instead of five:
![voronoi-test](https://github.com/scipy/scipy/assets/395952/d13cbf79-d158-4bcd-a70a-5c1a70938bed)

I confirm this with Mathematica:
```mathematica
(*Define the points*)
points = {{245.059986986012, 10.971011721360075}, {320.49044143557785,
     10.970258360366753}, {239.79023081978914, 
    13.108487516946218}, {263.38325791238833, 
    12.93241352743668}, {219.53334398353175, 13.346107628161008}};

(*Create the Voronoi tessellation*)
voronoiMesh = VoronoiMesh[points];

(*Plot the Voronoi diagram with the points overlaid and adjusted \
aspect ratio*)
Show[voronoiMesh, Graphics[{Red, PointSize[Medium], Point[points]}], 
 AspectRatio -> 1/GoldenRatio]
```
![voronoi-test-mathematica](https://github.com/scipy/scipy/assets/395952/fb07e4b2-c5e0-4a27-9837-1dbf55080ee0)




I think I'm still a bit skeptical of the bug claim at this point. For example, if I adjust the limits of the plot a bit, and include some of the finite Voronoi vertices with:

```diff
diff --git a/test.py b/test.py
index b7d5ace..c11f216 100644
--- a/test.py
+++ b/test.py
@@ -16,10 +16,14 @@ vor = Voronoi(points)
 
 # Plot the Voronoi tessellation
 fig, ax = plt.subplots()
-voronoi_plot_2d(vor, ax)
+voronoi_plot_2d(vor,
+                ax,
+                show_vertices=True)
 
 # Overlay the original points
 ax.plot(points[:, 0], points[:, 1], 'ko')
+ax.set_xlim(200, 350)
+ax.set_ylim(-600, 6000)
```

I can start to see small portions of the "dashed lines" that outline the Voronoi edges that head to infinity, which does indeed appear to encompass the five Voronoi regions, if you use your imagination a bit to extend the dashed lines outward:

![image](https://github.com/scipy/scipy/assets/7903078/c3e61760-b9cc-4ca4-be1a-f6fc6ee213af)

I suppose it might be helpful to have a way to auto-plot this more intuitively, but `vor.regions` also clearly has five regions, most of which include a Voronoi vertex at infinity:

`[[], [0, -1], [3, 1, 2], [3, 0, -1, 1], [2, -1, 1], [3, 0, -1, 2]]` (the empty list is just an internal implementation detail for a generator at infinity used by Qhull)

To really drive the point home, we can try extending the dashed lines a bit by modifying the SciPy `voronoi_plot_2d` source (maybe we could come up with a nice way to allow a parameter for this, but as far as this being a genuine bug in the computational geometry algorithm vs. a desired plotting enhancement, I remain skeptical):

```diff
--- a/scipy/spatial/_plotutils.py
+++ b/scipy/spatial/_plotutils.py
@@ -249,7 +249,7 @@ def voronoi_plot_2d(vor, ax=None, **kw):
             direction = np.sign(np.dot(midpoint - center, n)) * n
             if (vor.furthest_site):
                 direction = -direction
-            far_point = vor.vertices[i] + direction * ptp_bound.max()
+            far_point = vor.vertices[i] + direction * ptp_bound.max() * 100
 
             infinite_segments.append([vor.vertices[i], far_point])
```

![image](https://github.com/scipy/scipy/assets/7903078/81280d19-2fc1-46a7-ba7a-c3fd4a16b3cb)

Now it looks to me like our Voronoi diagram is actually very realistic, and this is mostly just a matter of plotting limits/aspect ratio type stuff. If I remove my modification to the SciPy source, and instead adjust the generators such that the their `x` and `y` coords have more similar ranges/magnitudes, we again get something pretty sensible looking (see below), so I'm feeling fairly confident that there's no evidence of a genuine bug and maybe just a possibility for a small plotting enhancement on the edges to infinity. Conversely, we probably can't really expect `voronoi_plot_2d` to handle every situation perfectly and automatically--there will always be some pathological data where the algorithm itself is correct, and we just have to ask the user to dig a bit deeper in terms of making the plot look nice.

```diff
--- a/test.py
+++ b/test.py
@@ -10,13 +10,16 @@ points = np.array([
     [263.38325791238833, 12.93241352743668],
     [219.53334398353175, 13.346107628161008]
 ])
+points[..., 0] /= 10
 
 # Create a Voronoi tessellation for these points
 vor = Voronoi(points)
 
 # Plot the Voronoi tessellation
 fig, ax = plt.subplots()
-voronoi_plot_2d(vor, ax)
+voronoi_plot_2d(vor,
+                ax,
+                show_vertices=True)
 
 # Overlay the original points
 ax.plot(points[:, 0], points[:, 1], 'ko')
```

![image](https://github.com/scipy/scipy/assets/7903078/ec32514a-30c8-4f82-8f73-e73001f8f8ea)



It's clear that this is a problem with the default visualization and not computing the tesselation. And the example data are definitely pathological; with that aspect ratio, it _literally_ is an edge case. 

With that said, in an ideal world it would be nice to handle this correctly, automatically. Yes, there is an opportunity cost and you don't want to make the plotting code worse in some way to handle this. Again, Mathematica handled it without problem. I don't see why scipy can't.
Well-tested PRs are always welcome. In the meanwhile, I'll switch the label/title to "enhancement," since there's no algorithmic bug I'd likely want to backport a fix for, just would be nice to have the plots more appealing in more cases.

I'll do what I can to patch this, the `ptp` approach to setting bounds on the generators alone probably isn't the best way to do it, but since it is plotting code defending against regressions is inevitably more annoying.
Happy to contribute something substantive, but not sure I'd trust my code. I'll defer to your taxonomy, but I still think it's a 'defect.' The plot is not unappealing, it is wrong.

Commit message: ENH, MAINT: voronoi_plot_2d nicer inf lines

* Fixes #19653

* we currently set the plotted length of infinitely-long
Voronoi edges as a function of the largest x-or-y Voronoi
generator range; in cases where the range of one axis is
substantially larger than the other, it seems that we
can produce a more appealing representation of the infinitely-long
line by accounting for the aspect ratio; otherwise, there's a danger
the Voronoi edge may not be visible with our default `voronoi_plot_2d`
settings

* note that we only have very crude testing for `voronoi_plot_2d`,
so there is no guarantee this doesn't introduce new problems
for other scenarios, and I'm also not super keen on extracting
line data from matplotlib like this, but I do also think this
probably is a step in the right testing direction if we're going to respond
to user feedback in a sane way; the function already takes 8 parameters,
and to some extent the user should be expected to make reasonable
adjustments, but since the Voronoi edge was invisible in this
default case I thought we might try to fix it

* while in the neighborhood, remove some local test
shims for `matplotlib` `1.x`; `matplotlib` is already
at `3.x` so I think we can drop those

[skip cirrus]

Patch to review:
diff --git a/scipy/spatial/_plotutils.py b/scipy/spatial/_plotutils.py
index 7aa99889330d..9e7e0998a5a7 100644
--- a/scipy/spatial/_plotutils.py
+++ b/scipy/spatial/_plotutils.py
@@ -249,7 +249,8 @@ def voronoi_plot_2d(vor, ax=None, **kw):
             direction = np.sign(np.dot(midpoint - center, n)) * n
             if (vor.furthest_site):
                 direction = -direction
-            far_point = vor.vertices[i] + direction * ptp_bound.max()
+            aspect_factor = abs(ptp_bound.max() / ptp_bound.min())
+            far_point = vor.vertices[i] + direction * ptp_bound.max() * aspect_factor
 
             infinite_segments.append([vor.vertices[i], far_point])
 
diff --git a/scipy/spatial/tests/test__plotutils.py b/scipy/spatial/tests/test__plotutils.py
index c85d218763bf..54e11a5a2d87 100644
--- a/scipy/spatial/tests/test__plotutils.py
+++ b/scipy/spatial/tests/test__plotutils.py
@@ -1,5 +1,7 @@
 import pytest
-from numpy.testing import assert_, assert_array_equal, suppress_warnings
+import numpy as np
+from numpy.testing import assert_, assert_array_equal, assert_allclose
+
 try:
     import matplotlib
     matplotlib.rcParams['backend'] = 'Agg'
@@ -22,10 +24,7 @@ def test_delaunay(self):
         fig = plt.figure()
         obj = Delaunay(self.points)
         s_before = obj.simplices.copy()
-        with suppress_warnings() as sup:
-            # filter can be removed when matplotlib 1.x is dropped
-            sup.filter(message="The ishold function was deprecated in version")
-            r = delaunay_plot_2d(obj, ax=fig.gca())
+        r = delaunay_plot_2d(obj, ax=fig.gca())
         assert_array_equal(obj.simplices, s_before)  # shouldn't modify
         assert_(r is fig)
         delaunay_plot_2d(obj, ax=fig.gca())
@@ -34,10 +33,7 @@ def test_voronoi(self):
         # Smoke test
         fig = plt.figure()
         obj = Voronoi(self.points)
-        with suppress_warnings() as sup:
-            # filter can be removed when matplotlib 1.x is dropped
-            sup.filter(message="The ishold function was deprecated in version")
-            r = voronoi_plot_2d(obj, ax=fig.gca())
+        r = voronoi_plot_2d(obj, ax=fig.gca())
         assert_(r is fig)
         voronoi_plot_2d(obj)
         voronoi_plot_2d(obj, show_vertices=False)
@@ -46,9 +42,29 @@ def test_convex_hull(self):
         # Smoke test
         fig = plt.figure()
         tri = ConvexHull(self.points)
-        with suppress_warnings() as sup:
-            # filter can be removed when matplotlib 1.x is dropped
-            sup.filter(message="The ishold function was deprecated in version")
-            r = convex_hull_plot_2d(tri, ax=fig.gca())
+        r = convex_hull_plot_2d(tri, ax=fig.gca())
         assert_(r is fig)
         convex_hull_plot_2d(tri)
+
+    def test_gh_19653(self):
+        # aspect ratio sensitivity of voronoi_plot_2d
+        # infinite Voronoi edges
+        points = np.array([[245.059986986012, 10.971011721360075],
+                           [320.49044143557785, 10.970258360366753],
+                           [239.79023081978914, 13.108487516946218],
+                           [263.38325791238833, 12.93241352743668],
+                           [219.53334398353175, 13.346107628161008]
+                          ])
+        vor = Voronoi(points)
+        fig = voronoi_plot_2d(vor)
+        ax = fig.gca()
+        infinite_segments = ax.collections[1].get_segments()
+        expected_segments = np.array([[[282.77256, -254.76904],
+                                      [282.729714, -4544.744698]],
+                                     [[282.77256014, -254.76904029],
+                                      [430.08561382, 4032.67658742]],
+                                     [[229.26733285,  -20.39957514],
+                                      [-168.17167404, -4291.92545966]],
+                                     [[289.93433364, 5151.40412217],
+                                      [330.40553385, 9441.18887532]]])
+        assert_allclose(infinite_segments, expected_segments)


Code review comments:
diff_hunk: @@ -46,9 +42,29 @@ def test_convex_hull(self):
         # Smoke test
         fig = plt.figure()
         tri = ConvexHull(self.points)
-        with suppress_warnings() as sup:
-            # filter can be removed when matplotlib 1.x is dropped
-            sup.filter(message="The ishold function was deprecated in version")
-            r = convex_hull_plot_2d(tri, ax=fig.gca())
+        r = convex_hull_plot_2d(tri, ax=fig.gca())
         assert_(r is fig)
         convex_hull_plot_2d(tri)
+
+    def test_gh_19653(self):
+        # aspect ratio sensitivity of voronoi_plot_2d
+        # infinite Voronoi edges
+        points = np.array([[245.059986986012, 10.971011721360075],
+                           [320.49044143557785, 10.970258360366753],
+                           [239.79023081978914, 13.108487516946218],
+                           [263.38325791238833, 12.93241352743668],
+                           [219.53334398353175, 13.346107628161008]
+                          ])
+        vor = Voronoi(points)
+        fig = voronoi_plot_2d(vor)
+        ax = fig.gca()
+        infinite_segments = ax.collections[1].get_segments()
+        expected_segments = np.array([[[282.77256, -254.76904],
+                                      [282.729714, -4544.744698]],
+                                     [[282.77256014, -254.76904029],
+                                      [430.08561382, 4032.67658742]],
+                                     [[229.26733285,  -20.39957514],
+                                      [-168.17167404, -4291.92545966]],
+                                     [[289.93433364, 5151.40412217],
+                                      [330.40553385, 9441.18887532]]])
path: scipy/spatial/tests/test__plotutils.py
text: @user1:
```suggestion
        expected_segments = np.array([[[282.77256, -254.76904],
                                       [282.729714, -4544.744698]],
                                      [[282.77256014, -254.76904029],
                                       [430.08561382, 4032.67658742]],
                                      [[229.26733285,  -20.39957514],
                                       [-168.17167404, -4291.92545966]],
                                      [[289.93433364, 5151.40412217],
                                       [330.40553385, 9441.18887532]]])
```

diff_hunk: @@ -46,9 +42,29 @@ def test_convex_hull(self):
         # Smoke test
         fig = plt.figure()
         tri = ConvexHull(self.points)
-        with suppress_warnings() as sup:
-            # filter can be removed when matplotlib 1.x is dropped
-            sup.filter(message="The ishold function was deprecated in version")
-            r = convex_hull_plot_2d(tri, ax=fig.gca())
+        r = convex_hull_plot_2d(tri, ax=fig.gca())
         assert_(r is fig)
         convex_hull_plot_2d(tri)
+
+    def test_gh_19653(self):
+        # aspect ratio sensitivity of voronoi_plot_2d
+        # infinite Voronoi edges
+        points = np.array([[245.059986986012, 10.971011721360075],
+                           [320.49044143557785, 10.970258360366753],
+                           [239.79023081978914, 13.108487516946218],
+                           [263.38325791238833, 12.93241352743668],
+                           [219.53334398353175, 13.346107628161008]
+                          ])
path: scipy/spatial/tests/test__plotutils.py
text: @user1:
```suggestion
        points = np.array([[245.059986986012, 10.971011721360075],
                           [320.49044143557785, 10.970258360366753],
                           [239.79023081978914, 13.108487516946218],
                           [263.38325791238833, 12.93241352743668],
                           [219.53334398353175, 13.346107628161008]])
```