function [ sk ] = spirit_kernel(mc, Nk)
%
%  sk = spirit_kernel(mc, Nk)
%
%  Compute the SPIRiT kernel for a calibration region mc
%
%  mc -- calibration region
%  Nk -- kernel size, should be odd (like 5)
%
%  sk -- SPIRiT kernel 
%        sk is [Nk Nk Nc Nc] in size
%

[Nx Ny Nc] = size(mc);

% set up a matrix with all of the calibration data 
nd = [];
for kk=1:Nc,
    nd = [nd im2col(mc(:,:,kk),[Nk Nk], 'sliding').'];
end;

% regularization parameter, from Miki Lustig's code
lambda = norm(nd'*nd,'fro')/size(nd,2)*0.01;

% initialize kernel
sk = zeros(Nk,Nk,Nc,Nc);

% offset of the synthesized sample for first coil
ndx1 = floor(Nk*Nk/2)+1;

% calculate weights for each channel
for k=1:Nc,
    
   % index of synthesized sample for coil k
   ndx = ndx1 + (k-1)*Nk*Nk;
   
   % find weights by removing corresponding column and fitting to it
   A = [nd(:,1:ndx-1) nd(:,ndx+1:end)];
   AtA = A'*A;
   temp = inv(AtA + eye(size(AtA))*lambda)*A'*nd(:,ndx);
   
   % add a zero placeholder at synthesized sample so we can implement as
   % convolution later
   temp = [temp(1:ndx-1).' 0 temp(ndx:end).'];
   
   % save kernel
   sk(:,:,:,k)  = reshape(temp,Nk,Nk,Nc);
   
end;

end