我正在尝试实现最小二乘圆拟合,参考 this 论文(很抱歉我无法公开它)。该论文指出,我们可以通过计算几何误差作为欧几里得距离(Xi'')来拟合一个圆,其中特定点(Xi)与圆上相应点(Xi')之间的距离。我们有三个参数:Xc(圆心坐标向量)和R(半径)。
然而,拟合并不太好:如果我从良好的参数向量开始,算法在第一步终止(因此存在应该在其中的局部最小值),但是如果我扰动起始点(用一个无噪声的圆),拟合会停止,产生非常大的误差。我确定我在实现中忽略了某些东西。
我想到了以下的MATLAB代码(请注意,我尝试拟合的是圆形,而不是图片上标注的球体):
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
然而,拟合并不太好:如果我从良好的参数向量开始,算法在第一步终止(因此存在应该在其中的局部最小值),但是如果我扰动起始点(用一个无噪声的圆),拟合会停止,产生非常大的误差。我确定我在实现中忽略了某些东西。