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

Question

我想实现以下最小二乘圆拟合这纸(抱歉,我不能发布).通过计算几何误差作为特定点(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);

赠送: