Florie - 2 months ago 8x
Python Question

# Counting points inside an ellipse

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?

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

'''
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:
ellips = xx/(rx*rx)+yy/(ry*ry)
if ellips < 1:
ithEllipse = count
break
count+=1
ithEllipse = None
break

return ithEllipse

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

# plot Ellipses
ellipse = patches.Ellipse((0,0),rx*2,ry*2,fc='none',ec='red')

x=3.0
y=1.0
ellipse = patches.Ellipse((0,0),rx*2,ry*2,fc='none',ec='blue')

if idx != None:
circle = patches.Circle((x,y),.1)
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.