function [data_full,cov_inds] = GRAPPA_recon(sz_full,data_redu,Kernel,ACS,NFirsty,NFirstz,ACSFirstx,ACSFirsty,ACSFirstz)
%
% function [data_full,cov_inds] = GRAPPA_recon(sz_full,data_redu,Kernel,ACS,NFirsty,NFirstz,ACSFirstx,ACSFirsty,ACSFirstz)
%
% Perform GRAPPA reconstruction using calibrated Kernel (see
% GRAPPA_construct_kernel.m) on reduced-FOV data. The first coordinate of
% data_redu corresponds to the coordinate (1,NFirsty,NFirstz,1) of
% full-FOV k-space data_full. Note that the output k-space needs to be
% combined using a post-processing method to form a single composite image.
%
% Inputs:
%   sz_full - four-dimensional size of full k-space
%   data_redu - 4-dimensional reduced-FOV k-space (4th dim. is coil
%       channels)
%   Kernel - calibrated kernel structure
%   ACS - ACS lines (can leave empty if we don't want to preserve them in
%       the output)
%   NFirsty, NFirstz - y and z coordinates of first acquired point in
%       k-space (can be empty to use default-centered coordinates)
%   ACSFirstx, ACSFirsty, ACSFirstz - full k-space x,y,z coordinates of
%       first point in ACS (can leave empty if we don't want to preserve
%       ACS in the output, or if they're default-centered)
%
% Outputs:
%   data_full - 4-dimensional reconstructed full-FOV k-space
%   cov_inds - for each k-space point, the index into the Lambda_G matrix
%       from calibration that represents noise amplification at that point
%
% Copyright 2012-10-09, Daniel Weller, University of Michigan

if isequal(sz_full,'test'), GRAPPA_recon_test(); return; end

trgy = 0:Kernel.Ry-1;
trgz = 0:Kernel.Rz-1;
[trgsy,trgsz] = ndgrid(trgy,trgz);
trgsy = trgsy(2:end);
trgsz = trgsz(2:end);

NFullx = sz_full(1);
NFully = sz_full(2);
NFullz = sz_full(3);
NCha = sz_full(4);
ind_centerx = floor(NFullx/2)+1;
ind_centery = floor(NFully/2)+1;
ind_centerz = floor(NFullz/2)+1;
if isempty(NFirsty), NFirsty = mod(ind_centery-1,Kernel.Ry)+1; end
if isempty(NFirstz), NFirstz = mod(ind_centerz-1,Kernel.Rz)+1; end
data_full = complex(zeros([NFullx,NFully,NFullz,NCha],class(data_redu)));
cov_inds = zeros([NFullx,NFully,NFullz],'uint8');

% do linear convolution with kernels
[NRedux,NReduy,NReduz,NCha] = size(data_redu);
Npadbeforey = ceil((NFirsty-1)/Kernel.Ry);
Npadbeforez = ceil((NFirstz-1)/Kernel.Rz);
Npadaftery = floor((NFully - NFirsty)/Kernel.Ry) - (NReduy-1);
Npadafterz = floor((NFullz - NFirstz)/Kernel.Rz) - (NReduz-1);
Npady = NReduy + Npadbeforey + Npadaftery;
Npadz = NReduz + Npadbeforez + Npadafterz;
data_pad = complex(zeros([NRedux,Npady,Npadz,NCha],class(data_redu)));
data_pad(:,Npadbeforey+1:end-Npadaftery,Npadbeforez+1:end-Npadafterz,:) = data_redu;

for ir = 1:numel(trgsy)
    kernels_r = Kernel.kernels(:,:,:,:,ir,:);
    data_indsy = NFirsty+trgsy(ir)+(-Npadbeforey:NReduy-1+Npadaftery)*Kernel.Ry;
    data_indsz = NFirstz+trgsz(ir)+(-Npadbeforez:NReduz-1+Npadafterz)*Kernel.Rz;
    out_validsy = data_indsy >= 1 & data_indsy <= NFully;
    out_validsz = data_indsz >= 1 & data_indsz <= NFullz;
    data_indsy = data_indsy(out_validsy);
    data_indsz = data_indsz(out_validsz);
    out_indsy = 1:Npady; out_indsz = 1:Npadz;
    cov_inds(:,data_indsy,data_indsz) = ir+1;
    for iq = 1:NCha
        kernels_q = kernels_r(:,:,:,:,:,iq);
        kernels_q = kernels_q(end:-1:1,end:-1:1,end:-1:1,:);
        out_conv = 0;
        for ip = 1:NCha
            out_conv = out_conv + convn(data_pad(:,:,:,ip),kernels_q(:,:,:,ip),'full');
        end
        out_conv = out_conv(Kernel.offx+(1:NRedux),Kernel.offy+(1:Npady),Kernel.offz+(1:Npadz)); % recenter output based on kernel
        data_full(:,data_indsy,data_indsz,iq) = out_conv(:,out_indsy(out_validsy),out_indsz(out_validsz));
    end
end

% fill acquired data
data_full(:,NFirsty+(0:NReduy-1)*Kernel.Ry,NFirstz+(0:NReduz-1)*Kernel.Rz,:) = data_redu;
cov_inds(:,NFirsty+(0:NReduy-1)*Kernel.Ry,NFirstz+(0:NReduz-1)*Kernel.Rz) = 1;

% fill ACS if provided
if ~isempty(ACS)
    [NACSx,NACSy,NACSz,NCha] = size(ACS);
    if isempty(ACSFirstx), ACSFirstx = ind_centerx - floor(NACSx/2); end
    if isempty(ACSFirsty), ACSFirsty = ind_centery - floor(NACSy/2); end
    if isempty(ACSFirstz), ACSFirstz = ind_centerz - floor(NACSz/2); end
    data_full(ACSFirstx:ACSFirstx+NACSx-1,ACSFirsty:ACSFirsty+NACSy-1,ACSFirstz:ACSFirstz+NACSz-1,:) = ACS;
    cov_inds(ACSFirstx:ACSFirstx+NACSx-1,ACSFirsty:ACSFirsty+NACSy-1,ACSFirstz:ACSFirstz+NACSz-1) = 1;
end

end

%
% function GRAPPA_recon_test()
%
% Perform test on 3-D brain phantom from BrainWeb and simulated coil array
function GRAPPA_recon_test()

% load 12-channel brain data
load('/y/wellerds/MRI_data/BrainWeb/T1w_sim12ch.mat');

figure; slicer3(im_ref,'clim',[0,1]);
set(gcf,'Name','Reference Image'); drawnow;

% make multi-channel k-space
[Nx,Ny,Nz,NCha] = size(smap);
ref = im_ref;
ref = fft3c(ref(:,:,:,ones(1,NCha)).*smap);
clear smap;

% create noise
SNR = 40; % dB
noisevar = Nx*Ny*Nz*max(im_ref(:).*smap_norm(:)).^2.*10^(-SNR/10)/NCha; % in multi-channel k-space
Lambda = noisevar.*noisecov;
noise = sqrtm(Lambda./2)*complex(randn([NCha,Nx*Ny*Nz],class(im_ref)),randn([NCha,Nx*Ny*Nz],class(im_ref)));
noise = reshape(noise.',[Nx,Ny,Nz,NCha]);
ref = ref + noise;
clear noise;
im_full = sos_combine(ifft3c(ref),noisecov,smap_norm);
figure; slicer3(abs(im_full-im_ref));
set(gcf,'Name',sprintf('Additive Noise on Full Magnitude Image (SNR: %g dB)',SNR)); drawnow;

% sample k-space
Ry = 3;
Rz = 2;
NACSx = Nx; NACSy = 20; NACSz = 20;
NFirsty = mod(floor(Ny/2),Ry)+1;
NFirstz = mod(floor(Nz/2),Rz)+1;
ACSFirstx = floor(Nx/2) + 1 - floor(NACSx/2);
ACSFirsty = floor(Ny/2) + 1 - floor(NACSy/2);
ACSFirstz = floor(Nz/2) + 1 - floor(NACSz/2);
data_redu = ref(:,NFirsty:Ry:end,NFirstz:Rz:end,:);
ACS = ref(ACSFirstx:ACSFirstx+NACSx-1,ACSFirsty:ACSFirsty+NACSy-1,ACSFirstz:ACSFirstz+NACSz-1,:);

% perform GRAPPA calibration and reconstruction
Bx = 3; By = 4; Bz = 3;
Kernel = GRAPPA_construct_kernel(ACS,Ry,Rz,Bx,By,Bz,Lambda);
data_full = GRAPPA_recon([Nx,Ny,Nz,NCha],data_redu,Kernel,ACS,NFirsty,NFirstz,ACSFirstx,ACSFirsty,ACSFirstz);

% combine coils and display output
im_recon = sos_combine(ifft3c(data_full),noisecov,smap_norm);
figure; slicer3(im_recon,'clim',[0,1]);
set(gcf,'Name','GRAPPA-Reconstructed Image'); drawnow;

figure; slicer3(abs(im_recon-im_ref));
set(gcf,'Name','GRAPPA-Reconstructed Magnitude Noise Image'); drawnow;

end

