Я делаю упражнение, чтобы ознакомиться с Python least_squares
из scipy.optimize
.
В упражнении попробуйте подогнать эллипс к списку 2D точек, минимизируя сумму квадратов расстояний между точками иэллипс.
Возможно, математический подход не верен, но давайте представим, что это хорошо, потому что я думаю, что мои трудности в другом месте.
Идея состоит в том, чтобы сначала написать функцию, которая вычисляетрасстояние между точкой и эллипсом, а затем использовать эту функцию в оптимизаторе.
Я запрограммировал также эту функцию расстояния как задачу минимизации: учитывая точку запроса и параметрическое уравнение эллипса, я ищу минимумДлина отрезка линии, соединяющего точку запроса и точку, принадлежащую эллипсу, его длина - это желаемое расстояние.
import math
import numpy as np
from scipy.optimize import least_squares
# I would like to fit an ellipse to these points (a point in each row):
p=[
[614.0471123474172,289.51195416538405],
[404.85868232180786,509.3183970173126],
[166.5322099316754,335.6006010213824],
[302.6076456817051,71.14357043842081],
[614.094939200562,285.48762845572804]
]
# This is the x of the parametric equation of an ellipse
# centered at (C_x,C_y), with axis R_x and R_y and angle
# of rotation theta. alpha is the parameter that describe
# the ellipse when going from 0 to 2pi.
def x_e(alpha,R_x,R_y,theta,C_x):
return R_x*math.cos(alpha)*math.cos(theta)-R_y*math.sin(alpha)*math.sin(theta)+C_x
# This is the y
def y_e(alpha,R_x,R_y,theta,C_y):
return R_x*math.cos(alpha)*math.sin(theta)+R_y*math.sin(alpha)*math.cos(theta)+C_y
points = np.array(p)
x=points[:,0]
y=points[:,1]
def residual_for_distance(params,x_q,y_q,R_x,R_y,theta,C_x,C_y):
alpha = params[0]
return (x_q-x_e(alpha,R_x,R_y,theta,C_x))**2+(y_q-y_e(alpha,R_x,R_y,theta,C_y))**2
def ellipse_point_distance(x_q,y_q,R_x,R_y,C_x,C_y,theta):
params_0 = np.array([math.atan2(y_q-C_y,x_q-C_x)])
result = least_squares(residual_for_distance,params_0,args=(x_q,y_q,R_x,R_y,theta,C_x,C_y))
d=math.sqrt(residual_for_distance(result.x,x_q,y_q,R_x,R_y,theta,C_x,C_y))
return d
Теперь я проверяю ellipse_point_distance
на простом случае:
x_q=1
y_q=1
R_x=1
R_y=1
C_x=0
C_y=0
theta=0
print(ellipse_point_distance(x_q,y_q,R_x,R_y,C_x,C_y,theta))
и я получаю 0.414213562373
, это кажется хорошим, поэтому давайте продолжим с проблемой минимизации для подгонки:
def residual_for_fit(params,x,y):
R_x = params[0]
R_y = params[1]
C_x = params[2]
C_y = params[3]
theta = params[4]
return ellipse_point_distance(x,y,R_x,R_y,C_x,C_y,theta)
params_0 = np.array([227,227,x.mean(),y.mean(),0])
result = least_squares(residual_for_fit,params_0,args=(x,y),verbose=1)
I gи эта ошибка:
Traceback (most recent call last):
File "fit_ellipse.py", line 57, in <module>
result = least_squares(residual_for_fit,params_0,args=(x,y),verbose=1)
File "/home/aj/anaconda2/lib/python2.7/site-packages/scipy/optimize/_lsq/least_squares.py", line 799, in least_squares
f0 = fun_wrapped(x0)
File "/home/aj/anaconda2/lib/python2.7/site-packages/scipy/optimize/_lsq/least_squares.py", line 794, in fun_wrapped
return np.atleast_1d(fun(x, *args, **kwargs))
File "fit_ellipse.py", line 54, in residual_for_fit
return ellipse_point_distance(x,y,R_x,R_y,C_x,C_y,theta)
File "fit_ellipse.py", line 33, in ellipse_point_distance
params_0 = np.array([math.atan2(y_q-C_y,x_q-C_x)])
TypeError: only size-1 arrays can be converted to Python scalars
Краткий обзор TypeError: только массивы длины 1 могут быть преобразованы в скаляры Python при попытке экспоненциального соответствия данных , и я решил, что решил проблему:
def ellipse_point_distance_2(x_q,y_q,R_x,R_y,C_x,C_y,theta):
params_0 = np.array([np.arctan2(y_q-C_y,x_q-C_x)])
result = least_squares(residual_for_distance,params_0,args=(x_q,y_q,R_x,R_y,theta,C_x,C_y))
d=math.sqrt(residual_for_distance(result.x,x_q,y_q,R_x,R_y,theta,C_x,C_y))
return d
Я только что заменил math.atan2
на np.arctan2
в надежде на лучшее:
print(ellipse_point_distance_2(x_q,y_q,R_x,R_y,C_x,C_y,theta))
ellipse_point_distance_2
все еще хорош (дает 0.414213562373
), поэтому здесь мыявляются:
def residual_for_fit_2(params,x,y):
R_x = params[0]
R_y = params[1]
C_x = params[2]
C_y = params[3]
theta = params[4]
return ellipse_point_distance_2(x,y,R_x,R_y,C_x,C_y,theta)
params_0 = np.array([227,227,x.mean(),y.mean(),0])
result = least_squares(residual_for_fit_2,params_0,args=(x,y),verbose=1)
но теперь я получаю другую ошибку:
Traceback (most recent call last):
File "fit_ellipse.py", line 76, in <module>
result = least_squares(residual_for_fit_2,params_0,args=(x,y),verbose=1)
File "/home/aj/anaconda2/lib/python2.7/site-packages/scipy/optimize/_lsq/least_squares.py", line 799, in least_squares
f0 = fun_wrapped(x0)
File "/home/aj/anaconda2/lib/python2.7/site-packages/scipy/optimize/_lsq/least_squares.py", line 794, in fun_wrapped
return np.atleast_1d(fun(x, *args, **kwargs))
File "fit_ellipse.py", line 73, in residual_for_fit_2
return ellipse_point_distance_2(x,y,R_x,R_y,C_x,C_y,theta)
File "fit_ellipse.py", line 61, in ellipse_point_distance_2
result = least_squares(residual_for_distance,params_0,args=(x_q,y_q,R_x,R_y,theta,C_x,C_y))
File "/home/aj/anaconda2/lib/python2.7/site-packages/scipy/optimize/_lsq/least_squares.py", line 772, in least_squares
raise ValueError("`x0` must have at most 1 dimension.")
ValueError: `x0` must have at most 1 dimension.
теперь я немного запутался ... Я думаю, что моя проблема связана с проблемой векторизации, но яне могу решить это.