Contando puntos dentro de una elipse.

Estoy tratando de contar los puntos de datos dados dentro de cada anillo de elipse:

introduzca la descripción de la imagen aquí

El problema es que tengo una función para verificar que: así, para cada elipse, para asegurarme de que haya un punto en él, se deben calcular tres entradas:

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 

Ahora mismo estoy comprobando si está o no en una elipse por:

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

Pero luego tendría que calcular cada par de puntos focales de los bucles externos … ¿hay alguna forma más eficiente e inteligente de hacer esto?

Esto puede ser algo similar a lo que estás haciendo. Solo estoy mirando para ver si f (x, y) = x ^ 2 / r1 ^ 2 + y ^ 2 / r2 ^ 2 = 1.

Cuando f (x, y) es mayor que 1, el punto x, y está fuera de la elipse. Cuando es más pequeño, entonces está dentro de la elipse. Recorro cada elipse para encontrar el que f (x, y) es menor que 1.

El código tampoco tiene en cuenta una elipse que está centrada en el origen. Es un pequeño cambio para incluir esta característica.

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

Este código produce la siguiente figura: introduzca la descripción de la imagen aquí

Ten en cuenta que esto es solo un punto de partida. Por un lado, puede cambiar en la que se acepte una lista de los cuadrados de r1 y r2, es decir (r1 * r1, r2 * r2) pares, y eso reduciría aún más el cálculo.

Se complican las cosas. No es necesario calcular los puntos focales y las distancias a los puntos focales, etc. de acuerdo con la definición geométrica de la elipse. Si conoce los ejes mayor y menor (lo hace), simplemente exprima un poco la pregunta completa (de modo que ambos sean 1.0, por ejemplo, dividiendo x-centerx y y-centery por xaxis y yaxis) y luego la pregunta de si el punto está dentro de la elipse es simplemente

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

PD: En general, un buen consejo en este campo: no hacer sqrt si puede hacer lo mismo al no calcular una distancia, sino permanecer cómodamente en el ámbito de su cuadrado.

Aquí hay algunas ideas para usted:

  • Tienes la idea correcta de mover el código para calcular los focos fuera del bucle.
  • Los cálculos de distancia se pueden acelerar eliminando las raíces cuadradas. En otras palabras, sabemos que a < b implica sqrt(a) < sqrt(b) por lo que no es necesario calcular la raíz cuadrada.
  • Si las elipsis son concéntricas y el eje mayor es paralelo al eje x, puede simplificar el problema de la elipse a un problema de círculo volviendo a escalar el valor de x.

Además, aquí hay una pequeña nit de encoding. No es necesario que una sentencia if devuelva Verdadero o Falso . En su lugar, puede devolver la expresión condicional en sí:

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