1

I am doing an exercise to familiarize with Python least_squares from scipy.optimize.

The exercise try to fit an ellipse to a list of 2D points minimizing the sum of the square distances between the points and the ellipse.

Maybe the mathematical approach is not the right one but let's pretend it is fine because I think that my difficulties are elsewhere.

The idea is first to write a function that computes the distance between a point and an ellipse and then use this function in the optimizer.

I programmed also this distance function as a minimization problem: given a query point and a parametric equation of the ellipse, I seek for the minimum length line segment connecting the query point and a point belonging to the ellipse, its length is the desired distance.

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

Now I test ellipse_point_distance on a simple case:

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

and I get 0.414213562373, it seems good, so let's continue with the minimization problem for the fit:

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 get this error:

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

A quick look at TypeError: only length-1 arrays can be converted to Python scalars while trying to exponentially fit data and I thought I solved the problem:

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

I just replaced math.atan2 with np.arctan2 hoping for the best:

print(ellipse_point_distance_2(x_q,y_q,R_x,R_y,C_x,C_y,theta))

the ellipse_point_distance_2 is still good (gives 0.414213562373), so here we are:

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)

but now I get a different error: 

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.

now I am a little bit confused... I think my problem is related to a vectorization issue but I can not solve it.

Alessandro Jacopson
  • 18,047
  • 15
  • 98
  • 153
  • From scipy documents "x0 : array_like with shape (n,) or float" so I ask to you if the `x.mean()` and `y.mean()` used in `params_0 = np.array([227,227,x.mean(),y.mean(),0])` return scalars in a way that `param_0` would have the right shape – Hemerson Tacon Mar 05 '19 at 21:17
  • @HemersonTacon `print(params_0)` gives `[227. 227. 420.4281179 298.21243022 0. ]`. – Alessandro Jacopson Mar 05 '19 at 21:21
  • 1
    Just to make sure, what `print(params_0.shape)` returns? – Hemerson Tacon Mar 05 '19 at 21:23
  • You got the first problem right. `math.atan` takes scalar values (same for `math.cos`), while the `np.atan` works with array values. – hpaulj Mar 05 '19 at 21:46
  • @HemersonTacon `print(params_0.shape)` gives `(5,)` – Alessandro Jacopson Mar 05 '19 at 21:49
  • 1
    the error is in the inner `least_squares` call, the same that gave you problem with `atan`. Evidently one or more of the parameters passed to `ellipse_point_distance_2` are arrays, give problems to the `atan` call, and making that `params_0` 2d. – hpaulj Mar 05 '19 at 21:52

1 Answers1

0

In this function you need to change two lines:

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)])
    params_0 = np.array(math.atan2(y_q-C_y,x_q-C_x)) # removed inner square brackets
    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))
    d=np.sqrt(residual_for_distance(result.x,x_q,y_q,R_x,R_y,theta,C_x,C_y)) # changed from math.sqrt to np.sqrt
    return d

I think your code still doesn't work, but now it runs for me without errors. You may want to post another question if you have trouble getting least_squares to do what you want.

aganders3
  • 5,838
  • 26
  • 30