
mex Ccopy.c;
mex sparseFRAME_SUM2.c

model = load(['./results/model.mat']);

category= model.category;
isLocalNormalize = model.isLocalNormalize;
isSeparate = model.isSeparate;

sigsq = model.sigsq;
nOrient = model.nOrient;

maxImArea=300*300;

localHalfx=20;
localHalfy=20;
thresholdFactor=0.01;
localNormScaleFactor=0.4;

toShowAndSaveMAX2SUM2map=true;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  parameters for detection task %%%%%%%%%%

flipOrNot=0;            % template flip Or not
rotateShiftLimit = 3;   % template rotation  from -rotateShiftLimit to rotateShiftLimit, eg. (1)-2:2 if rotateShiftLimit=2 (2)0 is without rotation
numResolution = 20;  % number of resolutions to search for in detection stage, the middle level is the original size
scaleStepSize = 0.1; % you can tune this either 0.1 or 0.2
initializedScaleFactor =200; %173;  % at the begining, we will rescale the testing image by factor initializedScaleFactor/ min([size(img,1), size(img,2)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

locationShiftLimit=model.locationShiftLimit;  % arg-max
orientShiftLimit=model.orientShiftLimit;    % arg-max

nScaleGabor = length(model.GarborScaleList);
nScaleDoG=length(model.DoGScaleList);
filters = model.filters;

%useDoG = model.useDoG;


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

%lambdaF = model.lambdaF;
sizeTemplatex = model.sx;
sizeTemplatey = model.sy;
halfTemplatex=floor(sizeTemplatex/2);
halfTemplatey=floor(sizeTemplatey/2);


inTestPath = ['./testingImage/' category '/'];



imgCell = cell(0);
files = dir([inTestPath '/*.jpg']);
%files=files(3);
numImage=length(files);
for iImg = 1:numImage
    img = imread(fullfile(inTestPath,files(iImg).name));
    %img = imresize(img,[sx, sy]);
    
    minSide=min([size(img,1), size(img,2)]);
    scaleFactor=initializedScaleFactor/minSide;
    %if size(img,1) * size(img,2) > maxImArea
        img = imresize(img, scaleFactor,'nearest');
    %end
    
    im=im2double(rgb2gray(img));
    im=im-mean(im(:));
    im=im/std(im(:))*sqrt(sigsq);
    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   


if ~exist(['detection_result/' category '/'],'dir')
       mkdir(['detection_result/' category '/']);
end
			

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 = applyfilterBank_MultiResolution_sparse(ImageMultiResolution, filters, halfFilterSizes, nOrient,...
        locationShiftLimit,orientShiftLimit,isLocalNormalize,isSeparate,localNormScaleFactor,thresholdFactor,nScaleGabor,nScaleDoG);
 
    mapName = ['detection_result/' category '/SUMMAXmap' 'image' num2str(img)];  
    save(mapName, 'ImageMultiResolution', 'SUM1mapFind','allSizex', 'allSizey');                
    disp(['filtering time: ' num2str(toc) ' seconds']);


end

template=model.template;

% for i=1:template.numSelected
%     lambdaF{template.selectedFilter(i)}(template.selectedRow(i),template.selectedCol(i)) = template.selectedLambdas(i);
% end

% 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 * nOrient *2  % Gabor
            
              iScale=floor((template.selectedFilter(i)-1)/nOrient); % index of iScale starts from 0;
              iOri= mod(template.selectedFilter(i)-1,nOrient);  % index of iOri starts from 0;
              
              if iOri>0
                  newOri = nOrient-iOri;
                  template.selectedFilter(i) = iScale*nOrient + newOri+1; % index of filter starts from 1; 
              end
          end
          
        end
    end
    
    for (rot = -rotateShiftLimit : rotateShiftLimit)
    
        theta = pi*rot/nOrient;
        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 * nOrient *2  % DoG
                    allSelectedFilter(r,i) = template.selectedFilter(i);
                
                else   % Gabor
                    iScale=floor((template.selectedFilter(i)-1)/nOrient); % index of iScale starts from 0;
                    iOri= mod(template.selectedFilter(i)-1,nOrient);  % index of iOri starts from 0;
                    newOri = iOri - rot;
                    if (newOri>=nOrient)
                        newOri = newOri - nOrient;
                    end
                    if (newOri<0)
                        newOri = newOri + nOrient;
                    end
                    allSelectedFilter(r,i) = iScale*nOrient + 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);

cropedImage=single(zeros(sizeTemplatex, sizeTemplatey));
for (img=1:numImage)
        
        disp(['======> start detecting in image ' num2str(img)]); tic
        
        mapName = ['detection_result/' category '/SUMMAXmap' 'image' num2str(img)]; 
        load(mapName);   % load SUM1mapFind 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}, SUM1mapFind);    
            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(['detection_result/' category '/' num2str(img) '.mat'],'allFx', 'allFy','MAX2score','MSUM2','Mrot','Mflip','MFx','MFy','Mind');
       
       Ccopy(cropedImage, single(ImageMultiResolution{Mind}), MFx, MFy, floor(sizeTemplatex/2), floor(sizeTemplatey/2), sizeTemplatex, sizeTemplatey, allSizex(Mind), allSizey(Mind), Mrot*pi/nOrient);
       if(Mflip > 0) 
           cropedImage=fliplr(cropedImage);
       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/nOrient;
       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,[ 'detection_result/' category '/detection-' num2str(img,'%04d') ],'png'); %name is a string
       
       
       %%%% save the croped aligned training images 
       gLow = min(cropedImage(:));
       gHigh = max(cropedImage(:));
       img_tem = (cropedImage-gLow)/(gHigh-gLow);
       %imwrite(img_tem,[ savingFolder 'training-' num2str(img,'%04d') '.png']);
       imwrite(img_tem,[ 'detection_result/' category '/cropped-' num2str(img,'%04d') '.png']);
       
       disp(['mex-C finding time: ' num2str(toc) ' seconds']);
      
       if toShowAndSaveMAX2SUM2map
       fig=figure;
       mesh(MapOfMAX2VsTemplateAndResolution);
       xlabel('resolution'); 
       ylabel('transformed template');
       zlabel('MAX2');
       title('MAX2 vs resolution and transformed template');
       saveas(fig,[ 'detection_result/' category '/MapOfMAX2VsTemplateAndResolution-' num2str(img,'%04d') ],'png'); %name is a string
       close;
       
       fig=figure;
       mesh(double(MSUM2{Mind}));
       colormap(gray);
       zlabel('SUM2');
       title('SUM2 map at the optimal pair of resolution and transformed template')
       shading interp;
       saveas(fig,[ 'detection_result/' category '/MaxSUM2-' num2str(img,'%04d') ],'png'); %name is a string
       close;
       end
       
 end

