HS и LK надежный оптический поток

Question

Я пытаюсь сделать исходный алгоритм LK и HS устойчивым, используя метод Black & Anandan Lorentzian

Я не знал, как применить эту концепцию в моем коде MATLAB!... и когда я пытался применить то, что я знаю, у меня были худшие результаты, чем у исходного HS

, я прикрепляю свой код MATLAB ... если кто-нибудь знает, как это сделать ... или гдеэто мои ошибки .. Я буду признателен за это

этот код считывает 4 серых изображения (i0, i1, i2, i3) фактического процесса устойчивости в execute_reg_iter (я пытался применить к данным и сглаживающему члену), но этодает мне худший поток, чем оригинальный метод

function [regularized_u,regularized_v]=regularization_robust (i0,i1,i2,i3)
  alpha=10.0;
  numpass=100; 
  ROBUST_FLAG = true;
  [pic_x, pic_y]=size(i0); 

  I=cell(1,2);
  I{1,1}=double(box_filter(i0,i1,i2));
  I{1,2}=double(box_filter(i1,i2,i3));

  %Simocelli 
  [Ix,Iy,It]=apply_Simoncelli_filters(I{1,1},I{1,2});
   [row, col]= size(Ix); 
   %%%caculate the flow %%%%%%
  [regularized_u,regularized_v]= calc_intensity_regularization_flow(Ix,Iy,It,alpha,numpass, ROBUST_FLAG); 

        %%%%% print output            

  end 


 function [seq]=box_filter(seq1,seq2,seq3)
   box(1)=0.25;
   box(2)=0.5;
   box(3)=0.25;

   seq=box(1)*seq1+box(2)*seq2+box(3)*seq3;
   temp=imfilter(seq,box,'same','conv','symmetric');
   seq=imfilter(temp,box','same','conv','symmetric');
end % box_filter

  function [Ix,Iy,It]=apply_Simoncelli_filters(I1,I2)
    It=I2-I1;
    gd=[1 -8 0 8 -1]/12;
    Ix=-imfilter(I1,gd,'conv','symmetric','same');
    Iy=-imfilter(I1,gd','conv','symmetric','same');
  end 

 % Robust function 
function [y]= robust_functions (type, x, sigma, order)

switch (type)

case 'lorentzian'

    switch (order)
        case 0
            y = log(1 + x.^2 / (2 * sigma^2));
        case 1
            y = 2 * x ./ (2 * sigma^2 + x.^2);
             if (2 * sigma^2 + x.^2) ==0
                y = 100; 
             end 
        case 2
            y = 2 ./ (2 * sigma^2 + x.^2);
    end


case 'quadratic'
    switch (order)
        case 0
            y = x.^2 / sigma^2;
        case 1
            y = 2 * x / sigma^2;
        case 2
            y = repmat(2 / sigma^2, size(x));
    end
otherwise
    error('Invalid robust function type.');

 end

 end

function [regularized_u,regularized_v]=...
        calc_intensity_regularization_flow(Ix,Iy,It,alpha,numpass,ROBUST_FLAG)

 regularized_u=zeros(size(Ix),'double');
 regularized_v=zeros(size(Iy),'double');

% Perform iterations 
for iter_no=0:2:numpass


 [u2,v2]=perform_reg_iter(regularized_u,regularized_v,...
Ix,Iy,It,alpha,ROBUST_FLAG);

[regularized_u,regularized_v]=perform_reg_iter(u2,v2,Ix,Iy,It,...
 alpha,ROBUST_FLAG);

end 
end

function [u,v]=perform_reg_iter(u,v,Ix,Iy,It,alpha,ROBUST_FLAG)

 kernel_1=[1/12 1/6 1/12;1/6 0 1/6;1/12 1/6 1/12];
 estimator_weight=0.0;
 robust_type= 'lorentzian'; % 'lorentzian'; % charbonnier
 % (0) original ... (1) first derivative ... (2) 2 ./ (2 * sigma^2 + x.^2)
  % in Sun paper 
 robust_order=2;
 uAvg=conv2(u,kernel_1,'same');
 vAvg=conv2(v,kernel_1,'same');
if (ROBUST_FLAG)
    % Smoothing term 

  %  robust_sigma =1;  %0.03 in Sun paper 
  %uAvg = feval('robust_functions',robust_type,uAvg,robust_sigma,robust_order);
  %vAvg = feval('robust_functions',robust_type,vAvg,robust_sigma,robust_order);

    %Data term 
    robust_sigma=0.01;  % 1.5 in Sun paper 
    value = Ix.*u + Iy.*v + It; 
    estimator_weight = 
    feval('robust_functions',robust_type,value,robust_sigma,robust_order);        
    Ix= Ix.* estimator_weight; 
    Iy= Iy.* estimator_weight; 
    It= It.* estimator_weight; 

end 


u= uAvg - ( Ix .* ( ( Ix .* uAvg ) + ( Iy .* vAvg ) + It ) ) ./ ...
 ( alpha^2 + Ix.^2 + Iy.^2); 
v= vAvg - ( Iy .* ( ( Ix .* uAvg ) + ( Iy .* vAvg ) + It ) ) ./ ...
 ( alpha^2 + Ix.^2 + Iy.^2);
  end

Без надежного расчета я получил этот результат

с надежным Я получил этот результат