function demoCodedAperturePair;
% This code is to demonstrate: defocus, depth estimation, defocus deblurring
% This code also compares the performance of our optimized coded aperture with the circular aperture
% pair

% written by Changyin Zhou @ Columbia U, September 2009

inFolder = 'input\';
outFolder = 'output\';

levelOfGT = 20;
levelOfET = 20;
levelOfStep = 0.5;

inSigma = 0.003; %Noise level
pattern = 'c';  %Choose the aperture pattern pair

%% Load Image and Depth map
f0 = im2double(imread([inFolder, 'texture.png']));
[hei, wid, ch] = size(f0);

[X, Y] = meshgrid(1:wid, 1:hei);
blurSizeMap = Y/max(Y(:))*6 + 1;

%Discretized the blursize map
level = levelOfGT;
blurSizeMap = floor(blurSizeMap/max(abs(blurSizeMap(:))) * level)* (max(abs(blurSizeMap(:)))/level); %40 levels
blurSizeMap(blurSizeMap<=1&blurSizeMap>=0) = 1; % |blursize| <= 1 just means in focus
blurSizeMap(blurSizeMap>=-1&blurSizeMap<0) = -1;

figure; imshow(f0); title('Focused Image: Rich Texture on left, Weak Texture on right');
figure; imagesc(blurSizeMap); title('Depth Map');

%% Load Aperture Pair
sz = 33; %Can be any number 
[fullKernel1, fullKernel2] = zKernel('c', sz, sz);
transmission = sum(fullKernel1(:));
[inKernel1, inKernel2] = zKernel(pattern, sz, sz);
inKernel1 = inKernel1/transmission;
inKernel2 = inKernel2/transmission;

figure; imagesc([inKernel1, zeros(sz, 10), inKernel2]); title('Coded Aperture Pair');

%% Synthesize defocused image
for ch = 1:3
    f1(:, :, ch) = Defocus(f0(:, :, ch), blurSizeMap, inKernel1, inSigma);
    f2(:, :, ch) = Defocus(f0(:, :, ch), blurSizeMap, inKernel2, inSigma);
end;

figure; imshow([f1, zeros(hei, 10, 3), f2]); title('Two Defocused Images');

imwrite(f1, [outFolder, 'f1.png']);
imwrite(f2, [outFolder, 'f2.png']);

%% Depth Estimation
inRange = [min(blurSizeMap(:)) - 0.5:levelOfStep:max(blurSizeMap(:))+ 0.5];
depth_map = zDFD(f1, f2, inKernel1, inKernel2, inSigma, inRange, 0);

figure; imagesc(depth_map); color bar;
title('Estimated depth map');
save 'estimatedDepth.mat' 'depth_map'
save 'groundtruthDepth.mat' 'blurSizeMap'

%% Defocus Deblurring
level = levelOfET;
depth_map = round(depth_map/max(abs(depth_map(:))) * level)* (max(abs(depth_map(:)))/level); %40 levels
for ch = 1:3
    deblurredImage(:, :, ch) = Deblurring(f1(:, :, ch), f2(:, :, ch), depth_map, inKernel1, inKernel2, inSigma);
end;

figure; imshow([f0, zeros(hei, 10, 3), deblurredImage]); title('Ground Truth Image (left) and Deblurred Image (right)');
imwrite(deblurredImage, [outFolder,  '\deblurredImage.png']);
