Нелинейный МНК для круговой сигмоидальной функции (CSF-LS)

Question

Я пытаюсь реализовать нелинейный метод наименьших квадратов для круговой сигмоидальной функции, чтобы уточнить положение центра окружности, описываемой массивом точек - пикселей изображения. Алгоритм беру из статьи https://www.google.ru/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&ved=2ahUKEwjms8G5v4j0AhXqAhAIHVcJCPoQFnoECCkQAQ&url=https%3A%2F%2Foaktrust.library.tamu.edu%2Fbitstream%2Fhandle%2F1969.1%2F191564%2FBORISSOV-DISSERTATION-2020.pdf%3Fsequence%3D1%26isAllowed%3Dy&usg=AOvVaw1xddaQRJ7uqUkdAOrx33dF (начиная со страницы 72). Мне нужно найти окружность, поэтому якобиан я ищу эксцентриситета. Также, для критерия останова я использую следующее правило: Моя реализация выглядит следующим образом:

function [r0,c0,r_e] = csf_ls(G,Lights,x0)
%%    [ix,iy] = find(G~=0);
%%    data = [iy';ix']';
%%    Lights = [];
%%    for i = 1:length(data(:,1))
%%        for j = 1:length(data(:,2))
%%            if [data(i,1),data(j,2)]==data(i,:)
%%                %disp([num2str(G(data(j,2),data(i,1))) ',' num2str(data(i,1)) ',' num2str(data(j,2))]);
%%                Lights=[Lights; G(data(j,2),data(i,1)),data(i,1),data(j,2)];
%%            end
%%        end
%%    end 
    Lights = unique(Lights,'rows');
    eps = 0.001;
    res = Inf;
    iter = 0;
    t = 0:pi/50:2*pi; 
    cond1 = 1;
    cond2 = 1;
    x0 = double(x0);
    r    = Lights(:,2);
    c    = Lights(:,3);
    fun = @(x,points)x(1)+(x(2)-x(1))/(1+exp(x(3)*(x(4)-sqrt((x(6)-points(:,2)).^2+(x(5)-points(:,1)).^2))));
    Niter = 50;
    while cond1*cond2==1
        iter = iter+1;
        
        r0   = x0(5);
        c0   = x0(6);
        r_e  = x0(4);
        k    = x0(3);
        ymax = x0(1);
        ymin = x0(2);
        d = sqrt((c0-c).^2+(r0-r).^2);
        alp = exp(k*(r_e-d));
        df_dr0 = k*(r0-r).*alp.*(ymax-ymin)./(d.*(1+alp).^2);
        df_dc0 = k*(c0-c).*alp.*(ymax-ymin)./(d.*(1+alp).^2);
        df_dre = -(ymax-ymin)*alp*k./((1+alp).^2);
        df_dk = -(ymin-ymax)*alp.*(r_e-d)./((1+alp).^2);
        df_dymin = -alp./(1+alp);
        df_dymax = 1./(1+alp);
        J = [df_dymin, df_dymax, df_dk, df_dre, df_dr0, df_dc0];
        f = fun(x0,[Lights(:,2),Lights(:,3)]);
        
        res = -inv(J'*J)*J'*(f' - Lights(:,1));
        x0 = x0+res';
%%        plot(x0(1),x0(2),'c+')
%%        plot(x0(1)+x0(3)*cos(t),x0(2)+x0(3)*sin(t),'c-')
        disp([num2str(x0(5)) ',' num2str(x0(6)) ',' num2str(x0(4))])
        disp(['iter = ' num2str(iter)])
        
        dXk = sqrt(abs(res(5)^2+res(6)^2-res(4)^2))
        cond1 = dXk>=eps/Niter
        cond2 = iter<50;
    end
    
    r0 = x0(5);
    c0 = x0(6);
    r_e = x0(4);
end

На входе подаю изображение I, по которому ищу окружность:

В массиве Lights записаны точки, отобранные предыдущими алгоритмами устранения выборосов:

Также на входе у меня начальный вектор значений x0, который я получаю после применения метода Таубина (вот тут я брала реализацию: https://uk.mathworks.com/matlabcentral/fileexchange/22678-circle-fit-taubin-method):

re = results_t(3);
r0 = results_t(1);
c0 = results_t(2);
ymin = unique(min(I(I~=0)));
ymax = unique(max(I(I~=0)));
k = 3;
x0 = [ymin,ymax,k,re,r0,c0];

results - это координаты пикселей с яркостью, не равной 0.

Моя проблема состоит в том, что dXk с каждой итерацией растет, поэтому я получаю недостаточно точные значения. Прошу помочь решить вопрос, как добиться сходимости алгоритма.

БЛОГ НА HUSL

Нелинейный МНК для круговой сигмоидальной функции (CSF-LS)

Ответы (0 шт):