Próbuję wprowadzić dopasowanie koła najmniejszych kwadratów po papierze this (przepraszam, nie mogę go opublikować). Artykuł stwierdza, że możemy dopasować okrąg, obliczając błąd geometryczny jako odległość euklidesową (Xi '') między określonym punktem (Xi) i odpowiednim punktem na okręgu (Xi '). Mamy trzy parametry: Xc (wektor współrzędnych środka koła) i R (promień).Dopasowywanie najmniejszego kwadratu koła przy użyciu narzędzia MATLAB Optymalizacja Toolbox
wpadłem na następujący kod MATLAB (zauważ, że staram się zmieścić koła, a nie sfery jak to jest wskazane na zdjęciach):
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
Oprawa jednak nie jest zbyt dobra: jeśli zacznę od dobrego wektora parametrów, algorytm kończy się w pierwszym kroku (więc tam gdzie powinien być lokalny minimalny), ale jeśli zaburzę punkt początkowy (z bezszumowym kółkiem), łącznik zatrzymuje się z bardzo duże błędy. Jestem pewien, że przeoczyłem coś w mojej implementacji.