MATLAB अनुकूलन टूलबॉक्स

Question

का उपयोग कर कम वर्ग सर्कल फिटिंग मैं this पेपर के बाद कम से कम वर्ग सर्कल फिटिंग को लागू करने की कोशिश कर रहा हूं (क्षमा करें मैं इसे प्रकाशित नहीं कर सकता)। पेपर कहता है, कि हम एक विशिष्ट बिंदु (Xi) और सर्कल (Xi ') के इसी बिंदु के बीच euclidean दूरी (Xi') के रूप में ज्यामितीय त्रुटि की गणना करके, एक सर्कल फिट कर सकते हैं। हमारे पास तीन पैरामीटर हैं: एक्ससी (सर्कल के केंद्र का समन्वय करने वाला एक वेक्टर), और आर (त्रिज्या)। तथापि

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 

फिटिंग

:

मैं निम्नलिखित MATLAB कोड (ध्यान दें कि मैं हलकों, क्षेत्रों नहीं फिट करने के लिए के रूप में यह छवियों से दिखाया जाता है कोशिश कर रहा हूँ) के साथ आया था बहुत अच्छा नहीं है: अगर मैं अच्छे पैरामीटर वेक्टर से शुरू करता हूं तो एल्गोरिदम पहले चरण में समाप्त हो जाता है (इसलिए एक स्थानीय मिनीमा जहां यह होना चाहिए), लेकिन अगर मैं शुरुआती बिंदु (एक निर्बाध सर्कल के साथ) को परेशान करता हूं तो फिटिंग बंद हो जाती है बहुत बड़ी त्रुटियां मुझे यकीन है कि मैंने अपने कार्यान्वयन में कुछ अनदेखा कर दिया है।

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

: