clear all;
close all;
mex MatchingPursuit.c;

%% set parameters 
initialScale = 1.4; %initial scale for Gabor wavelets
shrinkFactor = sqrt(2.0);
expandDoG = 14;

GarborScaleList=[initialScale, initialScale/shrinkFactor, initialScale/(shrinkFactor^2)];  % scales list for Gabor filters
DoGScaleList=[initialScale*expandDoG,  initialScale*expandDoG/shrinkFactor];  % scales list for DoG filters

numScale = length(GarborScaleList); %number of scales for Gabor wavelets
numDoG = length(DoGScaleList);
numOrient = 15;%number of orientations for Gabor wavelet


numSketch = 300; %number of sketches
locationShiftLimit = 3; %shift in normal direction
orientShiftLimit = 1; %shift in orientation

epsilon = .1; 
Upperbound = 6.; 
SHUTUP = 0.; 

sx=100;
sy=120;
%% read training images
imgFold = 'Image';
imgName = dir([imgFold '/*.jpg']);
numImg = size(imgName, 1);
clear I;
for (i = 1 : numImg)
    img = imread([imgFold '/' imgName(i).name]);
    img = imresize(img, [sx,sy],'nearest');
    img = im2double(rgb2gray(img));
    img=img-mean(img(:));
    img=img/std(img(:));
    I{i} = img;
end


%% filter images
tic;
[allFilterR, allFilterI, allSymbol] = MakeFilterBank(GarborScaleList, DoGScaleList, numOrient);
disp(['finished generating filter bank: '  num2str(toc) ' seconds']);

tic;
[Crr, Cri, Cir, Cii] = CorrFilter(allFilterR, allFilterI); %calculate correlation among filters
disp(['finished computing filters correlation: ' num2str(toc) ' seconds']);

tic;
[allFilteredR, allFilteredI, allFiltered] = ApplyFilter(I, allFilterR, allFilterI);
disp(['finished images filtering: ' num2str(toc) ' seconds']);

half = zeros(1, numOrient*numScale+numDoG); 
for (i=1:numOrient*numScale+numDoG)
    half(i) = (size(allFilterR{1, i}, 1)-1)/2;
end

%% prepare for storing experimental results
syma = zeros(sx, sy); 
sym = cell(1, numScale+1); 
for (i = 1 : numScale+1)
    sym{i} = zeros(sx, sy); 
end
Asym = cell(1, numImg); 
for (i = 1 : numImg)
    Asym{i} = zeros(sx, sy); 
end
Asyma = cell(1, numImg); 
for (i = 1 : numImg)
    Asyma{i} = zeros(sx, sy); 
end
Asym0 = cell(1, numImg*(numScale+1));
for (i = 1 : numImg*(numScale+1))
    Asym0{i} = zeros(sx, sy);  
end

filterSelected = zeros(1, numSketch);
xSelected = zeros(1, numSketch);
ySelected = zeros(1, numSketch);
numSketchSelected = 0;

%% do matching pursuit
tic;
MatchingPursuit(numImg, numOrient*numScale+numDoG, numOrient, numDoG, allFilteredR,...
                allFilteredI, allFiltered, Crr, Cri, Cir, Cii, half, sx, sy,...
                allFilterR(1, :), allFilterI(1, :), allSymbol(1, :), syma, sym,...
                Asym, Asym0, Asyma, locationShiftLimit, orientShiftLimit, numSketch,...
                epsilon, Upperbound, SHUTUP, filterSelected, xSelected, ySelected, numSketchSelected);
 
disp(['finished Matching Pursuit: ' num2str(toc) ' seconds']);
  
            
 filterSelected = filterSelected(1:numSketchSelected); 
 xSelected = xSelected(1:numSketchSelected);
 ySelected = ySelected(1:numSketchSelected);

 
 
 
%% show experimental results
output = 'output';
if ~exist(output,'dir')
   mkdir(output);
end
figure; 
col = 6;%each row has 6 images 
row = 1 + floor((numImg+1.)/2); 
for (i = 1:numScale+1) %show templates, the last one is for 
   subplot(row, col, i); imshow(-sym{i}, []);
end
subplot(row, col, numScale+2); imshow(-syma, []); 
showIndex = 6; 
for (i = 1 : numImg)
    showIndex = showIndex + 1;  subplot(row, col,  showIndex); 
    imshow(I{i}, []);       % all images
    showIndex = showIndex + 1;  subplot(row, col,  showIndex); 
    imshow(-Asyma{i}, []);  % sketch templates for all images
    showIndex = showIndex + 1;  subplot(row, col,  showIndex); 
    imshow(Asym{i}, []);    % reconstructed images
end

%% save experimental results
for (i = 1:numScale+1) %save learned templates
   towrite =  -sym{i}; 
   if max(towrite(:))-min(towrite(:))~=0
       towrite = (towrite-min(towrite(:)))/(max(towrite(:))-min(towrite(:)));
   end;
   imwrite(towrite, [output '/' sprintf('LearnedTemplate%d.png',i)],'PNG'); 
end
space = 4;
for (i = 1 : numImg)%save images, filters and recovered images
    h = size(I{i},1);
    w = size(I{i},2);
    towrite = ones(h, 4*w+3*space); % the ones make the background white, zeros make it black
    towrite(1:h,1:w) = (I{i}-min(I{i}(:)))/(max(I{i}(:))-min(I{i}(:)));%original image
    towrite(1:h,w+space+1:2*w+space) = (-Asyma{i}-min(-Asyma{i}(:)))/(max(-Asyma{i}(:))-min(-Asyma{i}(:)));%filters
    towrite(1:h,2*w+2*space+1:3*w+2*space) = (Asym{i}-min(Asym{i}(:)))/(max(Asym{i}(:))-min(Asym{i}(:)));%filters
    residual = I{i} - Asym{i};
    towrite(1:h,3*w+3*space+1:4*w+3*space) = (residual-min(residual(:)))/(max(residual(:))-min(residual(:)));%filters
    imwrite(towrite, [output '/' sprintf('ReconstructedImage%d.png',i)],'PNG');
end


%% generate html to show experimental results
GenerateHtml;
