Approximation d'un cercle avec la méthode des moindres carrés

Hypothèses

Cette page explique comment ajuster un cercle 2D à un nuage de points en minimisant les erreurs des moindres carrés.

Considérons un nuage de $n$ points ayant pour coordonnées $ {x_i,y_i} $. Le but est d'estimer $ x_c $ , $ y_c $ et $ r $, les paramètres du cercle qui épouse au mieux le nuage de points :

$ x_c $ est la coordonnées en x du centre du cercle
$ y_c $ est la coordonnées en y du centre du cercle
$ r $ est le rayon du cercle

L'équation du cercle idéal est donnée par :

$$ (x_i-x_c)^2 + (y_i-y_c)^2 = r^2 $$

Système matriciel

L'équation précédente peut être réécrite sous la forme :

$$ x_i^2 + x_c^2 - 2x_ix_c + y_i^2 + y_c^2 - 2y_iy_c = r^2 $$

ou encore :

$$ 2x_cx_i + 2y_cy_i + r^2 - x_c^2 - y_c^2 = x_i^2 + y_i^2 $$

et finalement :

$$ ax_i + by_i + c = x_i^2 + y_i^2 $$

avec :

$ a = 2x_c$
$ b = 2y_c$
$ c = r^2 - x_c^2 - y_c^2 $

Le système complet (pour tous les points) peut donc s'écrire :

$$ \begin{array}{rcl} ax_1 + by_1 + c & = & x_1^2 + y_1^2 \\ ax_2 + by_2 + c & = & x_2^2 + y_2^2 \\ & ... & \\ ax_n + by_n + c & = & x_n^2 + y_n^2 \\ \end{array} $$

La formulation matricielle du système est donnée par :

$$ \left[ \begin{matrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ ... & ... & ... \\ x_n & y_n & 1 \\ \end{matrix} \right]. \left[ \begin{matrix} a \\ b \\ c \end{matrix} \right] = \left[ \begin{matrix} x_1^2 + y_1^2 \\ x_2^2 + y_2^2 \\ ... \\ x_n^2 + y_n^2 \\ \end{matrix} \right] $$

Définissions $A$, $B$ and $X$:

$$ \begin{matrix} A=\left[ \begin{matrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ ... & ... & ... \\ x_n & y_n & 1 \\ \end{matrix} \right] & B=\left[ \begin{matrix} x_1^2 + y_1^2 \\ x_2^2 + y_2^2 \\... \\ x_n^2 + y_n^2 \\ \end{matrix} \right] & X=\left[ \begin{matrix} a \\ b \\ c \end{matrix} \right] \end{matrix} $$

Le système est maintenant donné par :

$$ A.X=B $$

Résolution du système

La solution optimale est donnée par :

$$ \hat{x}=A^{+}.B = (A^{T}.A)^{-1}A^{T}.B $$

Où $ A^{+} $ est la pseudoinverse de $ A $. $ A^{+} $ peut être calculée avec la formule suivante (Moore-Penrose inverse) :

$$ A^{+}=(A^{T}.A)^{-1}A^{T} $$

Comme nous avons effectué un changement de variable, il ne nous reste plus quà calculer $ x_c $ , $ y_c $ and $ r $ :

$$ x_c = \dfrac{a}{2} $$ $$ y_c = \dfrac{b}{2} $$ $$ r = \dfrac{ \sqrt{4c + a^2 + b^2 }}{2} $$

et nous obtenons finalement les paramètres du cercle :

Code source

Le code source MatLab ci-dessous a été utilisé pour générer la figure présentée plus haut :

close all;
clear all;
clc;

%% Parameters of the system
%% This is the parameters we need to estimate
% Center
C = [3 , 2];
% Radius
R = 1;

% Number of points
N = 1000;
% Noise
k = 0.050;

%% Create noisy circle
alpha = 2*pi*rand(1,N);
noise = k*2*randn(1,N)-1;
Points = C + [ R*noise.*cos(alpha) ; R*noise.*sin(alpha) ]';

%% Draw the noisy circle
plot (Points(:,1) , Points(:,2), 'b.');
axis square equal;
grid on;

%% Prepare matrices
A = [Points(:,1) Points(:,2) ones(N,1)];
B = [Points(:,1).*Points(:,1) + Points(:,2).*Points(:,2)];

%% Least square approximation
X=pinv(A)*B;
% or:
%X=inv(A'*A)*A'*B;

%% Calculate circle parameter
xc = X(1)/2
yc = X(2)/2
r = sqrt(4*X(3) + X(1)*X(1) + X(2)*X(2) )/2

%% Draw the fitted circle
Circle = [ xc + r*cos([0:pi/50:2*pi]) ; yc + r*sin([0:pi/50:2*pi]) ]';
hold on;
plot (Circle(:,1), Circle(:,2), 'r', 'LineWidth', 2);
legend ('Points', 'Least-square circle');

Téléchargement

Le code source MatLab ci-dessous a été utilisé pour générer la figure présentée plus haut :

least-squares-circle-fitting.m

Voir aussi

Dernière mise à jour : 22/11/2023