计算椭圆内的点数

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):
            .....

但我还需要计算外圈的每对焦点……有没有更高效或者更聪明的方法来做到这一点？

Answer 1

这里有几个建议给你：

• 你做得对，把计算焦点的代码放到循环外面。
• 计算距离的时候可以加快速度，去掉平方根的计算。换句话说，我们知道如果 a < b，那么 sqrt(a) < sqrt(b)，所以其实不需要算平方根。
• 如果这些椭圆是同心的，并且长轴和x轴平行，你可以通过调整x值，把椭圆的问题简化成圆的问题。

另外，还有一个小的编码建议。其实不需要用if语句来返回True或False。你可以直接返回条件表达式本身：

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

Answer 2

你把事情搞复杂了。其实没必要根据椭圆的几何定义去计算焦点和焦点之间的距离。如果你知道长轴和短轴的长度（你是知道的），那么只需要稍微调整一下问题（比如把长轴和短轴都变成1.0，可以通过把x中心和y中心分别除以长轴和短轴来实现），那么判断一个点是否在椭圆内部的问题就变得简单了。

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

补充一下：在这个领域，有个好建议就是：如果可以不计算实际的距离，而是直接用距离的平方来处理，就不要用sqrt函数。

Answer 3

这可能和你正在做的事情有点相似。我想看看这个公式：f(x,y) = x²/r1² + y²/r2² = 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) 这样的配对，这样可以进一步减少计算量。