clear
close all;

bands = 5;
sigmaPre = 0;  % Stdev of gaussian prefilter (set to 0, 8, 16 etc.)

im0 = pgmRead('orange.pgm');
im1 = pgmRead('apple.pgm');
mask = pgmRead('mask.pgm')/255;

%% Build Laplacian pyramids
[pyr0,pind0] = buildLpyr(im0,bands); 
showLpyr(pyr0,pind0);
pause(1);

[pyr1,pind1] = buildLpyr(im1,bands); 
showLpyr(pyr1,pind1);
pause(1);


if sigmaPre > 0
  sigmaSqr = sigmaPre*sigmaPre;
  gFiltSize = 2 * round(3.0 * sigmaPre) + 1;
  x = [1:gFiltSize] - round((gFiltSize+1)/2);
  gFilt = exp(- x .* x / (2.0*sigmaSqr));
  gFilt = gFilt/ sum(gFilt(:));

  mask0 = rconv2sep(mask, gFilt, gFilt);

else
  mask0 = mask;
end

%% Build Gaussian mask pyramid
binom5 = binomialFilter(5);
lo_filt = binom5*binom5';  %construct a separable 2D filter

maskPyr = zeros(size(pyr1));
tmp = mask0;
maskPyr(1:prod(pind1(1,:))) = tmp;
for b=2:bands
 tmp = corrDn(tmp,lo_filt,'reflect1',[2 2]);
 ind0 = sum(prod(pind1(1:(b-1), :),2));
 maskPyr((ind0+1):(ind0+prod(pind1(b,:)))) = tmp(:);
end

for b=1:bands
 showIm(pyrBand(maskPyr,pind1,b));
 pause(1);
end


blendPyr = zeros(size(pyr1));
for b=1:bands
 ind0 = sum(prod(pind1(1:(b-1), :),2));
 subband0 = pyrBand(pyr0, pind0, b);
 subband1 =  pyrBand(pyr1, pind1, b);
 maskBand = pyrBand(maskPyr, pind1, b);

 tmp = maskBand .* subband0 + (1-maskBand) .* subband1;
 blendPyr((ind0+1):(ind0+prod(pind1(b,:)))) = tmp(:);
end

showLpyr(blendPyr,pind1);
pause(1);


blendIm = reconLpyr(blendPyr,pind1);
displayImage(blendIm);