function siftFeatureSet = sift_descriptors_computing(imageSet, sift_fea_dir, pars)
%==========================================================================
% usage: calculate the sift descriptors given the image directory
%
% inputs
% imageSet      -image database
% sift_fea_dir  -feature database root path
% gridSpacing   -spacing for sampling dense descriptors
% patchSizes    -patch size for extracting sift feature
% maxImSize     -maximum size of the input image
% nrml_threshold    -low contrast normalization threshold
%
% outputs
% siftFeatureSet      -the set for the calculated sift features
%
% Lazebnik's SIFT code is used.
%
% modified by Jimei Yang @Sony, May 2011
%==========================================================================

gridSpacing = pars.gridSpacing;
patchSizes = pars.patchSizes;
maxImSize = pars.maxImSize;
numOBins = pars.numOBins;
numSBins = pars.numSBins;


siftFeatureSet = imageSet;
siftFeatureSet.path = {};
if ~isdir(sift_fea_dir)
    mkdir(sift_fea_dir);
end

disp('Extracting SIFT descriptors...');

for ii = 1:imageSet.imnum,
    tic;
    impath = imageSet.imgpath{ii};
    [pathstr, name, ext] = fileparts(impath);
    fprintf('Processing the image %s...', name);
    
    fpath = fullfile(sift_fea_dir, [name, '.mat']);
    siftFeatureSet.path = [siftFeatureSet.path; fpath];
    
    if exist(fpath,'file')
        fprintf('Skipped.\n');
        continue;
    end
    
    % read and rescale the image
    I = imread(impath);
    if ndims(I) == 3,
        I = im2double(rgb2gray(I));
    else
        I = im2double(I);
    end

    [im_h, im_w] = size(I);

    if max(im_h, im_w) > maxImSize,
        I = imresize(I, maxImSize/max(im_h, im_w), 'bicubic');
        [im_h, im_w] = size(I);
    end;

    feaSet.feaArr = [];
    feaSet.x = [];
    feaSet.y = [];
    for p = 1 : length(patchSizes)
        patchSize = patchSizes(p);
        
        % make grid sampling SIFT descriptors
        remX = mod(im_w-patchSize,gridSpacing);
        offsetX = floor(remX/2)+1;
        remY = mod(im_h-patchSize,gridSpacing);
        offsetY = floor(remY/2)+1;

        [gridX,gridY] = meshgrid(offsetX:gridSpacing:im_w-patchSize+1, offsetY:gridSpacing:im_h-patchSize+1);

        % find SIFT descriptors
        siftArr = sp_find_sift_grid(I, gridX, gridY, patchSize, numOBins, numSBins, 0.8);
        siftArr = sp_normalize_sift(siftArr, 1);

        feaSet.feaArr = [feaSet.feaArr, siftArr'];
        x = gridX(:) + patchSize/2 - 0.5;
        y = gridY(:) + patchSize/2 - 0.5;
        feaSet.x = [feaSet.x;x];
        feaSet.y = [feaSet.y;y];

    end

    feaSet.width = im_w;
    feaSet.height = im_h;


    save(fpath, 'feaSet');

    fprintf(' in %f seconds.\n', toc);
end;

disp('Done!');
