function output = multi_scale(input_c, example_c, IO, config)

addpath('../../../timeLapseCode/TimeLapseVideo/helpers/');
addpath('../AlignFace/');
addpath('tool/');

  nLevel = config.nLevel;
  color_space = config.color_space;

  % Multiscale warping 
  [vx_pyramid, vy_pyramid]=correspondence_pyramid(input_c, example_c, config);
  [input example range]=rgb_2_X(input_c, example_c, color_space);

  [h w nc] = size(input);

  % Multiscale matching
  output = input;

  for k = nc : -1 : 1

    if k == 1 
      matchLevel=[1:nLevel];
    else
      matchLevel = [1:nLevel];
    end

    radius_sizes = ceil(config.radius_scale*2.^(matchLevel-1));

    %input_edge = range*rgb2gray(input);
    %example_edge = range*rgb2gray(example);
    input_edge = input(:,:,1);
    example_edge = example(:,:,1);
    [output(:,:,k) input_pyramid output_pyramid example_pyramid gain_pyramid]...
      = multiscale_eq(input(:,:,k), example(:,:,k), nLevel, matchLevel, ...
                      radius_sizes, input_edge, example_edge, range,...
                      vx_pyramid, vy_pyramid, config);
  end

  output =  X_2_rgb(output, color_space) ;

  %% Post processing
  if config.transfer_eyes
    output = transfer_eyes(output, example_c);
  end
  if config.transfer_bg
    output = transfer_background(input_c, output, example_c);
  end

  if config.dump_debug
    debugger(input_pyramid, example_pyramid, output_pyramid, gain_pyramid,...
          IO.input_name, IO.style, range, matchLevel, 2*radius_sizes+1,...
          example_c, vx_pyramid, vy_pyramid);
  end


%% Dense correspondence 
%% ----------------------------------------------------------------------------
function [vx_pyramid, vy_pyramid]=...
                correspondence_pyramid(input, example, config)
addpath('../../../timeLapseCode/TimeLapseVideo/helpers/image_pyramids');
addpath('../../../timeLapseCode/TimeLapseVideo/helpers/matting');

% Downscale for matching 
scale = 0.25; 
[h,w,c]= size(input);
input = imresize(input, scale);
example = imresize(example, scale);

if config.dense_warp
  fprintf('Compute sift matching...\n');
  [vx, vy] = sift_matching(uint8(255*input), uint8(255*example));
  vx = imresize(vx, [h w])/scale;
  vy = imresize(vy, [h w])/scale;

  % Fill zeros for out-of-boundary
  [vx vy] = thresh_v(vx,vy);
  
  if config.user_warp

    addpath('../../../timeLapseCode/TimeLapseVideo/helpers/UserWarp/');
    [uvx uvy] = warp();  
    vx = vx + uvx;
    vy = vy + uvy;

  end

    

else
  vx = zeros(h,w);
  vy = zeros(h,w);
end

% Multi-scale dense warping
vx_pyramid = gaussian_pyramid(vx, config.nLevel, false);
vy_pyramid = gaussian_pyramid(vy, config.nLevel, false);
for i = 1 : config.nLevel
  [vx_pyramid{i}, vy_pyramid{i}] = thresh_v(vx_pyramid{i}, vy_pyramid{i});
end



function [vx, vy] = thresh_v(vx, vy)
[h w]=size(vx);
[~, mask] = warpImage(zeros(h,w,3), vx, vy);
vx(logical(1-mask)) = 0;
vy(logical(1-mask)) = 0;



function [vx, vy] = sift_matching(input, example)
sift_path='../../../timeLapseCode/TimeLapseVideo/helpers/SIFTflow/';
addpath(sift_path);
addpath(fullfile(sift_path,'mexDenseSIFT'));
addpath(fullfile(sift_path,'mexDiscreteFlow'));


cellsize= 7 ;
gridspacing=1;

tic;
sift1 = mexDenseSIFT(input,cellsize,gridspacing);
sift2 = mexDenseSIFT(example,cellsize,gridspacing);
toc

SIFTflowpara.alpha=500;
SIFTflowpara.gamma=10;
SIFTflowpara.d=1000000;
SIFTflowpara.nlevels=4;
SIFTflowpara.wsize=3;
SIFTflowpara.topwsize=4;
SIFTflowpara.nTopIterations=60;
SIFTflowpara.nIterations=30;


tic;[vx,vy,energylist]=SIFTflowc2f(sift1,sift2, SIFTflowpara);toc



%% Inpainting
%% --------------------------------------------------------------------
function output = parser_mask(input, input_name, ori_input)

ex_hair = double(imread('masks/gc/hair.png'))/255;
ex_skin = double(imread('masks/gc/skin.png'))/255;
ex_bg   = (1-ex_hair).*(1-ex_skin);

input_hair = double(imread(sprintf('masks/%s/hair.png', input_name)))/255;
input_skin = double(imread(sprintf('masks/%s/skin.png', input_name)))/255;
input_bg   = (1-input_hair).*(1-input_skin);

valid_mask = ex_hair.*input_hair + ex_bg.*input_bg + ex_skin.*input_skin;
v_mask_3 = cat(3, valid_mask, valid_mask, valid_mask);

output = inpainting( ori_input, input, 1-v_mask_3);

function output = inpainting(fg, bg, mask)

addpath('../../../timeLapseCode/TimeLapseVideo/helpers/poisson/');

[fh fv]= imgrad(fg);
[bh bv]= imgrad(bg);
ch = (1-mask).*bh + mask.*fh;
cv = (1-mask).*bv + mask.*fv;
cc = (1-mask).*bg + mask.*fg;
output = PoissonJacobi(cc, ch, cv, mask);

%% Transfer background 
%% ------------------------------------------------------------------------
function output = transfer_background(input, output, example)
addpath('../../../timeLapseCode/TimeLapseVideo/helpers/mexopencv/mexopencv/');
addpath('../../../timeLapseCode/TimeLapseVideo/helpers/matting/');
fprintf('Transferring background...\n');

import cv.*;
  alpha = extract_fg(input);
  bg    = reconstruct_bg(example);
  output = output.*alpha + (1-alpha).*bg;
  


%% Transfering eye highlights
%% --------------------------------------------------------------------------
function output = transfer_eyes(input, example)
addpath('../../../timeLapseCode/TimeLapseVideo/helpers/iris');

fprintf('Transferring highlight...\n');
% TODO: load from face template
output = transfer_eye(input, example, [240 490 220 170]);
output = transfer_eye(output, example, [540 490 220 170]);

function output = transfer_eye(input, example, roi)
  eye1 = getROI(example, roi);
  eye2 = getROI(input, roi);

  new_eye = iris_transfer(eye1, eye2);
  output = setROI(input, roi, new_eye);


function [output input_pyramid output_pyramid example_pyramid gain_pyramid]...
          = multiscale_eq(input, example, nLevel, matchLevel, ...
                          radius_sizes, input_edge, example_edge, range, ...
                          vx_pyramid, vy_pyramid, config)

  if nargin < 3 
    nLevel = 6 ; 
  end
  if nargin < 4
    matchLevel = [1:nLevel];
  end
  if nargin < 5
    radius_sizes = 3*2.^(matchLevel-1); 
  end

  addpath('../../../timeLapseCode/TimeLapseVideo/helpers/image_pyramids');
  addpath('../../../timeLapseCode/TimeLapseVideo/helpers/HistEqualization');

  if config.bilateral_decomp
    input_pyramid = bilateral_pyramid(input,nLevel, input_edge, range); 
    example_pyramid = bilateral_pyramid(example,nLevel, example_edge, range); 
  else
    input_pyramid = laplacian_pyramid(input,nLevel, false); 
    example_pyramid = laplacian_pyramid(example,nLevel, false); 
  end


  output_pyramid = input_pyramid;
  gain_pyramid   = input_pyramid; 

  fprintf('Multiscale matching...\n');
  for i = 1 : length(matchLevel)
    %output_pyramid{i} = HistTransfer(input_pyramid{i}, example_pyramid{i});
    level = matchLevel(i);
    fprintf('Layer %d...', level);
    if level == nLevel
      [output_pyramid{level} gain_pyramid{level}] = ...
                            LocalMeanVarMatch(input_pyramid{level}, ...
                            example_pyramid{level}, radius_sizes(i), ...
                            vx_pyramid{level}, vy_pyramid{level});
                            %LocalVarMatch(input_pyramid{level}, ...
    else
      [output_pyramid{level} gain_pyramid{level}] = ...
                            LocalVarMatch(input_pyramid{level}, ...
                            example_pyramid{level}, radius_sizes(i), ...
                            vx_pyramid{level}, vy_pyramid{level});
                            %LocalVarMatch(input_pyramid{level}, ...
    end
  end

  output = sum_pyramid(output_pyramid);
  fprintf('Done.\n');


function output = HistTransfer(input, example)

  input_plus = input >= 0 ; 
  input_minus = input < 0 ;
  output = HistTransferOneD(abs(input), abs(example));
  output(input_minus) = -1*output(input_minus);


function [output gain_map]= LocalVarMatchTwoSigns(input,example, radius, vx, vy)
input_plus = input.*(input>0);
input_minus = input.*(input<0);
example_plus = example.*(example>0);
example_minus = example.*(example<0);

[output_plus gain_plus] = LocalVarMatch(input_plus, example_plus, radius, vx, vy);
[output_minus gain_minus] = LocalVarMatch(input_minus, example_minus, radius, vx, vy);
%output = output_plus + output_minus; 
output = output_plus + output_minus; 
gain_map = gain_plus.*(input<0) + gain_plus.*(input>=0); 

function [output gain_map]= LocalVarMatch(input,example, radius, vx, vy)

  %epsilon = (max(abs(example(:)))/255)^2;
  epsilon = (max(input(:))/255)^2;
  %input_var =   epsilon*(2*radius+1)^2 +  boxfilter(input.^2, radius);
  %example_var =  boxfilter(example.^2, radius);

  h = fspecial('gaussian', ceil(6*radius*[1 1]), radius);
  epsilon = 0.001;
  input_var =   epsilon +  imfilter(input.^2, h);
  example_var =  imfilter(example.^2,h);

  warped_example_var = warpImage(example_var,vx,vy);

  gain_map = (warped_example_var./input_var).^0.5;
  
  output = input.*gain_map;
  %output = HistTransfer(output, example);

function [output gain_map]= LocalMeanVarMatch(input, example, radius, vx, vy)
  
  epsilon = (max(input(:))/255)^2;
  example_mean = boxfilter_n(example, radius);
  example_mean = warpImage(example_mean,vx,vy);
  input_contrast = input - boxfilter_n(input, radius);
  example_contrast = example - boxfilter_n(example, radius);

  input_var =   epsilon*(2*radius+1)^2 +  boxfilter(input_contrast.^2, radius);
  example_var =  boxfilter(example_contrast.^2, radius);

  example_var = warpImage(example_var,vx,vy);

  gain_map = (example_var./input_var).^0.5;

  %contrast_gain = (boxfilter(example_contrast.^2, radius)./boxfilter(input_contrast.^2, radius)).^0.5;
  output = example_mean + gain_map.*input_contrast;

function [output gain_map] = VarMatch(input, example, radius, vx, vy)
  [h w ] = size(input);
  output = input*0;
  gain_map = input*0;
  for i = 1 : 2*radius : h-2*radius
    for j = 1 : 2*radius : w-2*radius
      ii= i:i+2*radius-1;
      jj= j:j+2*radius-1;
      input_patch = input(ii,jj);
      ex_patch = example(ii,jj);

      [output_patch gain]= global_match(input_patch, ex_patch);

      gain_map(ii,jj)   =  gain; 
      output(ii,jj) = output_patch;
    end
  end
  

function [output gain] = global_match(input, example)
  mui     = mean(input(:)); 
  muo     = mean(example(:)); 
  gain   = std(example(:))/(std(input(:))+0.1);
  output = gain*(input - mui) + muo;

%% Debuging tool
%% ----------------------------------------------------------------------------
function debugger( input_pyramid, example_pyramid, output_pyramid, gain_pyramid,...
                   input_name, style,  max_, matchLevel, patch_sizes,...
                   example, vx_pyramid, vy_pyramid);

  if nargin < 6
    max_=1;
  end

    nLevel = length(input_pyramid);

    for i = 1 : nLevel
      sigma = 1*2^(i-1);

      if sum(matchLevel==i)>0
        patch_size = patch_sizes(find(matchLevel==i));
      else
        patch_size = 1;
      end

      if i == nLevel
        offset = 0;
        f=1;
      else
        offset = 0.5;
        f=2;
      end

      input_pyramid{i} = input_pyramid{i}/max_; 
      input_pyramid{i}(1:sigma, 1:sigma) = 1;
      input_pyramid{i}(end-patch_size-1:end, 1:patch_size) = 1;
      imwrite(f*input_pyramid{i}+offset, ...
      sprintf('results/%s_input_pyramid_%d_%s.jpg', input_name, i, style),... 
      'Quality', 100);

      example_pyramid{i} = example_pyramid{i}/max_; 
      example_pyramid{i}(1:sigma, 1:sigma) = 1;
      example_pyramid{i}(end-patch_size-1:end, 1:patch_size) = 1;
      imwrite(f*example_pyramid{i}+offset, ...
      sprintf('results/%s_example_pyramid_%d_%s.jpg', input_name, i, style),... 
      'Quality', 100);

      output_pyramid{i} = output_pyramid{i}/max_; 
      output_pyramid{i}(1:sigma, 1:sigma) = 1;
      output_pyramid{i}(end-patch_size-1:end, 1:patch_size) = 1;
      imwrite(f*output_pyramid{i}+offset, ...
      sprintf('results/%s_output_pyramid_%d_%s.jpg', input_name, i, style),... 
      'Quality', 100);

      imwrite(gain_pyramid{i}/5, ...
      sprintf('results/%s_gain_pyramid_%d_%s.jpg', input_name, i, style),... 
      'Quality', 100);

      imwrite(warpImage(example, vx_pyramid{i}, vy_pyramid{i}),...
      sprintf('results/%s_warp_pyramid_%d_%s.jpg', input_name, i, style),... 
      'Quality', 100);
    end

%% Utilities
%% --------------------------------------------------------------------
function output = setROI(input,roi, patch) 
  output = input;
  output(roi(2):roi(2)+roi(4)-1,... 
        roi(1):roi(1)+roi(3)-1, :) = patch;

function [input example range]=rgb_2_X(input_c, example, color_space)

  switch color_space
  case 'wb'
    input   = rgb2gray(input_c);
    example = rgb2gray(example);
    range = 1;
  case 'rgb'
    input = input_c;
    range = 1; 
  case 'lab'
    input   = RGB2Lab(input_c);
    example = RGB2Lab(example);
    range = 100;
  end

function output =  X_2_rgb(output, color_space) 

  switch color_space
  case 'lab'
    output = Lab2RGB(output);
  case 'wb'
    output = cat(3, output, output, output);
  otherwise
  end

