Кодирование FDS в MATLAB

Question

В контексте вариационной постановки оптического потока у меня есть следующие связанные эллиптические уравнения в частных производных, которые представляют собой уравнения Эйлера-Лагранжа минимизирующего функционала:

, где - неизвестное векторное поле.Вышеприведенные уравнения должны быть решены для u и v.

Я пытался закодировать вышеупомянутое в MATLAB, но я не могу получить схему конечных разностей правильно.

РЕДАКТИРОВАТЬ: Вот код MATLAB:

clear;

[file1,path1]=uigetfile('*.*','Select Image 1');
image1=fullfile(path1,file1);
I1=double(imread(image1));

[file2,path2]=uigetfile('*.*','Select Image 2');
image2=fullfile(path2,file2);
I2=double(imread(image2));

% Convert the images to grayscale
[m,n,p]=size(I1);
if p==3
    I1=rgb2gray(I1);
    I2=rgb2gray(I2);
end

% Compute Image derivatives
h1=0.25*[-1 1;-1 1];
I_x=imfilter(I1,h1,'conv','same') + imfilter(I2,h1,'conv','same');

h2=0.25*[-1 -1;1 1];
I_y=imfilter(I1,h2,'conv','same') + imfilter(I2,h2,'conv','same');

h3=-0.25*ones(2);h4=0.25*ones(2);
I_t=imfilter(I1,h4,'conv','same') + imfilter(I2,h3,'conv','same');

% Compute Mixed Partial Derivatives (spatial and temporal)
[I_xx, I_xy] = gradient(I_x);
[I_yx, I_yy] = gradient(I_y);
[I_tx, I_ty] = gradient(I_t);

% Choose regularization parameter
alpha=10;

% Choose the kernel mask
ker=[1/12,1/6,1/12;1/6,0,1/6;1/12,1/6,1/12];

% Set the initial plot values
u=0;v=0;I=I1;
u_x=0;u_y=0;u_xx=0;u_yy=0;u_xy=0;u_yx=0;
v_x=0;v_y=0;v_xx=0;v_yy=0;v_xy=0;v_yx=0;

% Number of iterations
n_iter=50;
for i=1:n_iter

    if i>1
        % Compute derivatives of u
        [u_x,u_y]=gradient(u);
        [u_xx,u_xy]=gradient(u_x);
        [u_yx,u_yy]=gradient(u_y);

        % Compute derivatives of v
        [v_x,v_y]=gradient(v);
        [v_xx,v_xy]=gradient(v_x);
        [v_yx,v_yy]=gradient(v_y);
    end

    % Standard Gaussian Filter (HS)
    u_avg=imfilter(u,ker,'conv','same');
    v_avg=imfilter(v,ker,'conv','same');

    % Compute the known terms
    c1=I.*(I_tx + I_xx.*u + I_yx.*v + 2*I_x.*u_x...
      + I_x.*v_y + I_y.*v_x + I.*(u_xx+v_yx));

    c2=I.*(I_ty + I_yy.*v + I_xy.*u + 2*I_y.*v_y...
      + I_x.*u_y + I_y.*u_x + I.*(u_xy+v_yy));

    % Compute the RHS
    R1=alpha*u_avg-c1;
    R2=alpha*v_avg-c2;

    % Note the coefficients
    a=alpha+I.*I_x;
    b=I.*I_xy;
    c=alpha+I.*I_y;

    % determinant
    det=a.*c-b.^2;

    % Solve for u and v
    u=(c.*R1-b.*R2)./det;
    v=(a.*R2-b.*R1)./det;
end

d=figure;
subplot(221)
imshow(I1,[])
g1=get(gca,'Position');
title('Image 1')
subplot(222)
imshow(I2,[])
g2=get(gca,'Position');
title('Image 2')
avg=(g1(1)+g2(1))/2;
pos=[avg+0.02,g1(1),g2(3)-0.08,g2(4)];
subplot('Position',pos)
axis([xlim(1),xlim(2),ylim(1),ylim(2)])
quiver(u,v,3,'color','b','linewidth',1)
set(gca,'YDir','reverse')
title('Computed Optical Flow')

Чтобы кратко объяснить, что я сделал в коде.1. Я расширил члены LHS в связанных уравнениях.2. Использовал дискретизацию лапласиана (u_avg - u). 3. Скинул все смешанные производные и перекрестные члены на RHS и сохранил u и v на LHS.Полученные уравнения представляют собой систему уравнений вида a u + b v = R1 & b u + c v = R2, где члены a, b, c, R1, R2как в коде.Затем я попытался решить эту систему алгебраических уравнений с помощью уравнений Гаусса-Якоби с начальными итерациями (u, v) = (0,0).Код работает для синтетических изображений, но не для каких-либо реальных наборов данных изображений.