椭圆内点的计数

Question

我正在尝试计算每个椭圆环内给定数据点的数量:

问题是我有一个用于检查以下内容的函数: 因此为了确保一个点是否在每个椭圆中，需要计算三个输入参数:

def get_focal_point(r1,r2,center_x):
    # f = square root of r1-squared - r2-squared
    focal_dist = sqrt((r1**2) - (r2**2))
    f1_x = center_x - focal_dist
    f2_x = center_x + focal_dist
    return f1_x, f2_x

def get_distance(f1,f2,center_y,t_x,t_y):
    d1 = sqrt(((f1-t_x)**2) + ((center_y - t_y)**2)) 
    d2 = sqrt(((f2-t_x)**2) + ((center_y - t_y)**2))
    return d1,d2

def in_ellipse(major_ax,d1,d2):
    if (d1+d2) <= 2*major_ax:
        return True
    else:
        return False

现在我正在检查它是否在椭圆内：

for i in range(len(data.latitude)):
    t_x = data.latitude[i] 
    t_y = data.longitude[i] 
    d1,d2 = get_distance(f1,f2,center_y,t_x,t_y)
    d1_array.append(d1)
    d2_array.append(d2)
    if in_ellipse(major_ax,d1,d2) == True:
        core_count += 1
        # if the point is not in core ellipse 
        # check the next ring up
    else:
        for i in range(loop):
            .....

但我需要计算外部循环的每一对焦点...有没有更有效或更聪明的方法来完成这个任务？

- Florie

---

我是否正确，这些椭圆具有相同的比例和中心？如果是这样，您可以通过例如 r1 来识别它们，并创建一个函数，为每个点分配最小的 r1，以使该点位于椭圆内。我是正确的还是我误解了什么？ - Tadeck

是的，它们是同心的！但是dist_to_point/size_per_ellipse这个公式只适用于圆形吗？ - Florie

@g.d.d.c：这与我所说的类似，只是“到点的距离”并不能确定它何时落下，更重要的是“纬度”和“经度”这一对，因为这些椭圆不是圆形。 - Tadeck

没错！在我看来，@g.d.d.c的解决方案只适用于圆形。 - Tadeck

@endolith：确实可以，但是考虑到距离某个点（中心？）的距离与椭圆大小之间的关系，你无法确定某个点是否落在某个椭圆内。因为即使在同一条椭圆上，距离也会发生变化，所以基本上你需要同时知道x和y才能确定某个点是在椭圆内还是外部。 - Tadeck

Answer 1

这可能与您正在做的事情类似。 我只是想看看如果 f(x,y) = x^2/r1^2 + y^2/r2^2 = 1，那么会发生什么。

当 f(x,y) 大于1时，点 x,y 在椭圆外部。 当它小于1时，则在椭圆内部。我循环遍历每个椭圆来找到 f(x,y) 小于1 的那个。

该代码还没有考虑到偏离原点的椭圆。 包括此功能是一个小改变。

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

def inWhichEllipse(x,y,rads):
    '''
    With a list of (r1,r2) pairs, rads, return the index of the pair in which
    the point x,y resides. Return None as the index if it is outside all 
    Ellipses.
    '''
    xx = x*x
    yy = y*y

    count = 0
    ithEllipse =0
    while True:
        rx,ry = rads[count]
        ellips = xx/(rx*rx)+yy/(ry*ry)
        if ellips < 1:
            ithEllipse = count
            break
        count+=1
        if count >= len(rads):
            ithEllipse = None
            break

    return ithEllipse

rads = zip(np.arange(.5,10,.5),np.arange(.125,2.5,.25))

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlim(-15,15)
ax.set_ylim(-15,15)

# plot Ellipses
for rx,ry in rads:
    ellipse = patches.Ellipse((0,0),rx*2,ry*2,fc='none',ec='red')    
    ax.add_patch(ellipse)

x=3.0
y=1.0
idx = inWhichEllipse(x,y,rads)
rx,ry = rads[idx]
ellipse = patches.Ellipse((0,0),rx*2,ry*2,fc='none',ec='blue')    
ax.add_patch(ellipse)

if idx != None:
    circle = patches.Circle((x,y),.1)
    ax.add_patch(circle)

plt.show()

这段代码生成以下图形：

请注意，这只是一个起点。你可以改变inWhichEllipse来接受r1和r2的平方的列表，即（r1*r1，r2*r2）对，这会进一步减少计算量。

- Yann

Answer 2

你把事情复杂化了。根据椭圆的几何定义，计算焦点和到焦点的距离等是没有必要的。如果你知道长轴和短轴（你知道的），只需稍微压缩整个问题（例如通过将x-centerx和y-centery除以xaxis和yaxis使两者均为1.0），那么问题是否在椭圆内部就变得非常简单。

xnormalized**2 + ynormalized**2 <= 1

P.S.: 一般来说，在这个领域里的一个好建议是：如果你可以不实际计算距离而保持在其平方的范围内，就不要使用sqrt。

- user1046334

---

+1 我认为 @herby 表达得比我好多了。这个问题可以大大简化。 - Raymond Hettinger

Answer 3

以下是一些建议：

• 将计算焦点的代码移出循环是正确的做法。
• 通过去掉平方根可以加速距离计算。换句话说，我们知道 a < b 意味着 sqrt(a) < sqrt(b)，因此没有必要计算平方根。
• 如果椭圆是同心的，并且长轴平行于 x 轴，则可以通过重新调整 x 值将椭圆问题简化为圆问题。

此外，这里有一个小的编码问题。没有必要使用 if 语句 来返回 True 或 False。相反，可以直接返回条件表达式本身：

def in_ellipse(major_ax,d1,d2):
    return (d1+d2) <= 2*major_ax:

- Raymond Hettinger

---

非常感谢，这非常有帮助。虽然主轴不总是x轴...我将在循环中处理许多数据集，因此主轴不一定是x轴。 - Florie

如果轴通过原点，您仍然可以使用单个复数乘法将椭圆问题旋转和缩放为圆问题。此外，如果将x、y表示为复数，则距离计算简化为“abs（p）”。 - Raymond Hettinger