使用 MATLAB 优化工具箱进行最小二乘圆拟合

Question

我正在尝试在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

但是拟合不太好：如果我从好的参数向量开始，算法会在第一步终止（所以应该有一个局部最小值），但是如果我扰乱了起点（无噪音圈）拟合以非常大的错误停止。我确信我在实现中忽略了一些东西。

Answer 1

不管怎样，我不久前在 MATLAB 中实现了这些方法。但是，我显然在知道lsqnonlin 等之前就这样做了，因为它使用了手动实现的回归。这可能很慢，但可能有助于与您的代码进行比较。

function [x, y, r, sq_error] = circFit ( P )
%# CIRCFIT fits a circle to a set of points using least sqaures
%#  P is a 2 x n matrix of points to be fitted

per_error = 0.1/100; % i.e. 0.1%

%# initial estimates
X  = mean(P, 2)';
r = sqrt(mean(sum((repmat(X', [1, length(P)]) - P).^2)));

v_cen2points = zeros(size(P));
niter = 0;

%# looping until convergence
while niter < 1 || per_diff > per_error

    %# vector from centre to each point
    v_cen2points(1, :) = P(1, :) - X(1);
    v_cen2points(2, :) = P(2, :) - X(2);  

    %# distacnes from centre to each point
    centre2points = sqrt(sum(v_cen2points.^2));

    %# distances from edge of circle to each point
    d = centre2points - r;

    %# computing 3x3 jacobean matrix J, and solvign matrix eqn.
    R = (v_cen2points ./ [centre2points; centre2points])';
    J = [ -ones(length(R), 1), -R ];
    D_rXY = -J\d';

    %# updating centre and radius
    r_old = r;    X_old = X;
    r = r + D_rXY(1);
    X = X + D_rXY(2:3)';

    %# calculating maximum percentage change in values
    per_diff = max(abs( [(r_old - r) / r, (X_old - X) ./ X ])) * 100;

    %# prevent endless looping
    niter = niter + 1;
    if niter > 1000
        error('Convergence not met in 1000 iterations!')
    end
end

x = X(1);
y = X(2);
sq_error = sum(d.^2);

然后运行：

X = [1 2 5 7 9 3];
Y = [7 6 8 7 5 7];
[x_centre, y_centre, r] = circFit( [X; Y] )

并绘制：

[X, Y] = cylinder(r, 100);
scatter(X, Y, 60, '+r'); axis equal
hold on
plot(X(1, :) + x_centre, Y(1, :) + y_centre, '-b', 'LineWidth', 1);

给予：