不使用Matplotlib如何测试点是否在椭圆内？

Question

我正在进行基于Python的数据分析工作。我有一些x-y数据点和一些椭圆形，希望确定点是否在任何一个椭圆形内。我已经实现了一种方法，但它比较笨拙。考虑到将我的软件分发给其他人使用，我想要更简洁的方法。

目前，我正在使用matplotlib.patches.Ellipse对象。Matplotlib中的椭圆形具有一个有用的方法contains_point()。你可以通过调用Axes.transData.transform()在Matplotlib Axes对象上以数据坐标的方式工作。

问题在于，我必须创建一个包含椭圆形的Figure和Axes对象。当程序运行时，Matplotlib Figure对象会得到渲染，并显示椭圆形，而我实际上不需要看到它们。我尝试了几种方法来抑制这个输出。我成功地从Axes中删除了椭圆形，使用Axes.clear()，导致出现一个空图表。但在调用pyplot.show()之前，我无法使用Matplolib的pyplot.close(fig_number)删除Figure本身。

感谢您的帮助！

- John Ladasky

---

省略号是象征性定义的还是由点集合组成的？ - William Miller

我可以提供椭圆的尺寸参数。我刚刚花了一个小时研究Matplotlib源代码。虽然Matplotlib不将椭圆视为光栅点集合，但它似乎将所有闭合多边形转换为“路径”，即周长上的点集合。所有闭合多边形都有一个共享的contains_point()方法。它比我意识到的更通用，也许速度更慢。我很高兴改变策略，不使用Matplotlib。 - John Ladasky

我认为Shapely库比Matplotlib更适合此任务，你可能需要研究一下。 - William Miller

请参考这里的解析解 https://math.stackexchange.com/questions/76457/check-if-a-point-is-within-an-ellipse。 - Diziet Asahi

请看我的答案，了解Shapely方法。 - William Miller

你有没有看过这个问题的发布答案的机会？ - William Miller

Answer 1

受木匠如何用两个钉子和一根绳子画椭圆的启发，这里提供了一个友好的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()

- JohanC

Answer 2

最简单的解决方案是使用shapely。如果您有一个包含一组顶点（xy）的Nx2形状数组，则轻松构建适当的shapely.geometry.polygon对象并检查是否包含任意点或一组点（points）。

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

或者，如果省略号是通过其参数定义的（即旋转角度、半长轴和半短轴），则必须构建包含顶点的数组，然后应用相同的过程。我建议使用相对于中心的极坐标形式，因为这与shapely构造多边形的方式最兼容。

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

- William Miller