Проверьте, находятся ли точки внутри эллипсов, без использования Matplotlib?

Question

Я работаю над анализом данных на основе Python. У меня есть несколько точек данных xy и несколько эллипсов, и я хочу определить, находятся ли точки внутри какого-либо из эллипсов. То, как я это делаю, работает, но это глупо. Размышляя о распространении своего программного обеспечения среди других людей, я обнаружил, что хочу найти более чистый способ.

Сейчас я использую matplotlib.patches.Ellipse объекты. Matplotlib Ellipses имеет полезный метод, который называется contains_point(). Вы можете работать с координатами данных на объекте Matplotlib Axes, вызывая Axes.transData.transform().

. Суть в том, что мне нужно создать объект Figure и Axes для хранения эллипсов. И когда моя программа запустится, будет выдан раздражающий объект Matplotlib Figure, показывающий эллипсы, которые мне на самом деле не нужны. Я пробовал несколько методов, чтобы подавить этот вывод. Мне удалось удалить эллипсы с осей, используя Axes.clear(), что привело к пустому графику. Но я не могу заставить pyplot.close(fig_number) Матплолиба удалить саму фигуру, прежде чем звонить pyplot.show().

Любой совет приветствуется, спасибо!

Answer 1

Самое простое решение здесь - использовать 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)):
        # ...

В качестве альтернативы, если эллипсы определяются их параметрами (т. е. углом поворота, большой полуосью и полуминорной осью), то массив, содержащий вершины, должен быть построен и затем применяется тот же процесс. Я бы порекомендовал использовать полярную форму относительно центра , так как это наиболее совместимо с тем, как стройно строятся полигоны.

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

Answer 2

Вдохновленный как плотник др aws эллипс , использующий два гвоздя и кусок струны, вот реализация, удобная для numpy, чтобы проверить, находятся ли точки l ie внутри данных эллипсов.

Одним из определений эллипса является то, что сумма расстояний до двух фокусов постоянна и равна ширине (или высоте, если она будет больше) эллипса. Расстояние между центром и фокусами составляет 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()