function vecs = siftVects(img, bbox, varargin)

% img - single-channel, double. it is the user's responsibility to
%       blur/rescale it if multi-scale analysis is being performed
% bbox - patches to find sift vector for
%        ex. bbox(i,:) = [ left right top bottom ] = [xmin xmax ymin ymax] 
%            inclusive dimensions (so the width is xmax-xmin+1)
% varargin - optional arguments
%
% vecs - output, MxN, where there were M bboxes, and N is SIFT vector dimension
%        currently output is of type double, normalized to length 1,
%        however, if desired, you can introduce *more* nonlinearities by
%        scaling to 0..255 and converting to UINT8
%
% daniel eaton, 2005, danieljameseaton@gmail.com
% notes: this should be programmed in c, because histogram indexing is a
% nightmare in matlab

[nSpatialBins nOriBins] = process_options(varargin, 'nSpatialBins', 3, ...
	'nOriBins', 4);

% compute the edge magnitude and orientation for each pixel
xPartial = conv2(img, [1 0 -1], 'same');
yPartial = conv2(img, [1 0 -1]', 'same');
ori = atan2(yPartial, xPartial);
mag = sqrt(xPartial.*xPartial + yPartial.*yPartial);

nBbox = size(bbox,1);
siftDimension = nSpatialBins^2 * nOriBins;
vecs = zeros(nBbox, siftDimension);

oriSpacing = 2*pi/nOriBins;
oriBins = [-pi:oriSpacing:(pi-oriSpacing)];

for bi=1:nBbox
	bb = bbox(bi,:);
	
	bbw = bb(2)-bb(1)+1;
	bbh = bb(4)-bb(3)+1;

	% discretize the box into nSpatialBins x nSpatialBins bins (each bin in
	% turn is split into nOriBins bins for orientation)
	xSpacing = bbw/(nSpatialBins);
	ySpacing = bbh/(nSpatialBins);
	xBins = xSpacing/2:xSpacing:bbw;
	yBins = ySpacing/2:ySpacing:bbh;

	oriEdgeHist = zeros(nSpatialBins, nSpatialBins, nOriBins);
	
	% compute the weights for histogram incrementation
	% point contributes to its 8 nearest neighbours
	for ci=1:bbw
		% compute which x bins to insert into (xw is 1xnSpatialBins, at most 2 elements
		% will be nonzero, and will be the weights for the gradient magnitude
		% 'vote' into the hist)
		% repmat xw s.t. it is nSpatialBins x nSpatialBins x nOriBins in order
		% to easily increment the histogram
		xw = repmat(max(1 - (abs(xBins - ci + 0.5)/xSpacing), 0), [nSpatialBins 1 nOriBins]);
		for ri=1:bbh
			% similarly for the y bins
			yw = repmat( max(1 - (abs(yBins - ri + 0.5)/ySpacing), 0)',  [1 nSpatialBins nOriBins]);
			
			% similarly for the orientation bins (now weight is determined by
			% distance to bin center modulo 2pi)
      ow = repmat( reshape(max(1 - abs( mod( ori( bb(3)+ri-1, bb(1)+ci-1 ) - oriBins + pi, 2*pi ) - pi )/oriSpacing, 0), ...
				[1 1 nOriBins]), [nSpatialBins nSpatialBins 1]);

			oriEdgeHist = oriEdgeHist + mag( bb(3)+ri-1, bb(1)+ci-1 )*xw.*yw.*ow;
		end
	end
	
	siftVector = reshape(oriEdgeHist, [1 siftDimension]);
	% normalize 
	siftVector = siftVector/sqrt(siftVector*siftVector');
	% threshold like lowe, then renormalize (apparently, to reduce the influence of large gradient magnitudes)
  siftVector = min(siftVector, 0.2);
	siftVector = siftVector/sqrt(siftVector*siftVector');
	
	vecs(bi, :) = siftVector;
end