Raccordo del cerchio dei minimi quadrati utilizzando Toolbox Tools Ottimizzazione Matlab

Question

Sto cercando di implementare il raccordo del cerchio dei minimi quadrati seguente Questa carta (scusa, non posso pubblicarlo). Gli stati cartacee, che potremmo montare un cerchio, calcolando l'errore geometrico come la distanza euclidea (XI ") tra un punto specifico (XI) e il punto corrispondente sul cerchio (XI"). Abbiamo tre parametri: XC (un vettore di coordina il centro del cerchio), e r (raggio).

Ho inventato il seguente codice MATLAB (nota che sto cercando di adattarsi a cerchi, non sfere come è indicato sulle immagini):

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

.

Il raccordo tuttavia non è troppo buono: se inizio con il buon parametro il vettore l'algoritmo termina al primo passo (quindi c'è un minima locale dove dovrebbe essere), ma se perturgo il punto di partenza (con un nobile cerchio) Il raccordo si ferma con errori molto grandi. Sono sicuro che ho trascurato qualcosa nella mia implementazione.

Answer 1

Per quello che vale, ho implementato questi metodi a Matlab a tempo fa.Tuttavia, l'ho fatto apparentemente prima che tu sapessi del lsqnonlin ecc., In quanto utilizza una regressione implementata a mano.Questo è probabilmente lento, ma può aiutare a confrontare con il tuo codice.

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);

.

Questo viene quindi eseguito con:

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

.

e tracciato con:

[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);

.

Dare: