% main function for the FRAME model for unsupervised image alignment
clear;
close all;

inPath = './positiveImage/100+/cat';

%% Step 0: prepare filter and training images
mex Ccopy.c; % copy around detect location

mex CgetMAX1.c;
mex CgetMAX1_DoG.c;
mex SharedMatchingPursuit.c;
mex grad_U_sparse.c;
mex sparseFeatureExtraction.c;
mex sparseFRAME_SUM2.c;
mex ComputeProjectionV2.c;
mex sparseFRAME_SUM2_partial.c;

output_sparseFrame='outputSparseFrame';
if exist(output_sparseFrame)
   rmdir(output_sparseFrame,'s');
   mkdir(output_sparseFrame); % working directory to store synthesized images files. 
else
   mkdir(output_sparseFrame); % working directory to store synthesized images files. 
end

output_CroppedImages='outputCroppedImages';
if exist(output_CroppedImages)
   rmdir(output_CroppedImages,'s');
   mkdir(output_CroppedImages); % working directory to store synthesized images files. 
else
   mkdir(output_CroppedImages); % working directory to store synthesized images files. 
end

working_add='working_matSaving';
if exist(working_add)
   rmdir(working_add,'s');
   mkdir(working_add); % working directory to store .mat files. 
else
   mkdir(working_add); % working directory to store .mat files. 
end

isToUseNewFilterBankAtLastIteration=true;
isOption1=false;
%% parameters for filters (Gabor and DoG)
if isOption1

  %% option 1: automatical filters generation
  initialScale = 1.4; %initial scale for Gabor wavelets
  shrinkFactor = sqrt(2.0);
  expandDoG = 13.5;
  GaborScaleList=[initialScale, initialScale/shrinkFactor, initialScale/(shrinkFactor^2), initialScale/(shrinkFactor^3)];  % scales list for Gabor filters
  DoGScaleList=[];
  
  % if isToUseNewFilterBankAtLastIteration=true, we will use the following new filter bank to learn the FRAME model in the last iteration.
  GaborScaleListLastIteration=[initialScale, initialScale/shrinkFactor, initialScale/(shrinkFactor^2), initialScale/(shrinkFactor^3)];  % scales list for Gabor filters
  DoGScaleListLastIteration=[initialScale*expandDoG, initialScale*expandDoG/shrinkFactor];  % scales list for DoG filters
  
else
    
  %% option 2: manually filters generation
  DoGScaleList=[];
  GaborScaleList=[0.7];   % add scales of Gabor filters here
  
  % if isToUseNewFilterBankAtLastIteration=true, we will use the following new filter bank to learn the FRAME model in the last iteration.
  GaborScaleListLastIteration=[1.4, 1, 0.7, 0.5];   % add scales of Gabor filters here
  DoGScaleListLastIteration=[18.90, 13.36];
  
  %GaborScaleListLastIteration=[0.7];   % add scales of Gabor filters here
  %DoGScaleListLastIteration=[];
  
end
  nScaleGabor=length(GaborScaleList);
  nScaleDoG=length(DoGScaleList);
  numOrient = 16;  % number of orientations
  numFilter = nScaleGabor*2*numOrient+ nScaleDoG;  % sin and cosin separately
  numSketch= 300;  %20 number of Gabors (and DoG) at the first scale
  numSketchListLastIteration =300;  % if isToUseNewFilterBankAtLastIteration=true, we will use numSketchListLastIteration in filter selection in the last iteration.
  
%% parameters for HMC
sigsq=10;
lambdaLearningRate = 0.1/sqrt(sigsq);
nIter = 80; % 40;     % the number of iterations for learning lambda
epsilon = 0.03; % step size for the leapfrog
L = 10;         % leaps for the leapfrog
numSample = 3;  % how many HMC calls for each learning iteration
nTileRow = 6; %6;  % nTileRow x nTileCol defines the number of paralle chains
nTileCol = 6; %6;

%%
sizeTemplatex=150;    % template size x
sizeTemplatey=150;    % template size y
numIterationEM=2;     % EM iteration
flipOrNot=0;            % template flip Or not
rotateShiftLimit = 2;   % template rotation  from -rotateShiftLimit to rotateShiftLimit, eg. -2:2 if rotateShiftLimit=2
numResolution = 1;  % number of resolutions to search for in detection stage
scaleStepSize = 0.1; % you can tune this either 0.1 or 0.2
locationShiftLimit=3;
orientShiftLimit=1;
RatioDisplacementSUM2=0.25;   % default=0. Compute all values in SUM2 map


%% parameters for normalization
isGlobalNormalization=true; % for the input training image. the images substracts the mean and is divided by standard deviation
isLocalNormalizeFRAME=false;      % local normalization for feature map (SUM1) FRAME model
isLocalNormalizeFilterSelection=false;      % local normalization for addictive model
isSeparate=false;
localNormScaleFactor=0.5;
thresholdFactor=0.01;       % threshold of local normalization for feature map (SUM1)


%% create filters bank for Frame model (this filters bank will be used in learning Frame model)
disp('===============> Generating filters bank for Frame model');
tic;
filters=[];
% Gabor filters [sin cos sin cos ......]
for iScale=1:nScaleGabor
    f = MakeFilter(GaborScaleList(iScale),numOrient);
    for i=1:numOrient
     f_r{i} =real(f{i});
     f_i{i} =imag(f{i});
    end
    
    filters = [filters f_r f_i];
end
% DoG filter
for iScale=1:nScaleDoG
    f0 = single(dog(DoGScaleList(iScale),0));
    filters = [filters f0];
end

halfFilterSizes = zeros(size(filters));
for iF = 1:numFilter
    filters{iF} = single(filters{iF});
    halfFilterSizes(iF)=(size(filters{iF},1)-1)/2;
end

Corr = CorrFilterFrame(filters); %calculate correlation among filters
disp(['finished creating filters bank: ' num2str(toc) ' seconds']);

%% create filters bank for sparse coding model (this filters bank will be used in filter selection)
disp('===============> Generating filters bank for sparse coding model');
tic;
[allFilterR, allFilterI, allSymbol] = MakeFilterBank(GaborScaleList, DoGScaleList, numOrient);
[Crr, Cri, Cir, Cii] = CorrFilter(allFilterR, allFilterI); %calculate correlation among filters

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

 % all necessary variables for filter selection will be stored in structure variable "FitersInformation", which will be sent into the filter selection function
FiltersInformation.allFilterR = allFilterR;  
FiltersInformation.allFilterI = allFilterI;
FiltersInformation.allSymbol = allSymbol;
FiltersInformation.Crr = Crr;
FiltersInformation.Cri = Cri;
FiltersInformation.Cir = Cir;
FiltersInformation.Cii = Cii;
FiltersInformation.half = half;
disp(['finished creating filters bank: ' num2str(toc) ' seconds']);


halfTemplatex=floor(sizeTemplatex/2);
halfTemplatey=floor(sizeTemplatey/2);


imgCell = cell(0);
files = dir([inPath '/*.jpg']);
numImage=length(files);

% % padding_x=round(halfTemplatex/3);
% % padding_y=round(halfTemplatey/3);

for iImg = 1:numImage
    
    img = imread(fullfile(inPath,files(iImg).name));      
        
    img = imresize(img,[sizeTemplatex, sizeTemplatey]);
    %img = imresize(img,[sizeTemplatex+padding_x, sizeTemplatey+padding_y]);
    im=im2double(rgb2gray(img));
                  %im=padarray(im,[padding_x padding_y],'replicate');
    if isGlobalNormalization
       im=im-mean(im(:));
       im=im/std(im(:))*sqrt(sigsq);
    end
    
    imgCell{iImg}=im;    
    
end

%% Generate multi-resolution images
originalResolution = round(numResolution/2);    %3; % original resolution is the one at which the imresize factor = 1, see 11th line beneath this line 
allSizex = zeros(1, numResolution); allSizey = zeros(1, numResolution); 
ImageMultiResolution = cell(1, numResolution); 
MAX2score = single(zeros(1, numResolution));  % maximum log-likelihood score at each resolution
allFx = zeros(1, numResolution); allFy = zeros(1, numResolution); % detected location at each resolution   

for (img = 1:numImage)
  
  for(resolution=1:numResolution)
    resizeFactor = 1.0 + (resolution - originalResolution)*scaleStepSize; 
    ImageMultiResolution{resolution} = imresize(imgCell{img}, resizeFactor, 'nearest');  % images at multiple resolutions   
    [sizex, sizey] = size(ImageMultiResolution{resolution}); 
    allSizex(resolution) = sizex; 
    allSizey(resolution) = sizey; 
    
  end
  disp(['======> start filtering and maxing image ' num2str(img)]); tic
     
  [SUM1mapFind, MAX1mapFind]= applyfilterBank_MultiResolution_sparse(ImageMultiResolution, filters, halfFilterSizes, numOrient,...
      locationShiftLimit,orientShiftLimit,isLocalNormalizeFRAME,isSeparate,localNormScaleFactor,thresholdFactor,nScaleGabor,nScaleDoG);  
  
  mapName = [working_add '/SUMMAXmap' 'image' num2str(img)];   
  save(mapName, 'ImageMultiResolution', 'SUM1mapFind','MAX1mapFind','allSizex', 'allSizey');                
  disp(['filtering time: ' num2str(toc) ' seconds']);
end


%% Initialize learning and compute training sample averages
savingFolder=[output_CroppedImages '/iteration0/'];
if ~exist(savingFolder,'dir')
       mkdir(savingFolder);
end


cropedImage=cell(1,numImage); 
SUM1mapLearn = cell(numImage,numFilter);
MAX1mapLearn = cell(numImage,numFilter);


cropedImage=cell(1,numImage);
for i=1:numImage
   cropedImage{1,i}=single(zeros(sizeTemplatex, sizeTemplatey));
end

disp(['Start initialization']);  
for (img = 1:numImage) 
    
  
   mapName = [ working_add '/SUMMAXmap' 'image' num2str(img)]; 
   load(mapName);   % load SUM1mapFind and ImageMultiResolution for imgage img     
 
  cropedImage{img}= imresize(single(ImageMultiResolution{originalResolution}), [sizeTemplatex, sizeTemplatey]);
  
  img_scaled = cropedImage{img};
  img_scaled = img_scaled - mean(img_scaled(:));
  img_scaled = img_scaled/std(img_scaled(:))*sqrt(sigsq);
  [SUM1mapLearn(img, :), MAX1mapLearn(img, :)]= applyfilterBank_MultiResolution_sparse({img_scaled}, filters, halfFilterSizes, numOrient,...
                      locationShiftLimit,orientShiftLimit,isLocalNormalizeFRAME,isSeparate,localNormScaleFactor,thresholdFactor,nScaleGabor,nScaleDoG);  
  
  %%%% save the croped aligned training images 
  imgSaved=cropedImage{img};
  gLow = min(imgSaved(:));
  gHigh = max(imgSaved(:));
  imgSaved = (imgSaved-gLow)/(gHigh-gLow);
  imwrite(imgSaved,[ savingFolder 'iteration-0-' 'training-' num2str(img,'%04d') '.png']);
end

for iFilter = 1:numFilter
    rHat{iFilter}=zeros(sizeTemplatex, sizeTemplatey,'single');
end
    
for img = 1:numImage
    for iFilter = 1:numFilter
       rHat{iFilter}= rHat{iFilter}+ MAX1mapLearn{img,iFilter}/numImage;
    end
end


savingFolder=[output_sparseFrame '/iteration0/'];
if exist(savingFolder)
   rmdir(savingFolder,'s'); 
   mkdir(savingFolder);
else
   mkdir(savingFolder);
end
               
 % filter selection
ticID1=tic;                        
[template, templateMP]= filtersSelection_codebookV2(cropedImage, FiltersInformation, nScaleGabor, nScaleDoG, numOrient, numSketch,...
       locationShiftLimit, orientShiftLimit, double(sizeTemplatex),double(sizeTemplatey), isLocalNormalizeFilterSelection, localNormScaleFactor, thresholdFactor, savingFolder);
template.selectedLambdas=single(zeros(1,template.numSelected));            
disp(['===============> finished filters selection: '  num2str(toc(ticID1)) ' seconds']);

% sparse FRAME learning
ticID2=tic;
isSaved=1;               
[template,currSample,logZ]=sparseFRAMElearn(template, nIter, filters, rHat, double(sizeTemplatex),double(sizeTemplatey), halfFilterSizes, ...
        locationShiftLimit, nTileRow,nTileCol,epsilon,L,lambdaLearningRate,numSample, isSaved,savingFolder);
disp(['===============> finished filters selection: '  num2str(toc(ticID2)) ' seconds']);
                              
%save([working_add sprintf( '/learnedmodel_iter%d.mat',iter)], 'template', 'logZ', 'currSample', 'templateMP'); 

%%% save the synthesized results
saveSynImages(currSample, sizeTemplatex, sizeTemplatey, nTileCol, nTileRow, Corr, template, single(halfFilterSizes), filters, savingFolder);
   

%%%%%%%%%%%%%%%% ratate and flip template (lambda) with multiple scales of Gabor filters 
disp(['===============> transforming template']);
%load([working_add sprintf( '/learnedmodel_iter%d.mat',iter)], 'template', 'logZ');           


% shift the origin to the center of the template
template.selectedRow = template.selectedRow - round(sizeTemplatex/2) +1;
template.selectedCol = template.selectedCol - round(sizeTemplatey/2) +1;

allSelectedFilter=[];
allSelectedRow=[];
allSelectedCol=[];
%%% rotate or flip lambda maps

numElement=template.numSelected;
for (flip = 0 : flipOrNot)
    
    if (flip > 0)
        
        for (i=1: numElement)
        
          template.selectedCol(i) = -template.selectedCol(i);
        
          if  template.selectedFilter(i)<= nScaleGabor * numOrient *2  % Gabor
            
              iScale=floor((template.selectedFilter(i)-1)/numOrient); % index of iScale starts from 0;
              iOri= mod(template.selectedFilter(i)-1,numOrient);  % index of iOri starts from 0;
              
              if iOri>0
                  newOri = numOrient-iOri;
                  template.selectedFilter(i) = iScale*numOrient + newOri+1; % index of filter starts from 1; 
              end
          end
          
        end
    end
    
    for (rot = -rotateShiftLimit : rotateShiftLimit)
    
        theta = pi*rot/numOrient;
        sintheta = sin(theta);
        costheta = cos(theta);
            
        r = rot+rotateShiftLimit+1  + (rotateShiftLimit*2+1)*flip;  
        
        if (rot == 0)
            for (i=1: numElement)
                allSelectedRow(r,i) = template.selectedRow(i);
                allSelectedCol(r,i) = template.selectedCol(i);
                allSelectedFilter(r,i) = template.selectedFilter(i);                      
            end
        else
            for (i=1: numElement)
                allSelectedRow(r,i) = template.selectedRow(i)*costheta + template.selectedCol(i)*sintheta;
                allSelectedCol(r,i) = -template.selectedRow(i)*sintheta + template.selectedCol(i)*costheta;
                
                if  template.selectedFilter(i)> nScaleGabor * numOrient *2  % DoG
                    allSelectedFilter(r,i) = template.selectedFilter(i);
                
                else   % Gabor
                    iScale=floor((template.selectedFilter(i)-1)/numOrient); % index of iScale starts from 0;
                    iOri= mod(template.selectedFilter(i)-1,numOrient);  % index of iOri starts from 0;
                    newOri = iOri - rot;
                    if (newOri>=numOrient)
                        newOri = newOri - numOrient;
                    end
                    if (newOri<0)
                        newOri = newOri + numOrient;
                    end
                    allSelectedFilter(r,i) = iScale*numOrient + newOri+1; % index of filter starts from 1;                        
                end              
            end       
        end       
    end

end

numTransform=(2*rotateShiftLimit+1)*(1+flipOrNot);
AllTemplate=cell(1,numTransform);
for r=1:numTransform
   AllTemplate{r}=struct('selectedRow',single(allSelectedRow(r,:)),'selectedCol',single(allSelectedCol(r,:)),'selectedFilter',single(allSelectedFilter(r,:)),'selectedLambdas',template.selectedLambdas); 
end


MapOfMAX2VsTemplateAndResolution=zeros(numTransform,numResolution);
%%%%% EM 
for(it = 1:numIterationEM)
    
    savingFolder=[output_CroppedImages '/iteration' num2str(it) '/'];
    if ~exist(savingFolder,'dir')
       mkdir(savingFolder);
    end
    
    disp(['**************** detection for iteration ' num2str(it)]);
    
    croppedImagesScaled=cell(1,numImage);
    for (img=1:numImage)
        disp(['======> start detecting in image ' num2str(img)]); tic
        
        mapName = [working_add '/SUMMAXmap' 'image' num2str(img)]; 
        load(mapName);   % load SUM1mapFind, MAX1mapFind and ImageMultiResolution for image img
  
        
        MMAX2 = -1e10; 
        
        for (flip = 0 : flipOrNot)
          for (rot = -rotateShiftLimit : rotateShiftLimit)
            
             r = rot+rotateShiftLimit+1 + (rotateShiftLimit*2+1)*flip;  
                      
            %[SUM2] = sparseFRAME_SUM2( single(numResolution), single(allSizex), single(allSizey), single(numFilter), AllTemplate{r}, MAX1mapFind);    
            [SUM2] = sparseFRAME_SUM2_partial( single(numResolution), single(allSizex), single(allSizey), single(numFilter), AllTemplate{r}, MAX1mapFind, single(RatioDisplacementSUM2));    
            for iResolution=1:numResolution
               [MAX2score(iResolution),ind]=max(SUM2{iResolution}(:));
               [allFx(iResolution),  allFy(iResolution)] = ind2sub([allSizex(iResolution), allSizey(iResolution)],ind);
            end
            
            MapOfMAX2VsTemplateAndResolution(r,:)=MAX2score;
            
            
            [maxOverResolution, ind] = max(MAX2score);   % most likely resolution
            if (MMAX2 < maxOverResolution) 
               MMAX2 = maxOverResolution; Mrot = rot; Mflip = flip; 
               Mind = ind; MFx = allFx(ind); MFy = allFy(ind); 
               MSUM2 = SUM2;
            end
          end
        end
        
        save([working_add '/' num2str(it) '-' num2str(img) '.mat'],'allFx', 'allFy','MAX2score','MSUM2','Mrot','Mflip','MFx','MFy','Mind');
        
      %%  save(['detection_result/' category '/' num2str(img) '.mat'],'allFx', 'allFy','MAX2score','MSUM2','Mrot','Mflip','MFx','MFy','Mind');
   
       Ccopy(cropedImage{img}, single(ImageMultiResolution{Mind}), MFx, MFy, floor(sizeTemplatex/2), floor(sizeTemplatey/2), sizeTemplatex, sizeTemplatey, allSizex(Mind), allSizey(Mind), Mrot*pi/numOrient);
       
       if(Mflip > 0) 
           cropedImage{img}=fliplr(cropedImage{img});
       end
       
       scaledImg=ImageMultiResolution{Mind};
       rect_Vertices=[-halfTemplatex -halfTemplatey; -halfTemplatex halfTemplatey; halfTemplatex halfTemplatey; halfTemplatex -halfTemplatey; -halfTemplatex -halfTemplatey];
       top_left_Vertices=[-halfTemplatex-0.1*halfTemplatex -halfTemplatey; -halfTemplatex-0.1*halfTemplatex -halfTemplatey+0.1*halfTemplatey; -halfTemplatex -halfTemplatey+0.1*halfTemplatey];
       combined_Vertices=[rect_Vertices; top_left_Vertices];
       if Mflip==1 
           combined_Vertices(:,2)=-1*combined_Vertices(:,2); 
       end
       
       theta = Mrot*pi/numOrient;
       R=[cos(theta) -sin(theta);sin(theta) cos(theta)];    %clockwise rotate theta
       
       new_rect_Vertices=combined_Vertices*R+ones(size(combined_Vertices,1),1)*[MFx, MFy];
       fig=figure; imshow(scaledImg,[]); hold on;
       plot(new_rect_Vertices(:,2),new_rect_Vertices(:,1),'LineWidth',1.5);
   
       saveas(fig,[savingFolder '/detected-iteration-' num2str(it) '-training-' num2str(img,'%04d') ],'png'); %name is a string
       
       
       %%%% save the croped aligned training images 
       imgSaved=cropedImage{img};
       gLow = min(imgSaved(:));
       gHigh = max(imgSaved(:));
       imgSaved = (imgSaved-gLow)/(gHigh-gLow);
       %imwrite(img_tem,[ savingFolder 'training-' num2str(img,'%04d') '.png']);
       imwrite(imgSaved,[ savingFolder '/iteration-' num2str(it) '-training-' num2str(img,'%04d') '.png']);
       
       disp(['finding time: ' num2str(toc) ' seconds']);
        
       
       img_scaled = cropedImage{img};
       img_scaled = img_scaled - mean(img_scaled(:));
       img_scaled = img_scaled/std(img_scaled(:))*sqrt(sigsq);
       
       croppedImagesScaled{img}=img_scaled;
       
    end
    close all; % close all the detected images
%% Iteration part 2: re-learn the template
      
 
    if(it == numIterationEM && isToUseNewFilterBankAtLastIteration)       
           useNewFilterBankLastIteration;       
    end
    for img=1:numImage
        [SUM1mapLearn(img, :), MAX1mapLearn(img, :)]= applyfilterBank_MultiResolution_sparse(croppedImagesScaled(img), filters, halfFilterSizes, numOrient,...
                      locationShiftLimit,orientShiftLimit,isLocalNormalizeFRAME,isSeparate,localNormScaleFactor,thresholdFactor,nScaleGabor,nScaleDoG); 
    end          
        
    disp(['**************** multiple image learning for iteration ' num2str(it)]); tic
    for iFilter = 1:numFilter
        rHat{iFilter}=zeros(sizeTemplatex, sizeTemplatey,'single');
    end
    
    for img = 1:numImage
        for iFilter = 1:numFilter
           rHat{iFilter}= rHat{iFilter}+ MAX1mapLearn{img,iFilter}/numImage;
        end
    end
    
    savingFolder=[output_sparseFrame '/iteration' num2str(it) '/'];
    if exist(savingFolder)
       rmdir(savingFolder,'s'); 
       mkdir(savingFolder);
    else
       mkdir(savingFolder);
    end
               
    % filter selection
    ticID1=tic;                        
    [template, templateMP]= filtersSelection_codebookV2(cropedImage, FiltersInformation, nScaleGabor, nScaleDoG, numOrient, numSketch,...
       locationShiftLimit, orientShiftLimit, double(sizeTemplatex),double(sizeTemplatey), isLocalNormalizeFilterSelection, localNormScaleFactor, thresholdFactor, savingFolder);
    template.selectedLambdas=single(zeros(1,template.numSelected));            
    disp(['===============> finished filters selection: '  num2str(toc(ticID1)) ' seconds']);
               
            
    % sparse FRAME learning
    ticID2=tic;
    isSaved=1;               
    [template,currSample,logZ]=sparseFRAMElearn(template, nIter, filters, rHat, double(sizeTemplatex),double(sizeTemplatey), halfFilterSizes, ...
        locationShiftLimit, nTileRow,nTileCol,epsilon,L,lambdaLearningRate,numSample, isSaved,savingFolder);
    disp(['===============> finished filters selection: '  num2str(toc(ticID2)) ' seconds']);   
    
    % save the synthesized images
    saveSynImages(currSample, sizeTemplatex, sizeTemplatey, nTileCol, nTileRow, Corr, template, single(halfFilterSizes), filters, savingFolder);   
    
    % shift the origin to the center of the template
    template.selectedRow = template.selectedRow - round(sizeTemplatex/2) +1;
    template.selectedCol = template.selectedCol - round(sizeTemplatey/2) +1;


    allSelectedFilter=[];
    allSelectedRow=[];
    allSelectedCol=[];
    %%% rotate or flip lambda maps

    numElement=template.numSelected;
    for (flip = 0 : flipOrNot)    
      if (flip > 0)        
        for (i=1: numElement)
        
          template.selectedCol(i) = -template.selectedCol(i);
        
          if  template.selectedFilter(i)<= nScaleGabor * numOrient *2  % Gabor
            
              iScale=floor((template.selectedFilter(i)-1)/numOrient); % index of iScale starts from 0;
              iOri= mod(template.selectedFilter(i)-1,numOrient);  % index of iOri starts from 0;
              
              if iOri>0
                  newOri = numOrient-iOri;
                  template.selectedFilter(i) = iScale*numOrient + newOri+1; % index of filter starts from 1; 
              end
          end
          
        end
      end
    
      for (rot = -rotateShiftLimit : rotateShiftLimit)
    
        theta = pi*rot/numOrient;
        sintheta = sin(theta);
        costheta = cos(theta);
            
        r = rot+rotateShiftLimit+1  + (rotateShiftLimit*2+1)*flip;  
        
        if (rot == 0)
            for (i=1: numElement)
                allSelectedRow(r,i) = template.selectedRow(i);
                allSelectedCol(r,i) = template.selectedCol(i);
                allSelectedFilter(r,i) = template.selectedFilter(i);                      
            end
        else
            for (i=1: numElement)
                allSelectedRow(r,i) = template.selectedRow(i)*costheta + template.selectedCol(i)*sintheta;
                allSelectedCol(r,i) = -template.selectedRow(i)*sintheta + template.selectedCol(i)*costheta;
                
                if  template.selectedFilter(i)> nScaleGabor * numOrient *2  % DoG
                    allSelectedFilter(r,i) = template.selectedFilter(i);
                
                else   % Gabor
                    iScale=floor((template.selectedFilter(i)-1)/numOrient); % index of iScale starts from 0;
                    iOri= mod(template.selectedFilter(i)-1,numOrient);  % index of iOri starts from 0;
                    newOri = iOri - rot;
                    if (newOri>=numOrient)
                        newOri = newOri - numOrient;
                    end
                    if (newOri<0)
                        newOri = newOri + numOrient;
                    end
                    allSelectedFilter(r,i) = iScale*numOrient + newOri+1; % index of filter starts from 1;                        
                end              
            end       
        end       
      end

   end

   numTransform=(2*rotateShiftLimit+1)*(1+flipOrNot);
   AllTemplate=cell(1,numTransform);
   for r=1:numTransform
       AllTemplate{r}=struct('selectedRow',single(allSelectedRow(r,:)),'selectedCol',single(allSelectedCol(r,:)),'selectedFilter',single(allSelectedFilter(r,:)),'selectedLambdas',template.selectedLambdas); 
   end

end

% save the results as html
fid = fopen([output_CroppedImages '/result.html'], 'wt'); 
for it=0:numIterationEM
  fprintf(fid, '%s\n', ['iteration-' num2str(it) '<br>'] ); 
 
  savingFolder=['iteration' num2str(it) '/'];
  for img=1:numImage
      fprintf(fid, '%s\n', ['<IMG SRC="' savingFolder 'iteration-' num2str(it) '-training-' num2str(img,'%04d') '.png"' ' height=70 ' 'width=70>'] ); 
  end

  fprintf(fid, '%s\n', ['<br>']);
  fprintf(fid, '%s\n', ['<hr>']);
end

fclose(fid);

