%% =========================================================================
% Demo codes for Texture-Synthesis using bi-directional similarity 
% Code Developed by Pulak Purkait : Date oct 10, 2012
% ----------------------------------------------------
% ECSU, Indian Statistical Institute Kolkata
% For any questions/bug related issues, send email to: pulak_r@isical.ac.in
% Referece 
%     P. Purkait and B chanda, Digital Restoration of Damaged Mural images,
%     ICVGIP 2012, Decembar 2012.
% =========================================================================
%%
% function [] = texture_synth3(image)
% image = 'cow';
% image = 'inpainting/bungee';
image = 'texture/25';
% image = 'inpainting/HillWeb';
% image = 'inpainting/Road_Trip';
% image = 'inpainting/man';

Aorg = imread([image '.jpg']);

close all;
% Aorg = repmat(Aorg, [1 1 3]);
cores = 4;    % Use more cores for more speed
algo = 'cputiled';
% algo = 'gpucpu';
tic

% nn_iters = 5; %defaut = 5
% rs_iters = 1; % default = 1;
verbose = true;

rs_factor2 = 2;
width = 23;

[m, n, rgb] = size(Aorg);
Borg = uint8(255*ones([rs_factor2*m, rs_factor2*n, rgb]));

[m2, n2, rgb] = size(Borg);

rs_factor = 1.5;
% Determinte starting scale
startscale = -ceil(log2(min(m2,n2))/log2(rs_factor))+rs_factor2+9;
scale = rs_factor^(startscale);

% Resize image to starting scale
A = imresize(Aorg,scale);
B = imresize(Borg,scale);
R = rand(size(B)) > 0.5;
Rnd = uint8(255*R);
B = Rnd;

[m, n, rgb] = size(A);
[m2, n2, rgb] = size(B);
sz = round([m2 - m, n2 - n]/2);
width2 = 2;%(width-1)/2-8;
iitt = 2;
iterations = max(ceil([m2 - m, n2 - n]/(2*width2)))+15;
B(sz(1)+1:sz(1)+m, sz(2)+1:sz(2)+n, :) = A;

%% Initial Approximation in lower scale
figure,
for iter = 1:iterations    
%     B(sz(1)+1:sz(1)+m, sz(2)+1:sz(2)+n, :) = A;
    if verbose
        fprintf('  Iteration %2d/%2d\n',iter,iterations);
        imshow(B);
        pause(0.001)
    end

    % Compute NN field
    m3 = min([sz(1)+m+iter*width2, m2]);
    n3 = min([sz(2)+n+iter*width2, n2]);
    
    m1 = max([sz(1)-iter*width2, 1]);
    n1 = max([sz(2)-iter*width2, 1]);
    
    for inner = 1:iitt
        PM_ann = nnmex(B(m1:m3, n1:n3, :), A, algo, width, [], [], [], [], [], cores);
        PM_bnn = nnmex(A, B(m1:m3, n1:n3, :), algo, width, [], [], [], [], [], cores, [], []);
        B(m1:m3, n1:n3, :) = votemex(A, PM_ann,  PM_bnn, algo, width, [], []);
%         B(m1:m3, n1:n3, :) = votemex(A, PM_ann,  [], algo, width, [], []);
    end
end

%% Refinement in Higher Scale
% iitt = 50;
for logscale = startscale+1:0

    scale = rs_factor^(logscale);
    %Upsample A for the next scale
    A = imresize(Aorg, scale);
    [m, n, rgb] = size(imresize(Borg, scale));
    B = imresize(B,[m n]);

    for iter = 1:iterations    
        if verbose
            fprintf('  Iteration %2d/%2d\n',iter,iterations);
            imshow(B);
            pause(0.001)
        end
        % Compute NN field
        PM_ann = nnmex(B, A, algo, width, [], [], [], [], [], cores, [], []);
        PM_bnn = nnmex(A, B, algo, width, [], [], [], [], [], cores, [], []);
        B = votemex(A, PM_ann,  PM_bnn, algo, width);
%         B = votemex(A, PM_ann,  [], algo, width, [], []);
    end

end
% PM_bnn = nnmex(A, B, algo, width, [], [], [], [], [], cores, [], []);
% B = votemex(A, PM_ann,  PM_bnn, algo, width);
toc

imshow(Aorg);
figure, imshow(B);

imwrite(B, [image '-texture.png']);
