function [m, p] = gridkb(d,k,w,n,osf,wg,opt)

% function m = gridkb(d,k,w,n,osf,kw,opt)
%
%     d -- k-space data
%     k -- k-trajectory, scaled -0.5 to 0.5
%     w -- k-space weighting
%     n -- image size (m will be osf*n X osf*n)
%     osf -- oversampling factor (usually between 1 and 2)
%     wg -- full kernel width in oversampled grid samples (usually 3 to 7)
%     opt -- 'k-space', 'image', defaults it 'image' if not specified
%
%     m -- gridded k-space data
%     p -- gridding kernel, optional
%
%  Uses optimum Kaiser-Bessel window for a given
%    oversampling factor and kernel size
%  Now uses Phil's numbers
%


%  Written by John Pauly, 2003, 2005, 2007, 2011
%  (c)Board of Trustees, Leland Stanford Junior University

if nargin < 7,
    opt = 'image';
end
    
% convert to single column
d = d(:);
k = k(:);
w = w(:);

% width of the kernel on the original grid
kw = wg/osf;

% preweight
dw = d.*w;

% compute kernel, assume e1 is 0.001, assuming nearest neighbor
kosf = floor(0.91/(osf*1e-3));

% half width in oversampled grid units
kwidth = osf*kw/2;

% beta from the Beatty paper
beta = pi*sqrt((kw*(osf-0.5)).^2-0.8);

% compute kernel
om = [0:kosf*kwidth]/(kosf*kwidth);
p = besseli(0,beta*sqrt(1-om.*om));
p = p./p(1);
% last sample is zero so we can use min() below for samples bigger than kwidth
p(end) = 0;

% convert k-space samples to matrix indices
nx = (n*osf/2+1) + osf*n*real(k);
ny = (n*osf/2+1) + osf*n*imag(k);

m = zeros(osf*n,osf*n);

% loop over samples in kernel at grid spacing
for lx = -kwidth:kwidth,
  for ly = -kwidth:kwidth,

    % find nearest samples
    nxt = round(nx+lx);
    nyt = round(ny+ly);

    % seperable kernel value
    kkx = min(round(kosf*abs(nx-nxt)+1), floor(kosf*kwidth)+1);
    kwx = p(kkx);
    kky = min(round(kosf*abs(ny-nyt)+1), floor(kosf*kwidth)+1);
    kwy = p(kky);

    % if data falls outside matrix, put it at the edge, zero out below
    nxt = max(nxt,1); nxt = min(nxt,osf*n);
    nyt = max(nyt,1); nyt = min(nyt,osf*n);

    % accumulate gridded data
    m = m+sparse(nxt,nyt,dw.*kwx'.*kwy',osf*n,osf*n);
  end;
end;

% zero out data at edges, which is probably due to data outside mtx
m(:,1) = 0; m(:,osf*n) = 0;
m(1,:) = 0; m(osf*n,:) = 0;

% stop here, if we just want the k-space data
if strcmp(opt,'k-space') return; end;

im = fftshift(fft2(fftshift(m)));


% compute deappodization function
x = [-osf*n/2:osf*n/2-1]/(n);
sqa = sqrt(pi*pi*kw*kw*x.*x-beta*beta);
dax = sin(sqa)./(sqa);
% normalize by DC value
dax = dax/dax(osf*n/2);
% make it a 2D array
da = dax'*dax;

% deappodize
im = im./da;

%return the result
m = im;