mcv-m2-w1

sol_Laplace_Equation_Axb.m

function [u] = sol_Laplace_Equation_Axb(f, dom2Inp, param)
%this code is not intended to be efficient. 

[ni, nj] = size(f);

%We add the ghost boundaries (for the boundary conditions)
f_ext = zeros(ni+2, nj+2);
f_ext(2:end-1, 2:end-1) = f;
dom2Inp_ext = zeros(ni+2, nj+2);
dom2Inp_ext (2:end-1, 2:end-1) = dom2Inp;

%Store memory for the A matrix and the b vector    
nPixels = (ni+2)*(nj+2); %Number of pixels

%We will create A sparse, this is the number of nonzero positions
nnz_A = nPixels + 4*nnz(dom2Inp) + 2*(ni+2) + 2*(nj+2);

%idx_Ai: Vector for the nonZero i index of matrix A
%idx_Aj: Vector for the nonZero j index of matrix A
%a_ij: Vector for the value at position ij of matrix A
idx_Ai = zeros(nnz_A,1);
idx_Aj = zeros(nnz_A,1);
a_ij = zeros(nnz_A,1);

b = zeros(nPixels,1);

%Vector counter
idx = 1;

%North side boundary conditions
i = 1;
for j = 1:nj+2
    %from image matrix (i,j) coordinates to vectorial (p) coordinate
    p = (j-1)*(ni+2)+i;
    
    %Fill Idx_Ai, idx_Aj and a_ij with the corresponding values and
    %vector b
    idx_Ai(idx)= p;
    idx_Aj(idx) = p;
    a_ij(idx) = 1;
    idx = idx+1;
    
    idx_Ai(idx) = p;
    idx_Aj(idx) = p+1;
    a_ij(idx) = -1;
    idx = idx+1;

    b(p) = 0;
end

%South side boundary conditions
i = ni+2;
for j = 1:nj+2
    %from image matrix (i,j) coordinates to vectorial (p) coordinate
    p = (j-1)*(ni+2)+i;
    
    %Fill Idx_Ai, idx_Aj and a_ij with the corresponding values and
    %vector b
    %TO COMPLETE 2
    idx_Ai(idx)= p;
    idx_Aj(idx) = p;
    a_ij(idx) = 1;
    idx = idx+1;
    
    idx_Ai(idx) = p;
    idx_Aj(idx) = p-1;
    a_ij(idx) = -1;
    idx = idx+1;

    b(p) = 0;
end

%West side boundary conditions
j = 1;
for i = 1:ni+2
    %from image matrix (i,j) coordinates to vectorial (p) coordinate
    p = (j-1)*(ni+2)+i;

    %Fill Idx_Ai, idx_Aj and a_ij with the corresponding values and
    %vector b
    %TO COMPLETE 3
    idx_Ai(idx)= p;
    idx_Aj(idx) = p;
    a_ij(idx) = 1;
    idx = idx+1;
    
    idx_Ai(idx) = p;
    idx_Aj(idx) = p+(ni+2);
    a_ij(idx) = -1;
    idx = idx+1;

    b(p) = 0;
end

%East side boundary conditions
j = nj+2;
for i = 1:ni+2
    %from image matrix (i,j) coordinates to vectorial (p) coordinate
    p = (j-1)*(ni+2)+i;
    
    %Fill Idx_Ai, idx_Aj and a_ij with the corresponding values and
    %vector b
    %TO COMPLETE 4
    idx_Ai(idx)= p;
    idx_Aj(idx) = p;
    a_ij(idx) = 1;
    idx = idx+1;
    
    idx_Ai(idx) = p;
    idx_Aj(idx) = p-(ni+2);
    a_ij(idx) = -1;
    idx = idx+1;

    b(p) = 0;
end

%Inner points
for j = 2:nj+1
    for i = 2:ni+1
     
        %from image matrix (i,j) coordinates to vectorial (p) coordinate
        p = (j-1)*(ni+2)+i;
                                            
        if (dom2Inp_ext(i,j)==1) %If we have to inpaint this pixel
            
            %Fill Idx_Ai, idx_Aj and a_ij with the corresponding values and
            %vector b
            %TO COMPLETE 5
            idx_Ai(idx)= p;
            idx_Aj(idx) = p;
            a_ij(idx) = -4;
            idx = idx+1;

            idx_Ai(idx)= p;
            idx_Aj(idx) = p+1;
            a_ij(idx) = 1;
            idx = idx+1;

            idx_Ai(idx)= p;
            idx_Aj(idx) = p-1;
            a_ij(idx) = 1;
            idx=idx+1;

            idx_Ai(idx)= p;
            idx_Aj(idx) = p+(ni+2);
            a_ij(idx) = 1;
            idx = idx+1;

            idx_Ai(idx)= p;
            idx_Aj(idx) = p-(ni+2);
            a_ij(idx) = 1;
            idx = idx+1;

            b(p) = 0;
        else %we do not have to inpaint this pixel 
            
            %Fill Idx_Ai, idx_Aj and a_ij with the corresponding values and
            %vector b
            %TO COMPLETE 6
            idx_Ai(idx)= p;
            idx_Aj(idx) = p;
            a_ij(idx) = 1;
            idx = idx+1;

            b(p) = f_ext(i, j);
        end
    end
end
    %A is a sparse matrix, so for memory requirements we create a sparse
    %matrix
    %TO COMPLETE 7
    A = sparse(idx_Ai, idx_Aj, a_ij, nPixels, nPixels); %??? and ???? is the size of matrix A
    
    %Solve the sistem of equations
    x = mldivide(A,b);
    
    %From vector to matrix
    u_ext = reshape(x, ni+2, nj+2);
    
    %Eliminate the ghost boundaries
    u = full(u_ext(2:end-1, 2:end-1));
    

start.m

%Example script: You should replace the beginning of each function ('sol')
%with the name of your group. i.e. if your gropu name is 'G8' you should
%call :
% G8_DualTV_Inpainting_GD(I, mask, paramInp, paramROF)

clearvars;
%There are several black and white images to test:
%  image1_toRestore.jpg
%  image2_toRestore.jpg
%  image3_toRestore.jpg
%  image4_toRestore.jpg
%  image5_toRestore.jpg

%name = 'image5';
name = 'image1';

I = double(imread([ name '_toRestore.jpg']));
%I=I(1:10,1:10);

%Number of pixels for each dimension, and number of channels
[ni, nj, nC] = size(I);

if nC==3
    I = mean(I,3); %Convert to b/w. If you load a color image you should comment this line
end

%Normalize values into [0,1]
I = I-min(I(:));
I = I/max(I(:));

%Load the mask
mask_img = double(imread([name '_mask.jpg']));
%mask_img =mask_img(1:10,1:10);
[ni, nj, nC] = size(mask_img);
if nC==3
    mask_img = mask_img(:,:,1); %Convert to b/w. If you load a color image you should comment this line
end
%We want to inpaint those areas in which mask == 1
mask = mask_img>128; %mask(i,j) == 1 means we have lost information in that pixel
                     %mask(i,j) == 0 means we have information in that pixel
                                                                    
%%%Parameters for gradient descent (you do not need for week1)
param.dt = 5*10^-7;
param.iterMax = 10^4;
param.tol = 10^-5;

%%Parameters 
param.hi = 1 / (ni-1);
param.hj = 1 / (nj-1);

figure(1)
imshow(I);
title('Before')

Iinp = sol_Laplace_Equation_Axb(I, mask, param);

figure(2)
imshow(Iinp)
title('After');

%% Challenge image. (We have lost 99% of information)
clearvars
I = double(imread('image6_toRestore.tif'));
%Normalize values into [0,1]
I = I/256;

%Number of pixels for each dimension, and number of channels
[ni, nj, nC] = size(I);

mask_img = double(imread('image6_mask.tif'));
mask = mask_img>128; %mask(i,j) == 1 means we have lost information in that pixel
                     %mask(i,j) == 0 means we have information in that pixel

param.hi = 1 / (ni-1);
param.hj = 1 / (nj-1);

figure(1)
imshow(I);
title('Before')

Iinp(:,:,1) = sol_Laplace_Equation_Axb(I(:,:,1), mask(:,:,1), param);
Iinp(:,:,2) = sol_Laplace_Equation_Axb(I(:,:,2), mask(:,:,2), param);
Iinp(:,:,3) = sol_Laplace_Equation_Axb(I(:,:,3), mask(:,:,3), param);

figure(2)
imshow(Iinp)
title('After');

%% Goal Image
clearvars;

%Read the image
I = double(imread('Image_to_Restore.png'));

%Number of pixels for each dimension, and number of channels
[ni, nj, nC] = size(I);

%Normalize values into [0,1]
I = I - min(I(:));
I = I / max(I(:));

%We want to inpaint those areas in which mask == 1 (red part of the image)
I_ch1 = I(:,:,1);
I_ch2 = I(:,:,2);
I_ch3 = I(:,:,3);

%TO COMPLETE 1
mask = I_ch1==1 & I_ch2==0 & I_ch3==0; %mask(i,j) == 1 means we have lost information in that pixel
                                       %mask(i,j) == 0 means we have information in that pixel

%%%Parameters for gradient descent (you do not need for week1)
%param.dt = 5*10^-7;
%param.iterMax = 10^4;
%param.tol = 10^-5;

%parameters
param.hi = 1 / (ni-1);
param.hj = 1 / (nj-1);

% for each channel 

figure(1)
imshow(I);
title('Before')

Iinp(:,:,1) = sol_Laplace_Equation_Axb(I_ch1, mask, param);
Iinp(:,:,2) = sol_Laplace_Equation_Axb(I_ch2, mask, param);
Iinp(:,:,3) = sol_Laplace_Equation_Axb(I_ch3, mask, param);
    
figure(2)
imshow(Iinp)
title('After');