Vòng tròn hình vuông nhỏ nhất phù hợp bằng cách sử dụng Hộp công cụ tối ưu hóa MATLAB

Question

Tôi đang cố gắng triển khai tối thiểu hình vuông phù hợp theo sau this giấy (xin lỗi tôi không thể xuất bản nó). Các trạng thái giấy, rằng chúng ta có thể phù hợp với một vòng tròn, bằng cách tính toán sai số hình học như khoảng cách euclide (Xi '') giữa một điểm cụ thể (Xi) và điểm tương ứng trên vòng tròn (Xi '). Chúng ta có ba parametres: Xc (vectơ tọa độ tâm của hình tròn), và R (bán kính).

tôi đã đưa ra đoạn mã sau MATLAB (lưu ý rằng tôi đang cố gắng để phù hợp với vòng tròn, không lĩnh vực như nó được ghi trên hình ảnh):

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 

Các phụ kiện tuy nhiên không phải là quá tốt: nếu tôi bắt đầu với vector tham số tốt thuật toán kết thúc ở bước đầu tiên (vì vậy có một minima địa phương nơi nó nên), nhưng nếu tôi làm rối loạn điểm khởi đầu (với vòng tròn không ồn) các lỗi rất lớn. Tôi chắc chắn rằng tôi đã bỏ qua một cái gì đó trong việc thực hiện của tôi.

Answer 1

Đối với những gì đáng giá, tôi đã triển khai các phương pháp này trong MATLAB một thời gian trước đây. Tuy nhiên, tôi đã làm điều đó rõ ràng trước khi tôi biết về lsqnonlin vv, vì nó sử dụng một hồi quy được thực hiện bằng tay. Điều này có thể chậm, nhưng có thể giúp so sánh với mã của bạn.

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

này sau đó chạy với:

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

Và âm mưu với:

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

Giving: