ENH: Voronoi diagram gives unexpected results from scipy.spatial about scipy HOT 6 CLOSED

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?

from scipy.

I can reproduce this bug with the following simplified code (scipy v1.11.4):

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:

I confirm this with 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]

from scipy.

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 --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:

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):

--- 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])

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.

--- 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')

from scipy.

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.

from scipy.

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.

from scipy.

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.

from scipy.