Fitting an implicit equation

Section 28.4 Fitting an implicit equation

This is an optional section, feel free to skip it.

Not all data points represent a function of the form \(y = f(x)\text{.}\) They may happen to accumulate along some closed curve such as an ellipse. General implicit equation of an ellipse is

\begin{equation} Ax^2 + Bxy + Cy^2 + Dx + Ey = 1\tag{28.4.1} \end{equation}

There are five unknown coefficients here. So, given a set of points \((x_k, y_k)\text{,}\) \(k=1, \dots, n\text{,}\) we get a linear system by plugging each point into (28.4.1). Its least squares solution describes an ellipse that fits the data best in the sense of having the smallest sum

\begin{equation*} \sum_{k=1}^n (Ax_k^2 + Bx_ky_k + Cy_k^2 + Dx_k + Ey_k - 1)^2 \end{equation*}

which should hopefully fit the ellipse.

🔗

Example 28.4.1. Fit an ellipse to given points.

Plot an ellipse that fits the given 7 data points:

  
x = [1 2 4 5 4 3 1]'; 
y = [2 0 -1 1 2 3 4]';

Then plot it together with the points, using fimplicit command. Its syntax is

fimplicit(@(x, y) ..., [xmin xmax ymin ymax])

where the first argument is the implicit equation (something equated to 0), and the second is the plot window.

🔗

Solution.

X = [x.^2 x.*y y.^2 x y];
b = X\ones(7, 1);
hold on
plot(x, y, 'r*')
fun = @(x, y) b(1)*x.^2 + b(2)*x.*y + b(3)*y.^2 + b(4)*x + b(5)*y - 1;
win = [min(x)-1 max(x)+1 min(y)-1 max(y)+1]; 
fimplicit(fun, win)
hold off

Here the vector of parameters is called b instead of beta to shorten the implicit formula. It is obtained by solving a system with 1 on the right hand side, according to (28.4.1). The plot window is calculated from the data to allow margins of 1 unit on all sides.

🔗

Note that while it’s safe to leave spaces around + signs in the formula fun (some find it improves readability), one has to be extra careful with such spaces within a vector. Indeed,

[min(x) -1 max(x) +1 min(y) -1 max(y) +1]

would be a vector of 8 numbers, not 4. Matlab does understand

[min(x) - 1 max(x) + 1 min(y) - 1 max(y) + 1];

as a vector of 4 numbers, however.

🔗

Question 28.4.2. Vectorizing a function for implicit plot.

It is tempting to write the anonymous function as fun = @(x, y) [x.^2 x.*y y.^2 x y]*b - 1 which is shorter and consistent with how we used parameters in the past. And this works but Matlab complains: “Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.” What is it complaining about?

🔗

Prev Top Next