function [ result, membrane ] = ConvolutionPyramid( source, target, mask )
%CONVOLUTIONPYRAMID Summary of this function goes here
%   Detailed explanation goes here

start = start_timer('Starting Convolution Pyramid... ');



% Copy pixels
if (false)
    
    % Choose upper-left corner and cut a part of the image
    sm = size(mask);
    
    
    temp = target;
    for i=1:3

        sr = source(:,:,i);
        cr = temp(:,:,i);
        cr(mask) = sr(mask);
        temp(:,:,i) = cr;

    end
    res_copy = temp;
    figure(); imshow(res_copy); title('Cut&Paste Pixels');
end


% Characteristic function: 1 on the boundary, 0 otherwise
h = fspecial('laplacian', 0);
chi = imfilter(double(mask), h);
chi(chi<0) = 0;
chi(chi>0) = 1;

erf = target - source;

% Filter set: h1, h2 and g
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w = [0.1507 0.6836 1.0334 0.0270 0.0312 0.7753];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h1 = w(1:3);
h1 = [h1 h1(end-1:-1:1)];
h1 = h1' * h1;

h2 = h1*w(4);

g = w(5:end);
g = [g g(end-1:-1:1)];
g = g' * g;

membrane = zeros(size(source));
result   = zeros(size(target));

for i=1:3
    
    tr = target(:,:,i);   
    
    a = erf(:,:,i);
    a(~chi) =  0;
    
    % Convolution Pyramid 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Ierf = evalf( a, h1, h2, g );
    Ichi = evalf( chi, h1, h2, g );

    membrane(:,:,i) = Ierf./Ichi;
    temp = membrane(:,:,i) + source(:,:,i);
    
    mem = (   abs(membrane(:,:,i)) ./ 2  )  .* mask;
    %im = imcrop(mem, [466 335 340 256]);
    X = 1:3:size(mem,1); 
    Y = 1:3:size(mem,2);
    
    %if (i==1) 
        %createfigure1(mem(X,Y));  
    %end;
    
    tr(mask) = temp(mask);
    result(:,:,i) = tr;
   
end

end_timer(start);

end


function t = start_timer(msg)
	fprintf(1, msg); t = tic;
end

function end_timer(timer)
	t = toc(timer); fprintf(1, '%.8f sec\n', t);
end

