python - 检查点是否在椭圆内比 contains_point 方法更快

Question

我使用 matplotlib 1.15.1 并尝试生成这样的散点图:

椭圆的大小是固定的，用中心坐标、宽度、高度和角度(从外部提供)绘制:我不知道它们的方程式是什么。

g_ell_center = (0.8882, 0.8882)
g_ell_width = 0.36401857095483
g_ell_height = 0.16928136341606
g_ellipse = patches.Ellipse(g_ell_center, g_ell_width, g_ell_height, angle=angle, fill=False, edgecolor='green', linewidth=2)

这个省略号应该在我的图上标记正常和半正常数据。然后，我有一个 ~500 点的数组，必须根据它们所属的椭圆进行着色。所以我尝试用 contains_point 方法检查每个点:

colors_array = []
colors_scheme = ['green', 'yellow', 'black']
for point in points_array:
    if g_ellipse.contains_point(point, radius=0):
        colors_array.append(0)
    elif y_ellipse.contains_point(point, radius=0):
        colors_array.append(1)
    else:
        colors_array.append(2)

最后绘制点:

plt.scatter(x_array, y_array, s=10, c=[colors_scheme[x] for x in colors_array], edgecolor="k", linewidths=0.3)

但是 contains_point 非常慢!它为 300 点散点图工作了 5 分钟，我必须并行生成数千个散点图。也许有更快的方法？附言整个项目绑定(bind)到 matplotlib，我不能使用其他库。

Answer 1

此方法应在给定椭圆的中心、宽度、高度和角度的情况下测试点是否在椭圆内。您找到该点相对于椭圆中心的 x 和 y 坐标，然后使用角度将其转换为沿长轴和短轴的坐标。最后，您会找到该点与像元中心的归一化距离，其中距离为 1 表示在椭圆上，小于 1 表示在内部，大于 1 表示在外部。

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np

fig,ax = plt.subplots(1)
ax.set_aspect('equal')

# Some test points
x = np.random.rand(500)*0.5+0.7
y = np.random.rand(500)*0.5+0.7

# The ellipse
g_ell_center = (0.8882, 0.8882)
g_ell_width = 0.36401857095483
g_ell_height = 0.16928136341606
angle = 30.

g_ellipse = patches.Ellipse(g_ell_center, g_ell_width, g_ell_height, angle=angle, fill=False, edgecolor='green', linewidth=2)
ax.add_patch(g_ellipse)

cos_angle = np.cos(np.radians(180.-angle))
sin_angle = np.sin(np.radians(180.-angle))

xc = x - g_ell_center[0]
yc = y - g_ell_center[1]

xct = xc * cos_angle - yc * sin_angle
yct = xc * sin_angle + yc * cos_angle 

rad_cc = (xct**2/(g_ell_width/2.)**2) + (yct**2/(g_ell_height/2.)**2)

# Set the colors. Black if outside the ellipse, green if inside
colors_array = np.array(['black'] * len(rad_cc))
colors_array[np.where(rad_cc <= 1.)[0]] = 'green'

ax.scatter(x,y,c=colors_array,linewidths=0.3)

plt.show()

请注意，整个脚本需要 0.6 秒才能运行并处理 500 个点。这包括创建和保存图形等。

使用上述np.where 方法设置colors_array 的过程需要0.00007s 500 点。

请注意，在下面显示的旧实现中，在循环中设置 colors_array 需要 0.00016 秒:

colors_array = []

for r in rad_cc:
    if r <= 1.:
        # point in ellipse
        colors_array.append('green')
    else:
        # point not in ellipse
        colors_array.append('black')