function k=FastInpaint(I,mask,liter)%(I,mask,liter,bar)
%Image inpainting by Manuel M. Oliveira's method Fast Digital Image
%Inpainting 
% 09-10-2007

%parameters:
%     I------image to be inpainted
%  mask------the noise mask which is a binary image
% liter------literation times
%   bar______diffusion barrier
tic %timer start
I=im2double(I);
     if islogical(mask)
       mask=mask;    
    else
        mask=im2bw(mask);
    end
diffker=[0.073235 0.176765 0.0732325; 0.17675 0 0.17675; 0.073235 0.176775 0.073235]
%There are two kinds of diffusion kernels in Olivera's article I use the
%first one[0.073235 0.176765 0.0732325; 0.17675 0 0.17675; 0.073235 0.176775 0.073235]. 
%Another is [0.125 0.125 0.125; 0.125 0 0.125; 0.125 0.125 0.125].
[r,c]=find(mask);
if (size(I, 3) == 3)  %Color image process
    r=FInpaintGray(I(:,:,1),mask,liter,diffker);
    g=FInpaintGray(I(:,:,2),mask,liter,diffker);
    b=FInpaintGray(I(:,:,3),mask,liter,diffker);
    k=cat(3,r,g,b);
else                  %Gray image process
    k=FInpaintGray(I,mask,liter,diffker);
end
toc %show the elapsed time.

%------------------------------------------------------------------------%
function g=FInpaintGray(I,mask,liter,diffker)
[r,c]=find(mask);
f=incrSize(I);
for n=1:liter
   for i=1:length(r)        
         x=r(i)+1;%+1
         y=c(i)+1;%+1
         f(x,y)=f(x-1,y-1)*0.073235 + f(x-1,y)*0.176765 + f(x-1,y+1)*0.073235 + f(x,y-1)*0.176765 +f(x,y+1)*0.176765 + f(x+1,y-1)*0.073235 + f(x+1,y)*0.176765 + f(x+1,y+1)*0.073235;
         
         %f(x,y)=f(x-1,y-1)*diffker(1,1)+f(x-1,y)*diffker(1,2)+f(x-1,y+1)*diffker(1,3)+f(x,y-1)*diffker(2,1)+f(x,y+1)*diffker(2,3)+f(x+1,y-1)*diffker(3,1)+f(x+1,y)*diffker(3,2)+f(x+1,y+1)*diffker(3,3);
    end
end
g=mat2gray(f(2:end-1,2:end-1));
%---------------------------------------------------------------------%
function u=incrSize(f)
A=im2double(zeros(size(f,1)+2,size(f,2)+2));
A(2:end-1,2:end-1)=f;
A(1:1,2:end-1)=f(1:1,1:end);
A(2:end-1,1:1)=f(1:end,1:1);
A(2:end-1,end:end)=f(1:end,end:end);
A(end:end,2:end-1)=f(end:end,1:end);
u=A;