Implicit orthogonal regression¶

Estimate the parameters of an ellipse from a set of coordinates.

$$\begin{align*} f(\bm{X}, \bm{\beta}) & = \frac{\left[(x-x_0)\cos\theta + (y-y_0)\sin\theta\right]^2}{a^2} \\ & + \frac{\left[(y-y_0)\cos\theta -(x-x_0)\sin\theta\right]^2}{b^2} - 1 = 0 \end{align*}$$

$$ \bm{X} = (x,y) $$

$$ \bm{\beta} = (x_0, y_0, a, b, \theta) $$

In [1]:

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Ellipse

from odrpack import odr_fit

import matplotlib.pyplot as plt import numpy as np from matplotlib.patches import Ellipse from odrpack import odr_fit

First, we define the observed data and the model function.

In [2]:

# Each row represents a point (x, y) of the ellipse
X = [[0.50, -0.12],
     [1.20, -0.60],
     [1.60, -1.00],
     [1.86, -1.40],
     [2.12, -2.54],
     [2.36, -3.36],
     [2.44, -4.00],
     [2.36, -4.75],
     [2.06, -5.25],
     [1.74, -5.64],
     [1.34, -5.97],
     [0.90, -6.32],
     [-0.28, -6.44],
     [-0.78, -6.44],
     [-1.36, -6.41],
     [-1.90, -6.25],
     [-2.50, -5.88],
     [-2.88, -5.50],
     [-3.18, -5.24],
     [-3.44, -4.86]]

# We need shape (2, n)
Xdata = np.array(X).T

# Ydata is not used, but is required
Ydata = np.zeros(Xdata.shape[-1])

# Each row represents a point (x, y) of the ellipse X = [[0.50, -0.12], [1.20, -0.60], [1.60, -1.00], [1.86, -1.40], [2.12, -2.54], [2.36, -3.36], [2.44, -4.00], [2.36, -4.75], [2.06, -5.25], [1.74, -5.64], [1.34, -5.97], [0.90, -6.32], [-0.28, -6.44], [-0.78, -6.44], [-1.36, -6.41], [-1.90, -6.25], [-2.50, -5.88], [-2.88, -5.50], [-3.18, -5.24], [-3.44, -4.86]] # We need shape (2, n) Xdata = np.array(X).T # Ydata is not used, but is required Ydata = np.zeros(Xdata.shape[-1])

In [3]:

def f(X: np.ndarray, beta: np.ndarray) -> np.ndarray:
    x0, y0, a, b, theta = beta
    x, y = X
    return ((x - x0)*np.cos(theta) + (y - y0)*np.sin(theta))**2 / a**2 \
         + ((y - y0)*np.cos(theta) - (x - x0)*np.sin(theta))**2 / b**2 \
         - 1

def f(X: np.ndarray, beta: np.ndarray) -> np.ndarray: x0, y0, a, b, theta = beta x, y = X return ((x - x0)*np.cos(theta) + (y - y0)*np.sin(theta))**2 / a**2 \ + ((y - y0)*np.cos(theta) - (x - x0)*np.sin(theta))**2 / b**2 \ - 1

Then, we define a plausible initial guess $\bm{\beta_0}$ for the model parameters, as well as the corresponding bounds.

In [4]:

beta0 = np.array([0.0, 0.0, 1.0, 1.0, 0.0])

lower = np.array([-1e2, -1e2, 0e0, 0e0, -np.pi/2])
upper = np.array([+1e2, +1e2, 1e2, 1e2, +np.pi/2])

beta0 = np.array([0.0, 0.0, 1.0, 1.0, 0.0]) lower = np.array([-1e2, -1e2, 0e0, 0e0, -np.pi/2]) upper = np.array([+1e2, +1e2, 1e2, 1e2, +np.pi/2])

Here, we expect the measurement error to be the same across both $\bm{X}$ coordinates, so a special weighting scheme is unnecessary.

In [5]:

weight_x = 1.0

weight_x = 1.0

We can now launch the regression! As the problem is implicit, we set task='implicit-ODR'. If you want to see a brief computation report, set report=short.

In [6]:

sol = odr_fit(f, Xdata, Ydata, beta0,
              bounds=(lower, upper),
              weight_x=weight_x,
              task='implicit-ODR')

sol = odr_fit(f, Xdata, Ydata, beta0, bounds=(lower, upper), weight_x=weight_x, task='implicit-ODR')

The result is packed in a OdrResult dataclass. Let's check the solution convergence and the estimated model parameters.

In [7]:

sol.stopreason

sol.stopreason

Out[7]:

'Parameter convergence.'

In [8]:

sol.beta

sol.beta

Out[8]:

array([-0.99938103, -2.93104802,  3.86422581,  3.15663797, -0.90359144])

All fine! Let's plot the solution.

In [9]:

fig, ax = plt.subplots()

# Plot observed data
ax.plot(*Xdata, 'o')

# Plot fitted ellipse
x0, y0, a, b, theta = sol.beta
ellipse = Ellipse((x0, y0), width=2*a, height=2*b, angle=np.degrees(theta),
                  facecolor='none', edgecolor='orange', linewidth=2)
ax.add_patch(ellipse)
ax.set_xlabel('x')
ax.set_ylabel('y')

fig, ax = plt.subplots() # Plot observed data ax.plot(*Xdata, 'o') # Plot fitted ellipse x0, y0, a, b, theta = sol.beta ellipse = Ellipse((x0, y0), width=2*a, height=2*b, angle=np.degrees(theta), facecolor='none', edgecolor='orange', linewidth=2) ax.add_patch(ellipse) ax.set_xlabel('x') ax.set_ylabel('y')

Out[9]:

Text(0, 0.5, 'y')