Подгонка окружности по наименьшим квадратам с помощью MATLAB Optimization Toolbox

Я пытаюсь реализовать подгонку круга по методу наименьших квадратов, следуя этой статье (жаль, что не могу ее опубликовать). В статье говорится, что мы можем подогнать окружность, вычислив геометрическую ошибку как евклидово расстояние (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

Однако подгонка не слишком хороша: если я начинаю с хорошего вектора параметров, алгоритм завершается на первом шаге (так что есть локальный минимум там, где он должен быть), но если я возмущаю начальную точку (с бесшумной окружностью), подгонка останавливается с очень большими ошибками. Я уверен, что я упустил что-то в своей реализации.