function feaSet = extract_sift(I, pars)

% parameter parsing
gridSpacing = pars.gridSpacing;
patchSize = pars.patchSize;
numOBins = pars.numOBins;
numSBins = pars.numSBins;
maxImSize = pars.maxImSize;

% read and rescale the image
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 = [];

[sift_arr, grid_x, grid_y] = dense_sift(I, numOBins, numSBins, patchSize, gridSpacing);
feaSet.feaArr = sift_arr;
feaSet.x = grid_x(:);
feaSet.y = grid_y(:);
% % 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];

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