% Third attempt at inpainting, now using the
% Poisson kernel. This is Laplace inpainting.
%
% Samuli Siltanen Nov 2014

% Parameters for plotting
lwidth = 3;

% Import the photo
im = imread('Photobomb.jpg','jpg');
im = double(im(:,:,2));
[row,col] = size(im);

% Crop the photo to square shape
im = im(:,end/5+[1:row]);
im = im(1:(end-1),1:(end-1));
im = MyDS2(im);
im = MyDS2(im);
[row,col] = size(im);

% Normalize the image
im = im-min(im(:));
im = im/max(im(:));

% Create coordinates for the image
D = 1/col;
t = D/2 + [0:(col-1)]*D;
[X,Y] = meshgrid(t);

% Construct angular variable for plotting the circle
% and for representing pixel values at the boundary
Nfii = 128;
fii  = [0:(Nfii-1)]/Nfii*2*pi;
fii  = fii-pi;
Dfii = fii(2)-fii(1);

% Choose a circle exactly surrounding the inpainting area,
% here in Matlab's image plotting coordinates
% determined by command "imagesc"
cenx = 120;
ceny = 80;
R    = 62;
% here the same in unit square coordinates
% 0<=x<=1 and 0<=y<=1:
cenx01 = cenx/col;
ceny01 = ceny/row;
R01    = R/row;
R_unit = R01+0.05;

% Evaluate pixel values along a circle slightly outside of the
% inpainting domain
f = interp2(X,Y,im,cenx01+R_unit*cos(fii),ceny01+R_unit*sin(fii));

% Perform Laplace inpainting. Loop over pixels located in the inpainting
% domain
im2          = im;
inpaint_ind  = sqrt((X-cenx01).^2+(Y-ceny01).^2)<R01;
xvec         = X(inpaint_ind);
yvec         = Y(inpaint_ind);
tmpvec       = zeros(size(im2(inpaint_ind)));
for iii = 1:length(xvec)
    x     = xvec(iii);
    y     = yvec(iii);
    z     = x-cenx01 + 1i*(y-ceny01);
    r     = abs(z)/R_unit;
    theta = angle(z);
    poisker     = (1-r^2)./(1-2*r*cos(theta-fii)+r^2);
    tmpvec(iii) = 1/(2*pi)*real(Dfii*sum(f(:).*poisker(:)));
    disp(tmpvec(iii))
    % Monitor the run
    if mod(iii,100)==0
       disp(['Done ', num2str(iii), ' out of ', num2str(length(xvec))]) 
    end
end
im2(inpaint_ind) = tmpvec;

% Take a look at the image and the two windows
figure(1)
clf
% Take a look at the image and the two windows
figure(1)
clf
subplot(1,2,1)
imagesc(im)
hold on
plot(cenx+R*cos(fii),ceny+R*sin(fii),'r','linewidth',lwidth)
colormap gray
axis equal
axis off
subplot(1,2,2)
imagesc(im2)
%hold on
%plot(cenx+R*cos(fii),ceny+R*sin(fii),'r','linewidth',lwidth)
colormap gray
axis equal
axis off


