structure_inpaint.m

this is a very very very nice code

%Inpaints the inpute image according to the mask, which is a binary image
%where 0 represents the regions to be inpainted. Parameters can be defined
%by param, or pass param=0 to use default values. An example of param can
%be seen below. T is the number of iterations, num_fill is the number of
%inpainting iterations to perform each T, and num_anis is the number of
%anisotropic diffusions to perform. num_initial is the number of initial
%"bleeding" iterations performed before inpainting to get a rough initial
%result. Delta is the update rate, and you typically want this value to be small
function out = structure_inpaint(in, mask, param)
if sum(class(param) == 'struct') < 6
    param = struct('T', 1000, 'num_fill', 15, 'num_anis', 2, 'num_initial', 300, 'delta', .1);
end
in = double(in);
out = in;

%We want the laplacian mask to be slightly larger than the input mask, in
%order for anisotropic diffusion to work
SE = strel('arbitrary',ones(3,3));
LapMask = imerode(mask,SE);

%Collect all the pixels to be inpainted and/or diffused
diff_set = get_pix_set(LapMask);
pix_set = get_pix_set(mask);

%Perform some initial diffusion. This just fills in the region to be
%inpainted with some color from the surrounding region. It serves as a good
%"first guess" to the inpainting result
for K = 1:param.num_initial
    out = out + param.delta.*comp_laplacian(out,diff_set,1);
end

for i = 1:param.T
    if mod(i,50) == 0
        disp([i, param.T])   %Update the user with the progress
    end
    
    %This is because normally you want intensities to be in the range [0,1]
    %but I choose to work in the range [0,255]. Thus I need to divide the
    %final result by 255 to keep it in the correct range.
    inpaint_delta = param.delta / 255;

    %Fill the region num_fill times
    for j = 1:param.num_fill
        %Compute gradients and the laplacian
        [Iy, Ix] = comp_backgrad(out,pix_set,2);
        [Iyb, Ixb] = comp_backgrad(out,pix_set,1);
        L = comp_laplacian(out,diff_set,1);
        nrm = get_norm(Iy,Ix,pix_set);
        
        %For every pixel to be inpainted...
        for k = 1:size(pix_set,1)
            y = pix_set(k,1);
            x = pix_set(k,2);
            if nrm(y,x) > 0
                %All the equations from the image inpainting paper
                dL = [L(y-1,x) - L(y+1,x), L(y,x+1) - L(y,x-1)];
                N = [Ix(y,x), -Iy(y,x)];
                beta = dL * N' / nrm(y,x);

                if beta > 0
                    Ixbm = min(0, Ixb(y,x));
                    IxfM = max(0, Ixb(y,x+1));
                    Iybm = min(0, Iyb(y,x));
                    IyfM = max(0, Iyb(y-1,x));
                    magI = sqrt(Ixbm^2 + IxfM^2 + Iybm^2 + IyfM^2);
                else
                    IxbM = max(0, Ixb(y,x));
                    Ixfm = min(0, Ixb(y,x+1));
                    IybM = max(0, Iyb(y,x));
                    Iyfm = min(0, Iyb(y-1,x));
                    magI = sqrt(IxbM^2 + Ixfm^2 + IybM^2 + Iyfm^2);
                end
                
                out(y,x) = out(y,x) + inpaint_delta * beta * magI;
            end
        end
    end
    
    %Do some anisotropic diffusion, the method used can be found in the report
    for j = 1:param.num_anis
        [Iy, Ix] = comp_backgrad(out,diff_set,2);
        Iyy = comp_laplacian(out,diff_set,2);
        Ixx = comp_laplacian(out,diff_set,3);
        Ixy = comp_laplacian(out,diff_set,4);
        Ix2 = Ix.^2;
        Iy2 = Iy.^2;
        norm2 = Ix2 + Iy2;
        for K = 1:size(diff_set,1)
            y = diff_set(K,1);
            x = diff_set(K,2);
            if norm2(y,x) > 0
                diff = (Ixx(y,x)*Iy2(y,x) - 2*Ixy(y,x)*Ix(y,x)*Iy(y,x) + Iyy(y,x)*Ix2(y,x)) / norm2(y,x);
                out(y,x) = out(y,x) + param.delta * diff;
            end
        end
    end
    
end