Counting points inside an ellipse

Question

I'm trying to count given data points inside each ring of ellipse:

The problem is that I have a function to check that: so for each ellipse, to make sure whether a point is in it, three inputs have to be calculated:

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

Right now I'm checking whether or not it's in an ellipse by:

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

But I would then have to calculate each pairs of focal points of the outside loops.. is there any more efficient and or clever way to do this?

Answer 1

You complicate things. It is no need to compute focal points and the distances to the focal points etc. according to the geometric definition of ellipse. If you know major and minor axis (you do), just squeeze the whole question a bit (so that both are 1.0, for example, by dividing x-centerx and y-centery by xaxis and yaxis) and then the question whether the point is inside ellipse is simply

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

P.S.: In general, good advice in this field: no sqrt if you can do the same thing by not actually compute a distance but staying comfortably in the realm of its square.

Answer 2

This may be something similar to what you are doing. I'm just looking to see if f(x,y) = x^2/r1^2 + y^2/r2^2 = 1.

When f(x,y) is larger than 1, the point x,y is outside the ellipse. When it is smaller, then it is inside the ellipse. I loop through each ellipse to find the one when f(x,y) is smaller than 1.

The code also does not take into account an ellipse that is centered off the origin. It's a small change to include this feature.

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

This code produces the following figure:

Keep in mind, this is just a starting point. For one thing, you can change inWhichEllipse to accept a list of the square of r1 and r2, ie (r1*r1,r2*r2) pairs, and that would cut the computation down even more.

Answer 3

Here are a few ideas for you:

• You've got the right idea moving the code for computing the foci outside of the loop.
• The distance calculations can be sped-up by removing the square roots. In other words, we know a < b implies sqrt(a) < sqrt(b) so there is no need to calculate the square root.
• If the ellipses are concentric and the major axis is parallel to the x-axis, you can simplify the ellipse problem to a circle problem by rescaling the x value.

Also, here's one minor coding nit. There is no need for an if-statement to return True or False. Instead, you can return the conditional expression itself:

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