function [d, kmpd]= igridkb(m,k,osf,kw)

% function d = igridkb(m,k,osf,kw)
%
%     m -- image domain data, square
%     k -- k-trajectory, scaled -0.5 to 0.5
%     osf -- oversampling factor (usually between 1 and 2)
%     kw -- full kernel width in 1x samples (usually 3 to 7)
%
%     d -- k-space data
%
%  Uses optimum Kaiser-Bessel window for a given
%    oversampling factor and kernel size
%

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

% convert to single column
k = k(:);

[n n2] = size(m);
nosf = n*osf;
mp = zeros(nosf,nosf);
mp((nosf-n)/2+[1:n],(nosf-n)/2+[1:n]) = m;
dapm = kb_deap(osf,kw,n);
mpd = mp./dapm;

kmpd = fft2c(mpd);

% compute kernel
kosf = 64;
kwidth = osf*kw/2;
om = [0:kosf*kwidth]/(kosf*kwidth);
p = besseli(0,(osf-1/2)*pi*kw*sqrt(1-om.*om));
p = p./p(1);
%last sample is zero so we can use max() below for samples bigger than kwidth
p(kosf*kwidth+1) = 0;

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

ld = length(k);
a = sparse(ld,nosf*nosf);

% loop over samples in kernel
for lx = -osf:osf,
  for ly = -osf:osf,

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

    % seperable kernel value
    kkx = min(round(kosf*abs(nx-nxt)+1), kosf*kwidth+1);
    kwx = p(kkx);
    kky = min(round(kosf*abs(ny-nyt)+1), 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);

    a = a+sparse([1:ld],nxt+nosf*(nyt-1),kwx.*kwy,ld,nosf*nosf);

  end;
end;

d = a*kmpd(:);
