Ajuste de círculo de mínimos cuadrados utilizando MATLAB Optimization Toolbox

Question

Estoy tratando de implementar el ajuste de círculo de mínimos cuadrados después del papel this (lo siento, no puedo publicarlo). El documento establece que podríamos caber un círculo, calculando el error geométrico como la distancia euclidiana (Xi '') entre un punto específico (Xi) y el punto correspondiente del círculo (Xi '). Tenemos tres parámetros: Xc (un vector de coordenadas del centro del círculo) y R (radio).

me ocurrió con el siguiente código de MATLAB (en cuenta que estoy tratando de adaptarse a los círculos, y no las esferas como se indica en las imágenes):

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 

El accesorio sin embargo no es demasiado bueno: si comienzo con el buen vector de parámetros, el algoritmo termina en el primer paso (por lo que hay un mínimo local donde debería estar), pero si perturbo el punto de inicio (con un círculo silencioso) la conexión se detiene con errores muy grandes Estoy seguro de que he pasado por alto algo en mi implementación.

Answer 1

Por lo que vale, implementé estos métodos en MATLAB hace un tiempo. Sin embargo, lo hice aparentemente antes de saber sobre lsqnonlin, etc., ya que usa una regresión implementada a mano. Esto probablemente sea lento, pero puede ayudar a compararlo con su código.

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

Esto se ejecute con:

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

y se representa 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); 

Dar: