%% This is the handout file for CSC487/2503, A2, Question 1, Fall 2008 

clear 
%close all

% Dependencies: 
%   Matlab iseToolbox and utvisToolbox
%   In particular, you will need the pyramid tools available
%   in iseToolbox/pyrTools.

%   You need to set the path to include these toolboxes.
%   This is done by running the startup.m file (see below).
if isunix
  % You need to change the following two directories
  % if you are running under unix or linux.
  if ~exist('matlabVisRoot', 'var')
    run('/u/jepson/pub/matlab/startup.m');
  end
  codeDir = './';
else
  % You need to change the following two directories and
  % make sure you have run the startup.m file that
  % came with the utvisToolbox.
  if ~exist('matlabVisRoot', 'var')
   run('C:/Matlab/startup.m');
  end
  codeDir = './';
end
dataDir = [codeDir, '../data/'];

% Parameters:
chatty = true;  % Pause and display intermediate results.
maxLevels = 7;     % Number of pyramid levels to use. You can change this.

%% Read in the apple, orange, and mask images
im1 = pgmRead([dataDir 'apple.pgm']);
im2 = pgmRead([dataDir 'orange.pgm']);
maskIm = pgmRead([dataDir 'mask.pgm']);
% Convert mask to a logical image
middle = (max(max(maskIm)) + min(min(maskIm)))/2;
mask = maskIm > middle;
clear maskIm

% Show the three images.
figure(1); clf;
showIm(im1);

figure(2); clf;
showIm(im2);

hardMerge = im1;
hardMerge(mask) = im2(mask);
figure(3); clf;
showIm(hardMerge);

if chatty
  fprintf(2, 'Press any key to continue ... ');
  pause; 
end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% You can change any code below here.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% I have provided some demo code to help you get started.

% Compute Laplacian pyramid
[pyr1, pind1] = buildLpyr(im1, maxLevels);

figure(3); clf; 
showLpyr(pyr1, pind1);

if chatty
  fprintf(2, 'Press any key to continue ... ');
  pause; 
end

% Note pyr1 is stored as a 1D vector
size(pyr1)

% The size of the pyramid images is in
pind1

% The pyramid is stored finest scale first, next finest,
% and so on. So to get only the second band, grab the following
% subvector out of pyr:
bnd = pyr1(prod(pind1(1,:)) + (1:prod(pind1(2,:))));
bnd = reshape(bnd, pind1(2,:));
figure(3); clf; showIm(bnd);

% You can get the same range of indices for a particular
% band more easily using pyrBandIndices.m.  The following
% max absolute difference should be zero:
max(abs(bnd(:) - pyr1(pyrBandIndices(pind1, 2))))

% You can also access this same band, and have it
% reshaped to the correct image shape, using pyrBand.m
max(max(abs(bnd - pyrBand(pyr1, pind1, 2))))

if chatty
  fprintf(2, 'Press any key to continue ... ');
  pause; 
end

% You will need to reset particular bands in the pyramid.
% For a bogus example, let me just amplify the second band
% by a factor of 4... this will increase the contrast of the 
% edges.
modPyr = pyr1;
modPyr(pyrBandIndices(pind1, 2)) = 4* bnd(:);

% You will need to rebuild the image from the Laplace pyramid
modRecon = reconLpyr(modPyr, pind1);

figure(3); clf; 
showIm(modRecon);

% If you undo the contast enhancement...
modPyr(pyrBandIndices(pind1, 2)) = bnd(:);

% ... and rebuild ...
modRecon = reconLpyr(modPyr, pind1);

% ... you get the original image back
max(max(abs(modRecon - im1)))

if chatty
  fprintf(2, 'Press any key to continue ... ');
  pause; 
end

% You also need to build a Gaussian pyramid of the mask image.
[maskPyr, maskPind] = buildgpyr(double(mask), maxLevels);

% Ughh.. note the lower case g.

% The pyramids are the identical sizes
max(abs(pind1 - maskPind))

% The last band is:
bnd = pyrBand(maskPyr, maskPind, maxLevels);

figure(3); clf; 
showIm(bnd);

% Note that this Gaussian pyramid band has actually been
% scaled up.  The maximum value of level lev is 2^(lev-1) not 1.
% Here the level is given by the variable maxLevels:
max(max(bnd))

% You will need to rescale the band down by a factor of 2^(1-lev)
% before you do the mixture between the two pyramids described in 
% the lecture notes.

if chatty
  fprintf(2, 'Press any key to continue ... ');
  pause; 
end

% To see what this lowest band of the Gaussian pyramid looks
% like when interpolated back up, we can use reconLpyr. 
% (Do you see why?)
blurPyr = zeros(size(pyr1));
blurPyr(pyrBandIndices(pind1, maxLevels)) = bnd(:);

blurMask = reconLpyr(blurPyr, pind1);
figure(3); clf;
showIm(blurMask);


