J'essaie d'implémenter l'ajustement de cercle des moindres carrés en suivant le papier this (désolé, je ne peux pas le publier). L'article indique que nous pourrions faire un cercle, en calculant l'erreur géométrique comme la distance euclidienne (Xi '') entre un point spécifique (Xi) et le point correspondant sur le cercle (Xi '). Nous avons trois paramètres: Xc (un vecteur de coordonnées au centre du cercle), et R (rayon).Ajustement du cercle des moindres carrés à l'aide de MATLAB Optimization Toolbox
je suis venu avec le code Matlab suivant (notez que je suis en train de tenir des cercles, et non des sphères comme il est indiqué sur les images):
function [ circle ] = fit_circle(X)
% Kör paraméterstruktúra inicializálása
% R - kör sugara
% Xc - kör középpontja
circle.R = NaN;
circle.Xc = [ NaN; NaN ];
% Kezdeti illesztés
% A köz középpontja legyen a súlypont
% A sugara legyen az átlagos négyzetes távolság a középponttól
circle.Xc = mean(X);
d = bsxfun(@minus, X, circle.Xc);
circle.R = mean(bsxfun(@hypot, d(:,1), d(:,2)));
circle.Xc = circle.Xc(1:2)+random('norm', 0, 1, size(circle.Xc));
% Optimalizáció
options = optimset('Jacobian', 'on');
out = lsqnonlin(@ort_error, [circle.Xc(1), circle.Xc(2), circle.R], [], [], options, X);
end
%% Cost function
function [ error, J ] = ort_error(P, X)
%% Calculate error
R = P(3);
a = P(1);
b = P(2);
d = bsxfun(@minus, X, P(1:2)); % X - Xc
n = bsxfun(@hypot, d(:,1), d(:,2)); % || X - Xc ||
res = d - R * bsxfun(@times,d,1./n);
error = zeros(2*size(X,1), 1);
error(1:2:2*size(X,1)) = res(:,1);
error(2:2:2*size(X,1)) = res(:,2);
%% Jacobian
xdR = d(:,1)./n;
ydR = d(:,2)./n;
xdx = bsxfun(@plus,-R./n+(d(:,1).^2*R)./n.^3,1);
ydy = bsxfun(@plus,-R./n+(d(:,2).^2*R)./n.^3,1);
xdy = (d(:,1).*d(:,2)*R)./n.^3;
ydx = xdy;
J = zeros(2*size(X,1), 3);
J(1:2:2*size(X,1),:) = [ xdR, xdx, xdy ];
J(2:2:2*size(X,1),:) = [ ydR, ydx, ydy ];
end
Le raccord cependant n'est pas trop bon: si je commence par le bon vecteur de param'etres, l'algorithme se termine'a la premi'ere étape (il y a donc un minimum local), mais si je perturbe le point de départ (avec un cercle silencieux) très grosses erreurs. Je suis sûr que j'ai oublié quelque chose dans ma mise en œuvre.