在不使用 Matplotlib 的情况下测试点是否在椭圆内？

Question

我正在从事基于 Python 的数据分析。我有一些 x-y 数据点和一些椭圆，我想确定点是否在任何椭圆内。我一直这样做的方式很有效，但它很笨拙。当我考虑将我的软件分发给其他人时，我发现自己想要一种更简洁的方式。

现在，我正在使用 matplotlib.patches.Ellipse 对象。 Matplotlib Ellipses 有一个有用的方法叫做contains_point()。您可以通过调用 Axes.transData.transform() 在 Matplotlib Axes 对象上处理数据坐标。

问题是我必须创建一个 Figure 和一个 Axes 对象来保存椭圆。当我的程序运行时，会渲染一个烦人的 Matplotlib Figure 对象，显示椭圆，我实际上不需要看到。我尝试了几种方法来抑制此输出。我已经使用Axes.clear() 成功地从轴中删除了椭圆，从而产生了一个空图。但我无法让 Matplolib 的 pyplot.close(fig_number) 在调用 pyplot.show() 之前删除图形本身。

感谢任何建议，谢谢！

Answer 1

受how a carpenter draws an ellipse 的启发，使用两个钉子和一根绳子，这是一个对 numpy 友好的实现，用于测试点是否位于给定的椭圆内。

椭圆的定义之一是到两个焦点的距离之和是恒定的，等于椭圆的宽度（或高度，如果它更大的话）。中心与焦点之间的距离为sqrt(a*a - b*b)，其中a 和b 是宽度和高度的一半。使用该距离和所需角度的旋转可以找到焦点的位置。 numpy.linalg.norm 可用于使用 numpy 的高效数组操作计算距离。

计算完成后，会生成一个图表以直观地检查一切是否正确。

import numpy as np
from numpy.linalg import norm # calculate the length of a vector

x = np.random.uniform(0, 40, 20000)
y = np.random.uniform(0, 20, 20000)
xy = np.dstack((x, y))
el_cent = np.array([20, 10])
el_width = 28
el_height = 17
el_angle = 20

# distance between the center and the foci
foc_dist = np.sqrt(np.abs(el_height * el_height - el_width * el_width) / 4)
# vector from center to one of the foci
foc_vect = np.array([foc_dist * np.cos(el_angle * np.pi / 180), foc_dist * np.sin(el_angle * np.pi / 180)])
# the two foci
el_foc1 = el_cent + foc_vect
el_foc2 = el_cent - foc_vect

# for each x,y: calculate z as the sum of the distances to the foci;
# np.ravel is needed to change the array of arrays (of 1 element) into a single array
z = np.ravel(norm(xy - el_foc1, axis=-1) + norm(xy - el_foc2, axis=-1) )
# points are exactly on the ellipse when the sum of distances is equal to the width
# z = np.where(z <= max(el_width, el_height), 1, 0)

# now create a plot to check whether everything makes sense
from matplotlib import pyplot as plt
from matplotlib import patches as mpatches

fig, ax = plt.subplots()
# show the foci as red dots
plt.plot(*el_foc1, 'ro')
plt.plot(*el_foc2, 'ro')
# create a filter to separate the points inside the ellipse
filter = z <= max(el_width, el_height)
# draw all the points inside the ellipse with the plasma colormap
ax.scatter(x[filter], y[filter], s=5, c=z[filter], cmap='plasma')
# draw all the points outside with the cool colormap
ax.scatter(x[~filter], y[~filter], s=5, c=z[~filter], cmap='cool')
# add the original ellipse to verify that the boundaries match
ellipse = mpatches.Ellipse(xy=el_cent, width=el_width, height=el_height, angle=el_angle,
                           facecolor='None', edgecolor='black', linewidth=2,
                           transform=ax.transData)
ax.add_patch(ellipse)
ax.set_aspect('equal', 'box')
ax.autoscale(enable=True, axis='both', tight=True)
plt.show()

Answer 2

这里最简单的解决方案是使用shapely。

如果你有一个形状为 Nx2 的数组，其中包含一组顶点 (xy)，那么构造适当的 shapely.geometry.polygon 对象并检查任意点或点集是否是微不足道的(points) 包含在 -

import shapely.geometry as geom
ellipse = geom.Polygon(xy)
for p in points:
    if ellipse.contains(geom.Point(p)):
        # ...

或者，如果椭圆由它们的参数（即旋转角、长半轴和短半轴）定义，则必须构造包含顶点的数组，然后应用相同的过程。我会推荐使用polar form relative to center，因为这与形状构造多边形的方式最兼容。

import shapely.geometry as geom
from shapely import affinity

n = 360
a = 2
b = 1
angle = 45

theta = np.linspace(0, np.pi*2, n)
r = a * b  / np.sqrt((b * np.cos(theta))**2 + (a * np.sin(theta))**2)
xy = np.stack([r * np.cos(theta), r * np.sin(theta)], 1)

ellipse = affinity.rotate(geom.Polygon(xy), angle, 'center')
for p in points:
    if ellipse.contains(geom.Point(p)):
        # ...

这种方法是有利的，因为它支持任何正确定义的多边形——不仅仅是椭圆，它不依赖于 matplotlib 方法来执行包含检查，并且它产生了一个非常易读的代码（这在“将 [一个人的] 软件分发给其他人”）。

这是一个完整的示例（添加了绘图以显示它的工作原理）

import shapely.geometry as geom
from shapely import affinity
import matplotlib.pyplot as plt
import numpy as np

n = 360
theta = np.linspace(0, np.pi*2, n)

a = 2
b = 1
angle = 45.0

r = a * b  / np.sqrt((b * np.cos(theta))**2 + (a * np.sin(theta))**2)
xy = np.stack([r * np.cos(theta), r * np.sin(theta)], 1)

ellipse = affinity.rotate(geom.Polygon(xy), angle, 'center')
x, y = ellipse.exterior.xy
# Create a Nx2 array of points at grid coordinates throughout
# the ellipse extent
rnd = np.array([[i,j] for i in np.linspace(min(x),max(x),50) 
                      for j in np.linspace(min(y),max(y),50)])
# Filter for points which are contained in the ellipse
res = np.array([p for p in rnd if ellipse.contains(geom.Point(p))])

plt.plot(x, y, lw = 1, color='k')
plt.scatter(rnd[:,0], rnd[:,1], s = 50, color=(0.68, 0.78, 0.91)
plt.scatter(res[:,0], res[:,1], s = 15, color=(0.12, 0.67, 0.71))
plt.show()