
% LOAD IN tilt RESULTS and MEX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; 
load 'nlf'; 
mex Ctransform.c; 
mex Clearn.c;   
mex Cfind.c;   
mex Ccopy.c; 
mex Cmax.c; 


% CONSTANTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsilon = .1;  % allowed correlation between selected Gabors
sub = 2;    % subsampling rate for Gabor
Lrange = 3; % allowed shifting in location (times sub pixels)
Orange = 1; % allowed shifting in orientation
SHUTUP = 0.; % SHUTUP is to neglect the positions with low average responses               
binnum = 50; % number of bins 
nstep = 3;  % number of iterations

% FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C = corr(allfilter, epsilon);    % generate the inhibition maps

% TRAINING IMAGES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fold = .8; 
Totalsketch = 60;  % allowed number of selected Gabors   
% TRAINING IMAGES, only the first one needs to be labeled with bounding box
Iname = {'dear1.jpg' 'dear2.jpg' 'dear3.jpg' 'dear4.jpg' 'dear5.jpg'};
Totalsketch = 60;  % allowed number of selected Gabors   
sx0 = 168; sy0 = 185;  % size of bounding box
xhere = 22;  yhere = 45; ihere = 0;  % position of the bounding box in image ihere
fold = .8; 
sx0 = floor(sx0*fold+.5); sy0 = floor(sy0*fold+.5); 
xhere = floor(xhere*fold+.5); yhere = floor(yhere*fold+.5); 

nimage = size(Iname, 2);
Sx = zeros(1, nimage); Sy = zeros(1, nimage); 
for (i = 1 : nimage)
    I{i} = randn(sx0, sy0); 
    J{i} = imresize(double(rgb2gray(imread(Iname{i}))), fold, 'nearest'); 
    [sx1, sy1] = size(J{i}); 
    Sx(i) = sx1; Sy(i) = sy1; % sizes of images
end
Sx0 = zeros(1, nimage) + sx0; 
Sy0 = zeros(1, nimage) + sy0; 

% VARIABLES FOR MEX-C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mi = zeros(1, Totalsketch);  % orientation of selected Gabor
Mx = zeros(1, Totalsketch);  % position of selected Gabor
My = zeros(1, Totalsketch); 
Mm = zeros(1, Totalsketch);  % lambda
Mm1 = zeros(1, Totalsketch);  % logZ
sym = zeros(size(I{1}));  % symbolic plot of selected Gabor
Asym = cell(1, nimage);   % symbolic plot for each image with active Gabor
for (i = 1 : nimage)
    Asym{i} = zeros(size(I{1}));  
end

% START FROM SINGLE IMAGE LEARNING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
J0 = J{1}; 
I0{1} = J0(xhere:(xhere+sx0-1), yhere:(yhere+sy0-1)); 
nimage0 = 1; Lrange0 = 0; Orange0 = 0; 
allfiltered0 = applyfilterfftsame(I0, allfilter); 
Asym0{1} = zeros(size(I{1}));
Ctransform(nimage0, norient, allfiltered0, Sx0, Sy0, Upperbound);
disp(['start from single image learning']);
tic
Clearn(nimage0, norient, allfiltered0, C, h, sx0, sy0, Mi, Mx, My, Mm, Mm1, allsymbol(1, :), sym, Asym0, Lrange0, Orange0, sub, Totalsketch, Upperbound, SHUTUP, binnum, M, lam, e, lz);
disp(['mex-C learning time: ' num2str(toc) ' seconds']);
% figure; imshow(-sym, []); 

% TRAINING IMAGES AT DIFFERENT SCALES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
totalScale = 9;  % total number of scales
ax = floor(sx0/2); ay = floor(sy0/2); % half size of bounding box
Sxx = zeros(nimage, totalScale); Sxa = zeros(1, totalScale*nimage); 
Syy = zeros(nimage, totalScale); Sya = zeros(1, totalScale*nimage); 
mscale = zeros(nimage, totalScale);  % log-likelihood score over scales
JJ = cell(1, totalScale*nimage); 
tss1 = cell(1, totalScale*nimage); 
sym1 = cell(1, totalScale*nimage); 

for (i = 1:nimage)
  for(j=1:totalScale)
    JJ{(i-1)*totalScale+j} = imresize(J{i}, .8+(j-1)*.05, 'nearest');  % scaled images
    [sx1, sy1] = size(JJ{(i-1)*totalScale+j}); 
    Sxx(i, j) = sx1; Syy(i, j) = sy1; % sizes of images
    Sxa((i-1)*totalScale+j) = sx1; Sya((i-1)*totalScale+j) = sy1; 
    tss1{(i-1)*totalScale+j} = zeros(sx1, sy1);  % log-likelihood score map for testing images
    sym1{(i-1)*totalScale+j} = zeros(sx1, sy1);  % symbolic sketch on testing image
  end
end

disp(['start filtering training images at all scales']);
tic
allfiltered1 = applyfilterfft(JJ, allfilter);  % filter training images
disp(['filtering time: ' num2str(toc) ' seconds']);

Ctransform(totalScale*nimage, norient, allfiltered1, Sxa, Sya, Upperbound);
allmaxed = cell(totalScale*nimage, norient);
for (i = 1:totalScale*nimage)
    for (o = 1:norient)
        allmaxed{i, o} = allfiltered1{i, o} + 0.; 
    end
end
sub1 = 1; Lrange1 = 6; 
Cmax(totalScale*nimage, norient, allmaxed, Sxa, Sya, Lrange1, Orange, sub1);

% ITERATE DETECTION AND LEARNING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for(it = 1:nstep)    
disp(['detection for step ' num2str(it)]);
for (i=1:nimage)
disp(['    start detecting in image ' num2str(i)]);
tic
msci = zeros(1, totalScale); Fx = zeros(1, totalScale); Fy = zeros(1, totalScale); 
Sxi = Sxx(i, :); Syi = Syy(i, :); 
Cfind(i, nimage, totalScale, norient, allfiltered1, allmaxed, h, ax, ay, Mi, Mx, My, Mm, Mm1, allsymbol(1, :), Lrange1, Orange, sub1, Totalsketch, Upperbound, tss1, sym1, Sxi, Syi, msci, Fx, Fy);
disp(['    mex-C detection time: ' num2str(toc) ' seconds']);

maxscale = max(msci);   % most likely scale
for(j=1:totalScale)
    if(msci(j)==maxscale)
       Jshow{i} = JJ{(i-1)*totalScale+j}+sym1{(i-1)*totalScale+j}*100;  % overlaid sketch 
       Ccopy(I{i}, JJ{(i-1)*totalScale+j}, floor(Fx(j)+.5), floor(Fy(j)+.5), sx0, sy0, Sxi(j), Syi(j)); 
    end
end
end

allfiltered = applyfilterfftsame(I, allfilter); 
Ctransform(nimage, norient, allfiltered, Sx0, Sy0, Upperbound);
disp(['multiple image learning for step ' num2str(it)]);
tic
Clearn(nimage, norient, allfiltered, C, h, sx0, sy0, Mi, Mx, My, Mm, Mm1, allsymbol(1, :), sym, Asym, Lrange, Orange, sub, Totalsketch, Upperbound, SHUTUP, binnum, M, lam, e, lz);
disp(['mex-C learning time: ' num2str(toc) ' seconds']);
end

outepsgif(-sym, 'ex5_deer_sym'); 
for(i=1:nimage)
       outepsgif(Jshow{i}, ['ex5_deer' num2str(i)]);  % overlaid sketch 
end

