function dictPCA = SubspacePCA(dictionary, numPcaDimensions, descrs, index)
% DictPCA = SubspacePCA(dictionary, numPcaDimensions, numWords, descrs, index)
% build subspaces of each atom in the dictionary.
% 
%            dictPCA: an ensemble of subspaces (e.g. bases and features).
%         dictionary: dictionary obtained by k-means. DxM where D is the dimension 
%                     of descriptors and M is dictionary size .
%   numPcaDimensions: subspace dimension.
%             descrs: sample image descriptors. DxN (N is #descriptors sampled from training set).
%              index: the nearest index of atom for each descriptor , Nx1.
% 
% -------------------------------------------------------------------------
% For theoretical and technical details, please refer to the following paper:
%
% Peihua Li, Xiaoxiao Lu, Qilong Wang, "From Dictionary of Visual 
% Words to Subspaces: Locality-constrained Affine Subspace Coding", CVPR 2015.
%
% Please cite the paper above if you use the code.
% For questions, please contact: Xiaoxiao Lu (Email: shaw@mail.dlut.edu.cn)
%                                Qilong Wang (Email: wangqilong.415@163.com)
%                                Peihua Li   (Email: peihuali at dlut dot edu dot cn)
% -------------------------------------------------------------------------

numWords = size(dictionary, 2);
dictPCA.proj = cell(numWords,1);
dictPCA.err= cell(numWords,1);
dictPCA.log_coef = cell(numWords,1);
dictPCA.pcaNum = zeros(numWords, 1);

for i = 1:numWords
    sift_set = descrs(:, index==i);
    sift_set = bsxfun(@minus, sift_set, dictionary(:,i)); 
    cov_matrix = cov(sift_set');
    [vectors,values] = eig(cov_matrix);
    [max_value,max_idx] = sort(diag(values),'descend');
    vectors = vectors(:,max_idx);
    values = diag(max_value);

    max_vectors = vectors(:, 1:numPcaDimensions); 
    max_values = values(1:numPcaDimensions, 1:numPcaDimensions);
    
    dictPCA.pcaNum(i) = numPcaDimensions;
    dictPCA.proj{i} = max_values^(-1/2) * max_vectors';                % whitening
    dictPCA.log_coef{i} = log(2*pi)*ones(numPcaDimensions,1)...
                          +log(max_value(1:numPcaDimensions));         % for d_p
    dictPCA.err{i} = eye(size(descrs,1)) - max_vectors * max_vectors'; % for d_r

end