% Remove the "photobomb" hippo from the image
%
% Samuli Siltanen Dec 2015

% Read in the image
im = imread('Photobomb.jpg','jpg');
im = im(1:1320,:,:);
imR = MyDS2(double(im(:,:,1)));
imG = MyDS2(double(im(:,:,2)));
imB = MyDS2(double(im(:,:,3)));

% Remove the hippo from the red channel
hippoR    = imR(80:240,310:570);
[tmp1,tmp2] = size(hippoR);
row = tmp1-2;
col = tmp2-2;

% Construct the system matrix
disp('Constructing system matrix')
A = FD_Laplace(row,col);
disp('System matrix ready')

% Construct the Dirichlet data
vec_t = hippoR(1,2:(end-1));
vec_b = hippoR(end,2:(end-1));
vec_l = hippoR(2:(end-1),1);
vec_r = hippoR(2:(end-1),end);



% Initialize the right-hand side vector
m = zeros(row*col,1);

% Construct the right-hand side vector
% by looping over elements
disp('Constructing rhs')
for iii = 1:row
    for jjj = 1:col
        ind = (jjj-1)*row+iii;
        if jjj==1
            m(ind) = m(ind)+vec_l(iii);
        end
        if jjj==col
            m(ind) = m(ind)+vec_r(iii);
        end
        if iii==1
            m(ind) = m(ind)+vec_t(jjj);
        end
        if iii==row
            m(ind) = m(ind)+vec_b(jjj);
        end
    end
end
disp('Rhs constructed')

% Solve the Poisson equation
disp('Solving')
x = gmres(A,m,30,1e-6,100);
x = reshape(x,row,col);

% Substitute the result to the 
% hippo image
hippoR2 = hippoR;
hippoR2(2:(end-1),2:(end-1)) = x;

% Substitute the result to the whole image
imR2 = imR;
imR2(80:240,310:570) = hippoR2;

% Show the original and inpainted image patch
% using the same colorbar
MIN = min(min(hippoR(:)),min(hippoR2(:)));
MAX = max(max(hippoR(:)),max(hippoR2(:)));
figure(1)
clf
subplot(1,2,1)
imagesc(hippoR,[MIN,MAX])
axis equal
axis off
colormap gray
subplot(1,2,2)
imagesc(hippoR2,[MIN,MAX])
axis equal
axis off
colormap gray


% Show the original and inpainted image (red channel)
% using the same colorbar
MIN = min(min(imR(:)),min(imR2(:)));
MAX = max(max(imR(:)),max(imR2(:)));
figure(2)
clf
subplot(1,2,1)
imagesc(imR,[MIN,MAX])
axis equal
axis off
colormap gray
subplot(1,2,2)
imagesc(imR2,[MIN,MAX])
axis equal
axis off
colormap gray

% Save the images to file
im8bitR  = imR;
im8bitR  = im8bitR-MIN;
im8bitR  = im8bitR/(MAX-MIN);
im8bitR  = uint8(255*im8bitR);
im8bitR2 = imR2;
im8bitR2 = im8bitR2-MIN;
im8bitR2 = im8bitR2/(MAX-MIN);
im8bitR2 = uint8(255*im8bitR2);
imwrite(im8bitR,'HippoR.png','png')
imwrite(im8bitR2,'HippoR2.png','png')



break

% Show the image
figure(3)
clf
subplot(2,2,1)
imshow(im)
axis equal
axis off
subplot(2,2,2)
imshow(uint8(imR))
axis equal
axis off
subplot(2,2,3)
imshow(uint8(imG))
axis equal
axis off
subplot(2,2,4)
imshow(uint8(imB))
axis equal
axis off