双边滤波

function I2 = BF_HO1(I,hin,win,Imax,Bmax,Wq,sigma_r)

% Initional
I = double(I) ;
hend=size(I,1);
wend=size(I,2);
I1 = double(I) ; I2 = I1;

% Added Edge
II = ones(hend+2*hin,wend+2*win) ;
II(hin+1:hend+hin,win+1:wend+win) = I1;
I1 = II ; I2 = I1 ;

% the number of bins
b = 0:Bmax-1 ;
b = b’ ;
b = b.*(Imax/Bmax) ;

HisI0 = zeros(Bmax,wend);

%%%—the first row—%%%
for k1=hin+1:1:hin+1
for k2=win+1:1:wend+win
%%%— speed 2—%%%
if k2 == win+1
I0 = I1(k1-hin:k1+hin,k2-win:k2+win);
HisI0(:,k2) = myhist(I0,Wq,Imax,Bmax);
HisI = HisI0(:,k2) ;
else
subI0 = I1(k1-hin:k1+hin,k2-win-1);
HissubI0 = myhist(subI0,Wq(1:2*hin+1,1),Imax,Bmax);
%%%
addI0 = I1(k1-hin:k1+hin,k2+win);
HisaddI0 = myhist(addI0,Wq(1:2*hin+1,2*win+1),Imax,Bmax);

       %%%
       HisI0(:,k2) = HisI0(:,k2-1) + HisaddI0' - HissubI0' ;  
       HisI = HisI0(:,k2) ;             
    end       

    % BF
     Rfg1 =  HisI ;  
     % Rfg2 = (1/(sigma_r*sqrt(2*pi)))*exp(-(b-I1(k1,k2)).^2/(2*sigma_r^2));
     Rfg2 = exp(-(b-I1(k1,k2)).^2/(2*sigma_r^2));
     Rfg12 =  Rfg1.*Rfg2 ; 
     Rfg22 =  Rfg12.*(b) ;
     Rfg3 = sum(sum(Rfg22(:))) ;
     I2(k1,k2) = Rfg3/(sum(Rfg12(:)) + eps) ;         
    % BF               

end
end

%%%—other rows—%%%
for k1=hin+2:1:hend+hin
for k2=win+1:1:wend+win

%%%--- no speed---%%%
  % I0 = I1(k1-hin:k1+hin,k2-win:k2+win);          
  % HisI0 = myhist(I0,Wq,Imax,Bmax);        

% %%%— speed 1—%%%
% if k2==win+1
% I0=I1(k1-hin:k1+hin,k2-win:k2+win);
% HisI0 = myhist(I0,Wq,Imax,Bmax);
% else
% subI0 = I1(k1-hin:k1+hin,k2-win-1);
% HissubI0 = myhist(subI0,Wq,Imax,Bmax);
% addI0 = I1(k1-hin:k1+hin,k2+win);
% HisaddI0 = myhist(addI0,Wq,Imax,Bmax);
% HisI0 = HisI0 + HisaddI0 - HissubI0 ;
% end

% %%%— speed 2—%%%
% subI00 = I1(k1-hin-1,k2-win:k2+win);
% HissubI00 = myhist(subI00,Wq(1,1:2*win+1),Imax,Bmax);
% %%%
% addI00 = I1(k1+hin,k2-win:k2+win);
% HisaddI00 = myhist(addI00,Wq(2*hin+1,1:2*win+1),Imax,Bmax);
% %%%
% HisI0(:,k2) = HisI0(:,k2) + HisaddI00’ - HissubI00’ ;
% HisI = HisI0(:,k2) ;

   %%%--- speed 3---%%%         
       subI00 = I1(k1-hin-1,k2-win:k2+win);
       HissubI00 = myhist(subI00,Wq(1,1:2*win+1),Imax,Bmax);           
   if k2==win+1
       %%%
       addI00 = I1(k1+hin,k2-win:k2+win);
       HisaddI00 = myhist(addI00,Wq(2*hin+1,1:2*win+1),Imax,Bmax);
       %%%            
   else
       %%%
       HisaddI000 = myhist(I1(k1+hin,k2+win),Wq(2*hin+1,2*win+1),Imax,Bmax);
       HissubI000 = myhist(I1(k1+hin,k2-win-1),Wq(2*hin+1,1),Imax,Bmax);           
       HisaddI00 = HisaddI00 + HisaddI000 - HissubI000 ;
       %%%            
   end       
       HisI0(:,k2) = HisI0(:,k2) + HisaddI00' - HissubI00' ;
       HisI = HisI0(:,k2) ;

   %%%--- speed 4---%%% ？



    % BF
     Rfg1 =  HisI ;  
     % Rfg2 = (1/(sigma_r*sqrt(2*pi)))*exp(-(b-I1(k1,k2)).^2/(2*sigma_r^2));
     Rfg2 = exp(-(b-I1(k1,k2)).^2/(2*sigma_r^2));
     Rfg12 =  Rfg1.*Rfg2 ; 
     Rfg22 =  Rfg12.*(b) ;
     Rfg3 = sum(sum(Rfg22(:))) ;
     I2(k1,k2) = Rfg3/(sum(Rfg12(:)) + eps) ;         
    % BF         

end
end
% Comeback
I1 = I2(hin+1:hend+hin,win+1:wend+win);
I2 = I1 ;

%%%%%%%%%%%%%%— function end —%%%%%%%%%%%%%
end