explore_recon_dataset.m

function reconImages = explore_recon_dataset(varargin)
% EXPLORE_RECON_DATASET Explore working with an existing ISMRMRD data set
% Adapted from test_recon_dataset supplied with ISMRMRD
%
% This is a simple example of how to reconstruct images from data
% acquired on a fully sampled cartesian grid
%
% Capabilities:
%   2D/3D
%   multiple slices/slabs
%   multiple contrasts, repetitions
%   
% Limitations:
%   only works with a single encoded space
%   fully sampled k-space (no partial fourier or undersampling)
%   multiple repetitions
%   doesn't handle averages, phases, segments and sets
%   ignores noise scans (no pre-whitening)
% 

% You can first create a data set using the example program like this:
%   ismrmrd_generate_cartesian_shepp_logan -r 5 -C -o shepp-logan.h5
% This will produce the file shepp-logan.h5 containing an ISMRMRD
% dataset sampled evenly on the k-space grid -128:0.5:127.5 x -128:127
% (i.e. oversampled by a factor of 2 in the readout direction)
% with 8 coils, 5 repetitions and a noise level of 0.5
% with a noise calibration scan at the beginning
%
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Loading an existing file %
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
filename = pref_uigetfile('explore_recon_dataset','filename',varargin{:}) ;

if exist(filename, 'file')
    dset = ismrmrd.Dataset(filename, 'dataset');
else
    error(['File ' filename ' does not exist.  Please generate it.'])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read some fields from the XML header %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We need to check if optional fields exists before trying to read them

hdr = ismrmrd.xml.deserialize(dset.readxml);

% Encoding and reconstruction information
% Matrix size
enc_Nx = hdr.encoding.encodedSpace.matrixSize.x;
enc_Ny = hdr.encoding.encodedSpace.matrixSize.y;
enc_Nz = hdr.encoding.encodedSpace.matrixSize.z;
rec_Nx = hdr.encoding.reconSpace.matrixSize.x;
rec_Ny = hdr.encoding.reconSpace.matrixSize.y;
rec_Nz = hdr.encoding.reconSpace.matrixSize.z;

% Field of View
enc_FOVx = hdr.encoding.encodedSpace.fieldOfView_mm.x;
enc_FOVy = hdr.encoding.encodedSpace.fieldOfView_mm.y;
enc_FOVz = hdr.encoding.encodedSpace.fieldOfView_mm.z;
rec_FOVx = hdr.encoding.reconSpace.fieldOfView_mm.x;
rec_FOVy = hdr.encoding.reconSpace.fieldOfView_mm.y;
rec_FOVz = hdr.encoding.reconSpace.fieldOfView_mm.z;

% Number of slices, coils, repetitions, contrasts etc.
% We have to wrap the following in a try/catch because a valid xml header may
% not have an entry for some of the parameters

try
  nSlices = hdr.encoding.encodingLimits.slice.maximum + 1;
catch
    nSlices = 1;
end

try 
    nCoils = hdr.acquisitionSystemInformation.receiverChannels;
catch
    nCoils = 1;
end

try
    nReps = hdr.encoding.encodingLimits.repetition.maximum + 1;
catch
    nReps = 1;
end

try
    nContrasts = hdr.encoding.encodingLimits.contrast.maximum + 1 + 1;
catch
    nContrasts = 1;
end

% TODO add the other possibilites

disp(filename)
disp( sprintf('nSlices: %d, nCoils: %d, nReps: %d, nContrasts: %d',...
    nSlices, nCoils, nReps, nContrasts)) ;

% Read all the data
% Reading can be done one acquisition (or chunk) at a time, 
% but this is much faster for data sets that fit into RAM.
D = dset.readAcquisition();

% Note: can select a single acquisition or header from the block, e.g.
% acq = D.select(5);
% hdr = D.head.select(5);
% or you can work with them all together

% Ignore noise scans
% TODO add a pre-whitening example
% Find the first non-noise scan
% This is how to check if a flag is set in the acquisition header
isNoise = D.head.flagIsSet('ACQ_IS_NOISE_MEASUREMENT');
disp([' There are ',num2str(sum(isNoise)),' noise measure(s).'])
firstScan = find(isNoise==0,1,'first');
if firstScan > 1
    noise = D.select(1:firstScan-1);
else
    noise = [];
end

%DA: 
if isempty(firstScan)
    disp(['No non-noise scan detected, assume firstScan is 1'])
    firstScan = 1 ;
end

meas  = D.select(firstScan:D.getNumber);

% from more recent recon_ismrmrd_hdf5
uquser = unique(meas.head.idx.user','rows') ;
if size(uquser,1) > 1
    warning(['More than one user entry,using first'])
end

clear D;

% Reconstruct images
% Since the entire file is in memory we can use random access
% Loop over repetitions, contrasts, slices
reconImages = {};
nimages = 0;
for rep = 1:nReps
    for contrast = 1:nContrasts
        for slice = 1:nSlices
            % Initialize the K-space storage array
            K = zeros(enc_Nx, enc_Ny, enc_Nz, nCoils);
            % Select the appropriate measurements from the data
            acqs = find(  (meas.head.idx.contrast==(contrast-1)) ...
                        & (meas.head.idx.repetition==(rep-1)) ...
                        & (meas.head.idx.slice==(slice-1)) ...
                        & (ismember(meas.head.idx.user', uquser(1,:),'rows')') ...  % DA was uquser(1,:) - why??
                        & (meas.head.idx.average==0) );  % DA
            for p = 1:length(acqs)
                ky = meas.head.idx.kspace_encode_step_1(acqs(p)) + 1;
                kz = meas.head.idx.kspace_encode_step_2(acqs(p)) + 1;
                K(:,ky,kz,:) = meas.data{acqs(p)};
            end
            % Reconstruct in x
            K = fftshift(ifft(fftshift(K,1),[],1),1);
            % Chop if needed
            if (enc_Nx == rec_Nx)
                im = K;
            else
                ind1 = floor((enc_Nx - rec_Nx)/2)+1;
                ind2 = floor((enc_Nx - rec_Nx)/2)+rec_Nx;
                im = K(ind1:ind2,:,:,:);
            end
            % Reconstruct in y then z
            im = fftshift(ifft(fftshift(im,2),[],2),2);
            if size(im,3)>1
                im = fftshift(ifft(fftshift(im,3),[],3),3);
            end
            
            % Combine SOS across coils
            im = sqrt(sum(abs(im).^2,4));
            
            % Append
            nimages = nimages + 1;
            reconImages{nimages} = im;
        end
    end
end

% Display an image
figure
colormap gray
imdisp = ceil(nimages/2) ;
imagesc(reconImages{imdisp}); axis image; axis off; colorbar;