
% REPRODUCING EXPERIMENT 1

clear; 
load 'nlf'; 

% CONSTANTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

epsilon = .1;  % allowed correlation between selected Gabors
Totalsketch = 60;  % allowed number of selected Gabors   
sub = 1;    % 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          


sizex = 100; sizey = 100;  % forced size of input image
imageFolder = 'positiveImage'; % folder of training images  
imageName = dir([imageFolder '\*.jpg']);
numImage = size(imageName, 1); % number of training images 
I = cell(1, numImage); 
for (img = 1 : numImage)
    tmpIm = imread([imageFolder '\' imageName(img).name]); 
    if size(tmpIm,3) == 3
        tmpIm = rgb2gray(tmpIm);
    end
    I{img} = imresize(double(tmpIm), [sizex sizey], 'nearest'); 
end
Sx = zeros(1, numImage) + sizex; 
Sy = zeros(1, numImage) + sizey; 
nimage = numImage; sx = sizex; sy = sizey; 
disp(['number of training images = ' num2str(numImage)]);

% FILTERING IMAGES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['start filtering']);
tic
[allfilter, allsymbol] = makefilter(scale, norient);  % generate Gabor filters 
h = (size(allfilter{1}, 1)-1)/2;  % half size of Gabor
C = corr(allfilter, epsilon);    % generate the inhibition maps
allfiltered = applyfilterfftsame(I, allfilter, localHalfx, localHalfy, thresholdFactor);  % filter training images
disp(['filtering time: ' num2str(toc) ' seconds']);

% the following three lines are the only things that need to be added
mex Ctransform.c; 
Ctransform(nimage, norient, allfiltered, Sx, Sy, Upperbound);
allfilteredback = mcopy(allfiltered, nimage, norient); 

% 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(sx, sy);  % symbolic plot of selected Gabor
Asym = cell(1, nimage);   % symbolic plot for each image with active Gabor
for (i = 1 : nimage)
    Asym{i} = zeros(sx, sy);  
end
tss = zeros(nimage, 1);  % total scores for testing images

binnum = 50; % number of bins 

numCluster = 2; 
prob = zeros(1, numCluster); 
tssAll = zeros(nimage, numCluster); 
dataWeightAll = rand(nimage,numCluster);
dataWeightAll = dataWeightAll./repmat(sum(dataWeightAll,2),1,numCluster); 
save dataWeightAll dataWeightAll; 
%load dataWeightAll; % remove this line if you do not want strict reproducibility
symAll = cell(1, numCluster); 
AsymAll = cell(nimage, numCluster); 

% MEX-C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mex Clearn.c;   % compile C code

numIt = 4; 
for (it = 1 : numIt)
disp(['M-step of iteration ' num2str(it)]);
for (c = 1:numCluster)
    tic
    dataWeight = dataWeightAll(:, c) + 0.; 
    Clearn(nimage, norient, allfiltered, allfilteredback, C, h, sx, sy, Mi, Mx, My, Mm, Mm1, allsymbol(1, :), sym, Asym, Lrange, Orange, sub, Totalsketch, Upperbound, tss, SHUTUP, binnum, M, lam, e, lz, dataWeight);
    disp(['   mex-C learning time for cluster ' num2str(c) ' takes ' num2str(toc) ' seconds']);
    tssAll(:, c) = tss + 0.; 
    symAll{c} = sym + 0.; 
    for (i = 1 : nimage)
        AsymAll{i, c} = Asym{i} + 0.; 
    end
end

disp(['E-step of iteration ' num2str(it)]);
for (c = 1 : numCluster)
    prob(c) = sum(dataWeightAll(:, c))/nimage; 
end
for (c = 1 : numCluster)
    for (i = 1 : nimage)
       dataWeightAll(i, c) = prob(c)/sum(prob(1, :).*exp(tssAll(i, :)-tssAll(i, c))); 
    end
end

if (it==numIt)
 disp(['Displaying the clusters']);
 for (cl = 1:numCluster)    
    K = nimage;
    tss = tssAll(:, cl); 
    [tsssort, ind] = sort(tss); 
    tt = 0; 
    for (i = 1 : K)
       this = ind(nimage-i+1); 
       if (dataWeightAll(this, cl)==max(dataWeightAll(this, :)))
          tt = tt + 1; 
          Iactu{tt} = I{this}; 
          Symactu{tt} = -(AsymAll{this, cl}); 
       end
    end

K = tt; bx =0; by = 0; cc = 3;
col = 6; row = floor(K/col)+1;
Iout = zeros(row*sx+(row-1)*(2*bx+cc), 2*col*sy+(col-1)*(2*by+cc))+255; 
towrite = -(symAll{cl}); 
towrite = (255 * (towrite-min(towrite(:)))/(max(towrite(:))-min(towrite(:))));
Iout(1:sx, 1:sy) = towrite;
i = -1;  
for (r = 1 : row)
    for (c = 1 : col)
        i = i+1; 
        if ((i>=1)&&(i<=K))
            towrite = Iactu{i}; 
            towrite = (255 * (towrite-min(towrite(:)))/(max(towrite(:))-min(towrite(:))));
            Iout((r-1)*sx+1+(r-1)*(2*bx+cc): r*sx+(r-1)*(2*bx+cc), (c-1)*2*sy+1+(c-1)*(2*by+cc): (c-1)*2*sy+sy+(c-1)*(2*by+cc)) = towrite; 
            towrite =  Symactu{i}; 
            towrite = (255 * (towrite-min(towrite(:)))/(max(towrite(:))-min(towrite(:))));
            Iout((r-1)*sx+1+(r-1)*(2*bx+cc) : r*sx+(r-1)*(2*bx+cc), (c-1)*2*sy+sy+(c-1)*(2*by+cc)+1 : (c-1)*2*sy+sy+(c-1)*(2*by+cc)+sy) =towrite;  
        end
    end
end
for (r = 1 : row-1)
  Iout(r*sx+(r-1)*(2*bx+cc)+bx+1:r*sx+(r-1)*(2*bx+cc)+bx+cc, :) = 0; 
end
for (c = 1 : col-1)
  Iout(:, c*(sy*2)+(c-1)*(2*by+cc)+by+1:c*(sy*2)+(c-1)*(2*by+cc)+by+cc) = 0; 
end
outepsgif(Iout, ['LNex4_cat_cluster' num2str(cl)]); 
 end

 for (c = 1:numCluster)   
     outepsgif(-(symAll{c}), ['LNex4_cat_clusterTemplates' num2str(c)]); 
 end
end
end
tssAll


