classdef PhowExperience < handle
    %USE THIS CLASS TO RUN EXPERIMENTS ON THE 15-SCENES AND ISR DATASETS.
    %
    %%USAGE EXAMPLE:
    % dtsetFile='/path/to/15ScenesDataset';
    %
    % maxImgSz=300;
    % patchSize=24;
    % spacing=32;
    % K=100;
    % L=0;
    % kern=1;
    % nTrainingImg=30;
    %
    % exp=PhowExperience(dtsetFile,maxImgSz,patchSize,spacing,K,L,kern,nTrainingImg);
    %
    % C=100;
    % splitNumber=1;
    %
    % [accuracy,model,dec_values,predict_labels,confusion]=exp.runExperiment(C,splitNumber);
    
    
    properties (Constant)
        VLFEATPATH='/path/to/vlfeat';           %path to the vl-feat library
        LIBSVMPATH='/path/to/libsvm/matlab';    %path to the lib-svm solver matlab interface
        KERNAMES={'LINEAR','INTERSECTION','CHI-SQUARE','GAUSSIAN'};
    end
    
    properties
        %GENERAL SETTINGS
        doFromScratch       % if is set to true, the experiment is run from scratch, discarding the precomputed data
        
        %EXPERIMENT PATHS
        datasetPath         % the path where the dataset is stored  (a directory containing another directory called images, where the images are stored)
        imgPath             % by default is datasetPath/images
        featPath            % by default is datasetPath/features
        kernelsPath         % by default is datasetPath/kernels
        dictionaryPath      % by default is datasetPath/dictionaries
        resultsPath         % by default is datasetPath/results
        splitFile           % by default is datasetPath/splits_NTR..
        
        %SIFT SETTINGS
        maxImageSize        % the maximum size allowed for an image
        patchSize           % the size of the SIFT patch
        spacing             % spacing between SIFT points (for dense sampling)
        
        %DICTIONARY LEARNING SETTINGS
        nVocFeats           % the maximum number of local features used to construct the vocabulary
        K                   % the number of visual words
        
        %PHOW ENCODING SETTINGS
        L                   % the pyramidLevel
        pyramidWeights      % boolean value stating whether the spatial pyramid must be normalized using the weights, as in the paper
        
        %KERNEL SETTINGS
        kernelType          % specifies the type of kernel: 1 for linear, 2 for intersection, 3 for chi square, 4 for gaussian
        
        %CROSS VALIDATION SETTINGS
        nTrain              % the number of training samples
        
    end
    
    methods
        
        function newExperiment=PhowExperience(dataset_path,max_image_size,patch_size,spacing,K,L,kernelType,n_train)
            %CONSTRUCTOR OF THE EXPERIMENT
            %
            %PARAMETERS
            % dataset_path is the path where the dataset is stored
            %
            % max_image_size is the maximum size allowed for an image
            %
            % patch_size is the size of the SIFT patch
            %
            % spacing is spacing between SIFT points (for dense sampling)
            %
            % K is the number of visual words
            %
            % L is the pyramidLevel
            %
            % kernelType specifies the type of kernel: 1 for linear, 2 for
            % intersection, 3 for chi square, 4 for gaussian
            %
            % n_train is the number of training samples
            
            n_voc_feats=5e5;
            pyramid_weights=false;
            
            img_path=sprintf('%s/images',dataset_path);
            feat_path=sprintf('%s/features',dataset_path);
            kernels_path=sprintf('%s/kernels',dataset_path);
            dictionary_path=sprintf('%s/dictionaries',dataset_path);
            results_path=sprintf('%s/results',dataset_path);
            split_file=sprintf('%s/splits',dataset_path);
            
            newExperiment.doFromScratch=false;
            
            newExperiment.datasetPath=dataset_path;
            newExperiment.imgPath=img_path;
            newExperiment.featPath=feat_path;
            newExperiment.kernelsPath=kernels_path;
            newExperiment.dictionaryPath=dictionary_path;
            newExperiment.resultsPath=results_path;
            newExperiment.splitFile=split_file;
            
            newExperiment.maxImageSize=max_image_size;
            newExperiment.patchSize=patch_size;
            newExperiment.spacing=spacing;
            newExperiment.K=K;
            newExperiment.L=L;
            newExperiment.pyramidWeights=pyramid_weights;
            
            newExperiment.kernelType=kernelType;
            
            newExperiment.nTrain=n_train;
            newExperiment.nVocFeats=n_voc_feats;
            
            %ADDS VL_FEAT TO THE PATH
            run(sprintf('%s/toolbox/vl_setup.m',PhowExperience.VLFEATPATH));
            addpath(PhowExperience.LIBSVMPATH);
            
            %CREATE DIRS
            newExperiment.createDirs();
        end
        
        
        function createDirs(obj)
            %CREATES THE DATASET DIRECTORIES FOR FEATURES, KERNELS,
            %DICTIONARIES AND RESULTS
            
            if ~exist(obj.featPath,'dir')
                mkdir(obj.featPath);
            end
            if ~exist(obj.kernelsPath,'dir')
                mkdir(obj.kernelsPath);
            end
            if ~exist(obj.dictionaryPath,'dir')
                mkdir(obj.dictionaryPath);
            end
            if ~exist(obj.resultsPath,'dir')
                mkdir(obj.resultsPath);
            end
        end
        
        
        function [accuracy,model,dec_values,predict_labels,confusion]=runExperiment(obj,C,i,kernelType,labels,split)
            %RUNS THE EXPERIMENT AND RETURNS: accuracy, model, decision
            %values, predicted labels and confusion matrix
            %
            %PARAMETERS:
            %C is the SVM regularization parameter (higher C means less
            %regularization)
            %
            %i indicates which of the splits is being used
            %
            % kernelType specifies the type of kernel: 1 for linear, 2 for
            % intersection, 3 for chi square, 4 for gaussian
            %
            %labels is the vector with the labels
            %
            %split is matrix containing the training/testing splits
            
            if ~exist('i','var')
                i=1;
            end
            
            if ~exist('kernelType','var')
                kernelType=obj.kernelType;
            end
            
            [Ktr,Kte,labels_tmp]=obj.computeKernels(i,kernelType);
            if ~exist('labels','var')
                labels=labels_tmp;
            end
            if ~exist('split','var')
                split=obj.createTrTeSplit(labels,i);
            end
            
            fprintf('RunningExperiment %d.. \n',i);
            save_file=sprintf('%s/SVM_C%d_%s_PHOW_IS%d_PS%d_SP%d_NTR%d_VF%d_K%d_L%d_PW%d_N%d.mat',obj.resultsPath,C,obj.KERNAMES{kernelType},obj.maxImageSize,obj.patchSize,obj.spacing,obj.nTrain,obj.nVocFeats,obj.K,obj.L,obj.pyramidWeights,i);
            ytr=labels(split(:,i));
            yte=labels(~split(:,i));
            model=PhowExperience.trainOvaSVM(Ktr,ytr',C);
            save(save_file,'model');
            [dec_values,predict_labels,accuracy,confusion]=PhowExperience.testOvaSVM(Kte,yte',model);
            save(save_file,'model','C','accuracy','dec_values','predict_labels','confusion');
        end
        
        
        function [Ktr,Kte,labels]=computeKernels(obj,i,kernelType,phow,labels,split)
            %COMPUTES THE TRAINING Ktr AND TESTING Kte KERNELS
            %
            %PARAMETERS:
            %i indicates which of the splits is being used
            %
            % kernelType specifies the type of kernel: 1 for linear, 2 for
            % intersection, 3 for chi square, 4 for gaussian
            %
            %phow is the matrix containing all the phow features
            %
            %split is matrix containing the training/testing splits
            
            if ~exist('i','var')
                i=1;
            end
            
            if ~exist('kernelType','var')
                kernelType=obj.kernelType;
            end
            
            feat_file=sprintf('%s/%s_PHOW_IS%d_PS%d_SP%d_NTR%d_VF%d_K%d_L%d_PW%d_N%d.mat',obj.kernelsPath,obj.KERNAMES{kernelType},obj.maxImageSize,obj.patchSize,obj.spacing,obj.nTrain,obj.nVocFeats,obj.K,obj.L,obj.pyramidWeights,i);
            try
                fprintf('Trying to load %s Kernels %d..\n',obj.KERNAMES{kernelType},i);
                if obj.doFromScratch
                    throw(MException('i:d','DoFromScratch'));
                end
                load(feat_file);
                Ktr;
                Kte;
                labels;
                fprintf('Loaded %s Kernels %d.\n',obj.KERNAMES{kernelType},i);
            catch
                if ~exist('phow','var')
                    [phow,labels_tmp]=obj.featurePooling(i);
                end
                if ~exist('labels','var')
                    labels=labels_tmp;
                end
                if ~exist('split','var')
                    split=obj.createTrTeSplit(labels,i);
                end
                fprintf('Computing %s Kernels %d.. \n',obj.KERNAMES{kernelType},i);
                if kernelType==1
                    Ktr=PhowExperience.computeLinearKernel(phow(:,split(:,i)),phow(:,split(:,i)));
                    Kte=PhowExperience.computeLinearKernel(phow(:,~split(:,i)),phow(:,split(:,i)));
                elseif kernelType==2
                    Ktr=PhowExperience.computeIntersectionKernel(phow(:,split(:,i)),phow(:,split(:,i)));
                    Kte=PhowExperience.computeIntersectionKernel(phow(:,~split(:,i)),phow(:,split(:,i)));
                elseif kernelType==3
                    [Ktr,gamma]=PhowExperience.computeChiSquareKernel(phow(:,split(:,i)),phow(:,split(:,i)));
                    Kte=PhowExperience.computeChiSquareKernel(phow(:,~split(:,i)),phow(:,split(:,i)),gamma);
                elseif kernelType==4
                    [Ktr,gamma]=PhowExperience.computeGaussianKernel(phow(:,split(:,i)),phow(:,split(:,i)));
                    Kte=PhowExperience.computeGaussianKernel(phow(:,~split(:,i)),phow(:,split(:,i)),gamma);
                end
                save(feat_file, 'Ktr','Kte','labels','-v7.3');
            end
        end
        
        
        function [phow,labels]=featurePooling(obj,i,codes,labels)
            %PERFORMS THE SPATIAL PYRAMID POOLING OF THE FEATURES
            %phow=featurePooling(obj,codes,labels,i) returns a matrix of
            %size obj.K*sum(2.^(2*[0:L])) x number_of_features
            %
            %PARAMETERS:
            %i indicates which of the splits is being used
            %
            %codes is the structure containing the Bag of Visual Words
            %(hard) encoding
            %
            %labels is the vector with the labels
            %
            %split is matrix containing the training/testing splits
            
            if ~exist('i','var')
                i=1;
            end
            
            feat_file=sprintf('%s/PHOW_IS%d_PS%d_SP%d_NTR%d_VF%d_K%d_L%d_PW%d_N%d.mat',obj.featPath,obj.maxImageSize,obj.patchSize,obj.spacing,obj.nTrain,obj.nVocFeats,obj.K,obj.L,obj.pyramidWeights,i);
            try
                fprintf('Trying load PHOW %d..\n',i);
                if obj.doFromScratch
                    throw(MException('i:d','DoFromScratch'));
                end
                load(feat_file);
                phow;
                labels;
                fprintf('Loaded PHOW %d.',i);
            catch
                if ~exist('codes','var')
                    [codes,labels_tmp]=obj.featureEncoding(i);
                end
                if ~exist('labels','var')
                    labels=labels_tmp;
                end
                fprintf('Pooling visual words %d.. ',i);
                n_feats=numel(codes);
                p_o=0;
                for i=1:n_feats
                    prog=round(i/n_feats*100);
                    if ~rem(prog,10) && (prog-p_o)
                        fprintf('%d,',prog);
                        p_o=prog;
                    end
                    bw=zeros(codes(i).h,codes(i).w);
                    bw(sub2ind(size(bw),codes(i).pos(2,:),codes(i).pos(1,:)))=codes(i).feats;
                    phow(:,i)=PhowExperience.spatialPyramidPooling(bw,2,obj.L,obj.pyramidWeights,obj.K);
                end
                save(feat_file, 'phow','labels','-v7.3');
            end
            fprintf('\n');
        end
        
        
        function [bow,labels]=featureEncoding(obj,i,sifts,labels,dictionary)
            %PERFORMS THE SPATIAL PYRAMID POOLING OF THE FEATURES
            %bow=featureEncoding(obj,sifts,labels,dictionary,i) returns a
            %structure containing the fields: feats, pos, sq_dist, h, w
            %where: - feats is a matrix of size obj.K x number_of_features
            %       with the visual words encoding of the SIFT features,
            %       - pos is a 2d-matrix containing the position of the
            %       extracted features,
            %       - sq_dist contains the distances from the centers of
            %       the visual words,
            %       - h and w are the height and width of the original
            %       image
            %
            %PARAMETERS:
            %i indicates which of the splits is being used
            %
            %sifts is the structure containing the extracted SIFT features
            %
            %labels is the vector with the labels
            %
            %dictionary is the dictionary that has to be used to perform
            %the encoding
            
            if ~exist('i','var')
                i=1;
            end
            
            feat_file=sprintf('%s/Codes_SIFT_IS%d_PS%d_SP%d_NTR%d_VF%d_K%d_N%d.mat',obj.featPath,obj.maxImageSize,obj.patchSize,obj.spacing,obj.nTrain,obj.nVocFeats,obj.K,i);
            try
                fprintf('Trying to load BoW %d..\n',i);
                if obj.doFromScratch
                    throw(MException('i:d','DoFromScratch'));
                end
                load(feat_file);
                bow;
                labels;
                fprintf('Loaded BoW %d.\n',i);
            catch
                if ~exist('sifts','var')
                    [sifts,labels_tmp]=obj.featureExtraction();
                end
                if ~exist('labels','var')
                    labels=labels_tmp;
                end
                if ~exist('dictionary','var')
                    dictionary=obj.dictionaryLearning(i,sifts,labels);
                end
                fprintf('Encoding features %d.. ',i);
                n_feats=numel(sifts);
                p_o=0;
                for i=1:n_feats
                    prog=round(i/n_feats*100);
                    if ~rem(prog,10) && (prog-p_o)
                        fprintf('%d,',prog);
                        p_o=prog;
                    end
                    bow(i)=PhowExperience.hardCodeFeats(sifts(i),dictionary);
                end
                save(feat_file, 'bow','labels','-v7.3');
            end
            fprintf('\n');
        end
        
        
        function dictionary=dictionaryLearning(obj,i,sifts,labels,split)
            %CONSTRUCTS THE VOCABULARY
            %dictionary=dictionaryLearning(obj,sifts,split,i) returns a
            %matrix of size 128 x obj.K
            %structure containing the fields: feats, pos, sq_dist, h, w
            %
            %PARAMETERS:
            %i indicates which of the splits is being used
            %
            %sifts is the structure containing the extracted SIFT features
            %
            %split is matrix containing the training/testing splits
            
            if ~exist('i','var')
                i=1;
            end
            
            try
                rng('default');
                rng(i);
            catch
                rand('state',i);
            end
            feat_file=sprintf('%s/SIFT_IS%d_PS%d_SP%d_NTR%d_VF%d_N%d',obj.dictionaryPath,obj.maxImageSize,obj.patchSize,obj.spacing,obj.nTrain,obj.nVocFeats,i);
            voc_file=sprintf('%s/Voc_SIFT_IS%d_PS%d_SP%d_NTR%d_VF%d_K%d_N%d',obj.dictionaryPath,obj.maxImageSize,obj.patchSize,obj.spacing,obj.nTrain,obj.nVocFeats,obj.K,i);
            try
                fprintf('Trying to load dictionary..\n');
                if obj.doFromScratch
                    throw(MException('i:d','DoFromScratch'));
                end
                load(voc_file);
                dictionary;
                fprintf('Loaded dictionary.\n');
            catch
                try
                    if obj.doFromScratch
                        throw(MException('i:d','DoFromScratch'));
                    end
                    load(feat_file);
                catch
                    if ~exist('sifts','var')
                        [sifts,labels_tmp]=obj.featureExtraction();
                    end
                    if ~exist('labels','var')
                        labels=labels_tmp;
                    end
                    if ~exist('split','var')
                        split=obj.createTrTeSplit(labels,i);
                    end
                    fprintf('Building dictionary %d..\n',i);
                    %MAKES USE OF THE TRAINING FEATURES TO LEARN THE
                    %DICTIONARY
                    Ftr=sifts(split(:,i));
                    n_tr=numel(Ftr);
                    [X{1:n_tr}] = deal(Ftr.feats);
                    Ftr=cell2mat(X);
                    clear X;
                    %SELECTS A RANDOM SUBSET (WITH n_voc_feats LOCAL DESCRIPTORS)
                    %FROM THE TRAINING FEATURES
                    n_voc_feats=min(obj.nVocFeats,size(Ftr,2));
                    p=randperm(n_voc_feats);
                    %SAVES THE FEATURES THAT ARE TO BE USED TO LEARN THE
                    %DICTIONARY
                    VocFeats=Ftr(:,p);
                    save(feat_file,'VocFeats','-v7.3');
                end
                dictionary = sprintf('Complete the code to run the experiment:\nTip: use vl_kmeans on VocFeats, with obj.K visual words..\nhttp://www.vlfeat.org/mdoc/vl_kmeans.html');
                error(dictionary);
                save(voc_file,'dictionary');
            end
        end
        
        
        function [sifts,labels,classes]=featureExtraction(obj)
            %PERFORMS SIFT FEATURE EXTRACTION
            %[sifts,labels,classes]=featureExtraction(obj) returns an array
            %sift, containing the fields: feats, pos, h, w
            %where: - feats is a matrix of size 128 x number_of_local_features
            %       - pos is a 2d-matrix containing the position of the
            %       extracted features,
            %       - h and w are the height and width of the original
            %       image
            %it returns also the labels and the name of the classes
            
            %EXTRACT THE SIFT FEATURES FROM THE DATASET
            feat_file=sprintf('%s/SIFT_IS%d_PS%d_SP%d.mat',obj.featPath,obj.maxImageSize,obj.patchSize,obj.spacing);
            try
                fprintf('Trying to load SIFT features..\n');
                if obj.doFromScratch
                    throw(MException('i:d','DoFromScratch'));
                end
                load(feat_file);
                sifts;
                labels;
                classes;
                fprintf('Loaded features.\n');
            catch
                fprintf('Extracting SIFT features..\n');
                d=dir(obj.imgPath);
                k=0;
                c=0;
                for i=3:numel(d)
                    if d(i).isdir
                        c=c+1;
                        classes{c}=d(i).name;
                        fprintf('%s: ',classes{c});
                        cat_dir=sprintf('%s/%s',obj.imgPath,classes{c});
                        d2=dir(cat_dir);
                        
                        p_o=0;
                        for j=3:numel(d2)
                            prog=round((j-2)/(numel(d2)-2)*100);
                            if ~rem(prog,10) && (prog-p_o)
                                fprintf('%d,',prog);
                                p_o=prog;
                            end
                            f_name=d2(j).name;
                            if ~strcmp(f_name,'.dircksum')
                                k=k+1;
                                imfile=sprintf('%s/%s',cat_dir,f_name);
                                sifts(k)=PhowExperience.computeSIFT(imfile,obj.maxImageSize,obj.patchSize,obj.spacing);
                                labels(k)=c;
                            end
                        end
                        fprintf('\n');
                    end
                end
                save(feat_file, 'sifts','labels','classes','-v7.3');
            end
        end
        
        
        function split=createTrTeSplit(obj,labels,n_splits)
            %CREATES THE TRAINING/TEST SPLITS
            %[split,splitFile]=createTrTeSplit(obj,labels,n_splits) returns
            %a boolean matrix splits which is true for the training samples
            %and false for the testing ones
            %
            %PARAMETERS
            %
            %labels is the vector with the labels
            %
            %n_splits is the number of splits to be created
            
            classes=unique(labels);
            n_images=numel(labels);
            split=false(n_images,n_splits);
            sF=sprintf('%s_NTR%d.mat',obj.splitFile,obj.nTrain);
            try
                rng('default');
            catch
            end
            for i=1:n_splits
                try
                    fprintf('Trying to load training/testing splits..\n');
                    if obj.doFromScratch
                        throw(MException('i:d','DoFromScratch'));
                    end
                    load(sF);
                    split(:,i);
                    if sum(split(:,i)) ~= numel(classes)*obj.nTrain
                        error('split: %d need to be redone',i);
                    end
                    fprintf('Loaded training/testing splits.\n');
                catch
                    fprintf('Creating training/testing splits..\n');
                    try
                        rng(i-1);
                    catch
                        rand('state',i-1);
                    end
                    for c=classes
                        l_cla=find(labels==c);
                        r=randperm(numel(l_cla));
                        split(l_cla(r(1:obj.nTrain)),i)=true;
                    end
                    save(sF,'split');
                end
            end
        end
    end
    
    methods (Static)
        
        function sift=computeSIFT(imfile,max_image_size,patch_size,spacing)
            %COMPUTES A SIFT DESCRIPTOR FOR AN IMAGE
            %sift=computeSIFT(imfile,max_image_size,patch_size,spacing)
            %returns a structure sift, containing the fields: feats, pos, h, w
            %where: - feats is a matrix of size 128 x number_of_local_features
            %       - pos is a 2d-matrix containing the position of the
            %       extracted features,
            %       - h and w are the height and width of the original
            %       image
            %
            %PARAMETES
            % max_image_size is the maximum size allowed for an image
            %
            % patch_size is the size of the side of a square SIFT patch (with 4x4 bins)
            %
            % spacing is the spacing between SIFT points (for dense sampling)
            
            %READS AN IMAGE
            I=imread(imfile);
            if ndims(I) == 3
                I = im2single(rgb2gray(I));
            else
                I = im2single(I);
            end
            if max(size(I)) > max_image_size
                % if the longest edge of the image is longer than max_image_size pixels, scale it to
                % max_image_size pixels.
                I = imresize(I,max_image_size/max(size(I)));
            end
            [h,w]=size(I);
            bin_size=sprintf('Complete the code to run the experiment:\nTip: bin_size is a function of patch_size..\nhttp://www.vlfeat.org/mdoc/vl_dsift.html');
            error(bin_size);
            [pos, feats] = vl_dsift(I,'step',spacing,'size',bin_size);
            sift.pos=pos;
            sift.feats=feats;
            sift.h=h;
            sift.w=w;
        end
        
        function bow=hardCodeFeats(sift,dictionary)
            %GIVEN AN ARRAY OF SIFT FEATURES COMPUTES FOR EACH OF THEM, THE
            %CLOSEST NEIGHBOR AND RETURNS THE CORRESPONDING ID
            
            codes=sprintf('Complete the code to run the experiment\nTip: avoid using for-loops..');
            error(codes);
            
            bow.feats=codes;
            bow.sq_dist=sq_dist;
            bow.pos=sift.pos;
            bow.h=sift.h;
            bow.w=sift.w;
        end
        
        function phow=spatialPyramidPooling(wf,b,L,pyramidWeights,K)
            %phow=spatialPyramidPooling(wf,b,L,pyramidWeights,K)
            %GIVEN A BOW ENCODED IMAGE wf, PERFORMS THE POOLING OF THE
            %FEATURES, IN A SPATIAL PYRAMID FASHION, WHERE:
            %
            %PARAMETERS
            % sum(b.^(2*[0:L])) IS THE NUMBER OF HISTOGRAMS CONCATENATED
            %
            % pyramidWeights is a boolean value stating whether the spatial
            % pyramid must be normalized using the weights, as in the paper
            %
            % K is the number of visual words in the dictionary
            
            %if true visualizes the patches over which the histograms are
            %being collected
            visualize=false;
            
            xs=size(wf,2);
            ys=size(wf,1);
            
            if xs< b^L
                error('Image width is %d, while level %d requires %d pixels',xs,L,b^L);
            end
            if ys< b^L
                error('Image height is %d, while level %d requires %d pixels',ys,L,b^L);
            end
            
            phow = [];
            l_list=0:L;
            
            for l=l_list
                if pyramidWeights
                    if l==0
                        w=2^(-L);
                    else
                        w=2^(-L+l-1);
                    end
                end
                
                n_splits=b^l;
                
                x_wd=xs/(n_splits);
                y_wd=ys/(n_splits);
                tmp=round((1:n_splits).*x_wd);
                x_int=[0 tmp];
                tmp=round((1:n_splits).*y_wd);
                y_int=[0 tmp];
                
                cx=0;
                for x=x_int(1:end-1)
                    cx=cx+1;
                    cy=0;
                    for y=y_int(1:end-1)
                        cy=cy+1;
                        y_min=y+1;
                        y_max=y_int(cy+1);
                        x_min=x+1;
                        x_max=x_int(cx+1);
                        
                        if visualize
                            %visualizes the patches over which the histograms are
                            %being collected
                            vis=zeros(size(wf));
                            vis(y_min:y_max,x_min:x_max)=255;
                            colormap(gray);
                            image(vis);
                        end
                        
                        wf_cella = wf(y_min:y_max,x_min:x_max,:);
                        h=zeros(K,1);
                        for bn=1:K
                            ind3 = wf_cella==bn;
                            h(bn)=sum(sum(ind3));
                        end
                        if pyramidWeights
                            phow = [phow;h*w/sum(h)];
                        else
                            phow = [phow;h];
                        end
                    end
                end
            end
            phow=phow/sum(phow);
        end
        
        function K=computeLinearKernel(X1,X2)
            %COMPUTES THE LINEAR KERNEL BETWEEN TWO FEATURES:
            %X1 AND X2
            
            K=X1'*X2;
        end
        
        function K=computeIntersectionKernel(X1,X2)
            %COMPUTES THE INTERSECTION KERNEL BETWEEN TWO FEATURES:
            %X1 AND X2
            
            K=zeros(size(X1,2),size(X2,2));
            for i=1:size(X1,2)
                h1=X1(:,i);
                for j=1:size(X2,2)
                    h2=X2(:,j);
                    tmp = sprintf('Complete the code to run the experiment:\nTip: equation 1) in Lazebnik''s paper');
                    error(tmp);
                    K(i,j)=tmp;
                end
            end
        end
        
        function [K,gamma]=computeChiSquareKernel(X1,X2,gamma)
            %COMPUTES THE CHI-SQUARE KERNEL BETWEEN TWO FEATURES:
            %X1 AND X2
            
            %%% COMPUTES THE DISTANCE MATRIX
            D=zeros(size(X1,2),size(X2,2));
            for i=1:size(X1,2)
                h1=X1(:,i);
                for j=1:size(X2,2)
                    h2=X2(:,j);
                    D(i,j)=sum(((h1-h2).^2)./(h1+h2+eps));
                end
            end
            clear X1 X2;
            
            %%% TUNES THE GAMMA
            if ~exist('gamma','var')
                gamma=PhowExperience.tuneGamma(D);
            end
            
            %%% COMPUTES THE KERNEL
            K=exp(-gamma*D);
        end
        
        function [K,gamma]=computeGaussianKernel(X1,X2,gamma)
            %COMPUTES THE GAUSSIAN KERNEL BETWEEN TWO FEATURES:
            %X1 AND X2
            
            nsamp1=size(X1,2);
            nsamp2=size(X2,2);
            
            %%% COMPUTES REQUIRED QUANTITIES
            normX = full(sum(X1.^2,1));
            normY = full(sum(X2.^2,1));
            dotprod=full(X1'*X2);
            clear X1 X2;
            
            %%% COMPUTES THE DISTANCE MATRIX
            D=repmat(normX' ,1,nsamp2) + repmat(normY,nsamp1,1) -  2*dotprod;
            
            %%% TUNES THE GAMMA
            if ~exist('gamma','var')
                gamma=PhowExperience.tuneGamma(D);
            end
            
            %%% COMPUTES THE KERNEL
            K=exp(-gamma/2*D);
        end
        
        function gamma=tuneGamma(D)
            %PROVIDED A DISTANCE MATRIX D, IT COMPUTES THE INVERSE OF THE
            %AVERAGE DISTANCE
            n_samp=size(D,1);
            s=sum(sum(D));
            gamma=(n_samp^2)/s;
        end
        
        function model=trainOvaSVM(kernel,labels,C)
            %TRAINS A ONE VS ALL CLASSIFIER
            
            l=unique(labels);
            fprintf('Training using One-VS-All scheme\n');
            for c=1:max(l)
                fprintf('\t class %d\n',c);
                y_tr=(labels==c)*2-1;
                model.model(c)=PhowExperience.trainSVM(kernel,y_tr,C);
            end
        end
        
        function [dec_values,predict_labels,accuracy,confusion]=testOvaSVM(kernel,labels,model)
            %TESTS THE ONE VS ALL CLASSIFIER AND RETURNS: decision
            %values, predicted labels, accuracy and confusion matrix
            
            nBest=1;
            nCla=numel(model.model);
            nSamp=size(kernel,1);
            fprintf('Testing using One-VS-All scheme\n');
            dec_values=zeros(nSamp,nCla);
            for c=1:nCla
                fprintf('\t class %d\n',c);
                y_te=(labels==c)*2-1;
                tmp_dec_values= PhowExperience.testSVM(kernel, y_te, model.model(c));
                if isempty(tmp_dec_values)
                    tmp_dec_values=-inf(size(tmp_dec_values,1),1);
                end
                try
                    dec_values(:,c)=tmp_dec_values;
                catch
                    error('cannot process decision values\n');
                end
            end
            [slack,predict_labels_M]=sort(dec_values,2,'descend');
            
            tmp0=predict_labels_M(:,1:nBest);
            tmp=tmp0==repmat(labels,1,nBest);
            [slack,i]=max(tmp,[],2);
            predict_labels=diag((tmp0(1:numel(i),i)));
            
            accuracy=sum(predict_labels==labels)/numel(labels)*100;
            confusion=PhowExperience.confmat(uint8(labels),uint8(predict_labels));
        end
        
        
        function model=trainSVM(kernel,labels,C)
            %TRAINS AN SVM AND RETURNS THE TRAINED model
            %
            %PARAMETERS:
            %kernel is the training kernel matrix
            %            %
            %labels is the vector with the labels
            %
            %C is the SVM regularization parameter (higher C means less
            %regularization)
            
            n =numel(labels);
            param='-t 4';
            if exist('C','var')
                param=sprintf('%s -c %d',param,C);
            end
            model=svmtrain(double(labels), double([(1:n)', kernel]),param );
            clear kernel labels;
        end
        
        function [dec_values,predict_labels,accuracy]=testSVM(kernel,labels,model)
            %TESTS THE SVM CLASSIFIER AND RETURNS: decision
            %values, predicted labels and accuracy
            %
            %PARAMETERS:
            %kernel is the testing kernel matrix
            %
            %labels is the vector with the labels
            %
            %model is a trained SVM model
            
            n =numel(labels);
            [predict_labels, accuracy, dec_values] = svmpredict(double(labels), double([(1:n)',kernel]), model);
            if numel(model.Label)==2
                dec_values=dec_values*sign(model.Label(1));
            end
            accuracy=accuracy(1);
        end
        
        function c = confmat(label,prediction)
            %COMPUTES THE CONFUSION MATRIX
            %label is the deired output
            %prediction is the predicted output
            
            miny = min(label);
            maxy = max(label);
            c = zeros(maxy-miny);
            for i = miny:maxy
                index = label == i;
                for j = miny:maxy
                    z = prediction(index);
                    c(i-miny+1,j-miny+1) = length(find(z == j));
                end
            end
        end
    end
end
