Skip to content
Permalink
master
Switch branches/tags

aedes/aedes_readfid.m

Go to file
 
 
Cannot retrieve contributors at this time
1865 lines (1661 sloc) 59.4 KB
Raw Blame
function [DATA,msg_out] = aedes_readfid(filename,varargin)
% AEDES_READFID - Read VNMR (Varian) FID-files
%
%
% Synopsis:
% [DATA,msg]=aedes_readfid(filename,'PropertyName1',value1,'PropertyName2',value2,...)
% [DATA,msg]=aedes_readfid(filename,'header')
% [DATA,msg]=aedes_readfid(filename,output_filename)
%
% Description:
% The function reads VNMR (Varian) FID-file and return a structure
% DATA with fields DataFormat, HDR, FTDATA, KSPACE, PROCPAR, and
% PHASETABLE. The fields of the DATA structure are constructed as
% follows:
%
% DATA
% |-> DataFormat (string identifier for data format 'vnmr')
% |-> HDR
% |-> FileHeader (data file header)
% |-> BlockHeaders (data block headers, not stored by default)
% |-> fname (file name)
% |-> fpath (file path)
% |-> Param (parameter values used by AEDES_READFID to read the data)
% |-> FTDATA (real fourier transformed image data)
% |-> KSPACE (complex k-space, empty by default)
% |-> PROCPAR (parameters from procpar file)
% |-> PHASETABLE (phasetable)
%
% The DATA structure is returned as empty (without HDR, FTDATA,
% KSPACE, and PROCPAR fields) if an error has occurred during
% reading. The error message is returned in the second output
% argument msg. If AEDES_READFID is called with only one output argument
% (i.e. without MSG), the error message is echoed to the workspace.
%
% The first input argument is either a string containing full path to
% the FID-file or the header structure HDR. Only the data file header
% can be read by giving a string 'header' as the second input
% argument.
%
% By default the k-space is not returned, i.e. the field KSPACE is
% empty. The returned data can be adjusted by using the 'return'
% property and values 1, 2, or 3 (see below for more information).
%
% The supported property-value pairs in AEDES_READFID (property strings
% are not case sensitive):
%
% Property: Value: Description:
% ********* ****** ************
% 'Return' : [ {1} | 2 | 3 | 4 ] % 1=return only ftdata (default)
% % 2=return only k-space
% % 3=return both ftdata & kspace
% % 4=return raw kspace
%
% 'FastRead' : [{'off'} | 'on' ] % Enables an experimental
% % method for very fast reading
% % of fid-files. This is not as
% % robust as the older
% % method and can consume a lot
% % of memory.
%
% 'FileChunkSize' : [ megabytes ] % Chunk size in megabytes for
% % reading fid-files
% % (default=500 MB).
% % Lowering the chunk size
% % helps to conserve memory
% % when reading large files,
% % but is slower. This
% % property has no effect if
% % FastRead='off'.
%
%
% 'OutputFile' : filename % Writes the slices into
% % NIfTI-format files
% % using the given filename.
%
% 'DCcorrection': [ {'off'} | 'on' ] % 'on'=use DC correction
% % 'off'=don't use DC correction
% % (default)
%
% 'Procpar' : (procpar-structure) % Input procpar
% % structure. If this
% % property is not set the
% % procpar structure
% % will be read using
% % AEDES_READPROCPAR
%
% 'ZeroPadding' : ['off'|'on'|{'auto'}] % 'off' = force
% % zeropadding off
% % 'on' = force
% % zeropadding on (force
% % images to be square)
% % 'auto' = autoselect
% % using relative FOV
% % dimensions (default)
%
% 'PhaseTable' : (custom phase table) % User-specified
% % line order for
% % k-space. The phase table
% % is obtained from the file
% % specified in
% % procpar.petable if
% % necessary.
%
% 'sorting' : [ 'off' | {'on'} ] % 'off' =disable k-space
% % sorting, i.e. ignore the
% % phase table and/or arrays.
% % 'on' =sort k-space using
% % phase table and/or array
% % information if necessary
% % (default)
%
% 'wbar' : [ {'on'} | 'off' ] % 'on' = show waitbar
% % 'off' = don't show waitbar
%
% 'FlipKspace' : [ {'off'} | 'LR' |
% 'UD' | 'LRUD' ] % 'off' = don't flip (default)
% % 'LR' = left/right
% % 'UD' = up/down
% % 'LRUD' = left/right and
% % up/down
%
% 'FlipInd' : [ {'all'} | 'alt' |
% (custom vector) ] % 'all' = flip all slices
% % (default)
% % 'alt' = flip every
% % second slice
% % (custom vector) = A
% % custom vector containing
% % indices to the flipped slices
%
% 'Precision' : ['single'|{'double'}] % Precision of the
% % outputted data.
% % 'single' = 32-bit float
% % 'double' = 64-bit float
%
% 'OrientImages': [ {'on'} | 'off' ] % Orient FTDATA after
% % reading the data using
% % PROCPAR.orient property.
%
% 'RemoveEPIphaseIm' : [{'on'}|'off'] % Remove the phase image
% % from EPI data. This
% % option only affects EPI
% % images.
%
% 'EPIBlockSize' : [integer value] % Block size (number of
% % volumes) used for
% % processing multireceiver
% % EPI data. Default=100
%
% 'EPIPhasedArrayData' : ['on'|{'off'}] % Return data from
% % individual receivers from
% % phased array EPI files.
%
% Examples:
% [DATA,msg]=aedes_readfid(filename) % Read image data from 'filename'
% [DATA,msg]=aedes_readfid(filename,'header') % Read only data file header
%
% [DATA,msg]=aedes_readfid(filename,'return',1) % Return only image data (default)
% [DATA,msg]=aedes_readfid(filename,'return',2) % Return only k-space
% [DATA,msg]=aedes_readfid(filename,'return',3) % Return both image data and k-space
%
% % Read first data file header and then image data and k-space
% [DATA,msg]=aedes_readfid(filename,'header')
% [DATA,msg]=aedes_readfid(DATA.HDR,'return',3)
%
% % Read VNMR-data with default options and save slices in NIfTI
% % format
% DATA=aedes_readfid('example.fid','saved_slices.nii');
%
% % If you want to use non-default options and also write
% % NIfTI-files, use the property/value formalism, for example
% DATA=aedes_readfid('example.fid','ZeroPadding','off',...
% 'OutputFile','saved_slices.nii');
%
% See also:
% AEDES_READFIDPREFS, AEDES_READPROCPAR, AEDES_READPHASETABLE,
% AEDES_DATA_READ, AEDES_WRITE_NIFTI, AEDES
% This function is a part of Aedes - A graphical tool for analyzing
% medical images
%
% Copyright (C) 2006 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uku.fi>
%
% Department of Physics, Department of Neurobiology
% University of Kuopio, FINLAND
%
% This program may be used under the terms of the GNU General Public
% License version 2.0 as published by the Free Software Foundation
% and appearing in the file LICENSE.TXT included in the packaging of
% this program.
%
% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
if nargout==2
Dat.ShowErrors = false;
msg_out='';
else
Dat.ShowErrors = true;
end
%% ---------------------
% Defaults
Dat.ReturnKSpace = false;
Dat.ReturnFTData = true;
Dat.DCcorrection = false;
Dat.ZeroPadding = 2;
Dat.Sorting = true;
Dat.UsePhaseTable = true;
Dat.FastDataRead = true;
Dat.precision = 'single';
Dat.FileChunkSize = 500; % Chunk size (in MB) for FastRead
%% Other defaults
Dat.ShowWaitbar = true;
procpar=[];
Dat.phasetable=[];
Dat.FlipKspace = 0;
Dat.FlipInd = 'all';
Dat.OutputFile = false;
Dat.ReturnRawKspace = false;
Dat.ReorientEPI = false;
Dat.RemoveEPIphaseIm = true;
Dat.EPIBlockSize = 100;
Dat.EPIPhasedArrayData = false;
Dat.OrientImages = true;
%% -------------------------------------------------------------------
%% Set data format label
DATA.DataFormat = 'vnmr';
%% Parse input arguments
if nargin==0 || isempty(filename)
%% Get default directory
try
defaultdir = getpref('Aedes','GetDataFileDir');
catch
defaultdir = [pwd,filesep];
end
%% If no input arguments are given, ask for a file
[f_name, f_path, f_index] = uigetfile( ...
{'fid;FID','Varian FID-files (fid)'; ...
'*.*', 'All Files (*.*)'}, ...
'Select VNMR (Varian) FID file',defaultdir);
if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed
DATA=[];
msg_out='Canceled';
return
end
ReadHdr = true;
ReadData = true;
filename=[f_path,f_name];
%% Set default directory to preferences
setpref('Aedes','GetDataFileDir',f_path)
end
if nargin==1
if isstruct(filename)
hdr = filename;
filename = [hdr.fpath,hdr.fname];
ReadHdr = false;
ReadData = true;
% $$$ ReturnKSpace = false;
% $$$ ReturnFTData = true;
elseif ischar(filename)
ReadHdr = true;
ReadData = true;
% $$$ ReturnKSpace = false;
% $$$ ReturnFTData = true;
end
elseif nargin==2
if strcmpi(varargin{1},'header')
ReadHdr = true;
ReadData = false;
% $$$ ReturnKSpace = false;
elseif ischar(varargin{1})
ReadHdr = true;
ReadData = true;
Dat.OutputFile = varargin{1};
else
if ~Dat.ShowErrors
DATA=[];
msg_out=sprintf('Invalid second input argument "%s".',varargin{1});
return
else
error('Invalid second input argument "%s".',varargin{1})
end
end
else
if isstruct(filename)
hdr = filename;
filename = [hdr.fpath,hdr.fname];
ReadHdr = false;
ReadData = true;
elseif isempty(filename)
[f_name, f_path, f_index] = uigetfile( ...
{'fid;FID','Varian FID-files (fid)'; ...
'*.*', 'All Files (*.*)'}, ...
'Select VNMR (Varian) FID file');
if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed
DATA=[];
msg_out='Canceled';
return
end
ReadHdr = true;
ReadData = true;
filename=[f_path,f_name];
else
ReadHdr = true;
ReadData = true;
end
for ii=1:2:length(varargin)
switch lower(varargin{ii})
case 'return'
if length(varargin)<2
DATA=[];
if ~Dat.ShowErrors
msg_out='"Return" value not specified!';
else
error('"Return" value not specified!')
end
return
else
if varargin{ii+1}==1
Dat.ReturnKSpace = false;
Dat.ReturnFTData = true;
elseif varargin{ii+1}==2
Dat.ReturnKSpace = true;
Dat.ReturnFTData = false;
elseif varargin{ii+1}==3
Dat.ReturnKSpace = true;
Dat.ReturnFTData = true;
elseif varargin{ii+1}==4
Dat.ReturnRawKspace = true;
Dat.ReturnKSpace = true;
Dat.ReturnFTData = false;
end
end
case {'dc','dccorrection','dccorr'}
if strcmpi(varargin{ii+1},'on')
Dat.DCcorrection = true;
else
Dat.DCcorrection = false;
end
case 'procpar'
procpar=varargin{ii+1};
case 'zeropadding'
if ischar(varargin{ii+1})
if strcmpi(varargin{ii+1},'on')
Dat.ZeroPadding = 1; % on
elseif strcmpi(varargin{ii+1},'off')
Dat.ZeroPadding = 0; % off
else
Dat.ZeroPadding = 2; % auto
end
else
% Undocumented
Dat.ZeroPadding = 3; % Custom
Dat.CustomZeroPadding = varargin{ii+1};
end
case 'phasetable'
Dat.phasetable = varargin{ii+1};
case 'sorting'
if strcmpi(varargin{ii+1},'on')
Dat.UsePhaseTable = true;
Dat.Sorting = true;
else
Dat.UsePhaseTable = false;
Dat.Sorting = false;
end
case 'fastread'
if strcmpi(varargin{ii+1},'on')
Dat.FastDataRead = true;
else
Dat.FastDataRead = false;
end
case 'filechunksize'
Dat.FileChunkSize = round(varargin{ii+1});
case 'wbar'
if strcmpi(varargin{ii+1},'on')
Dat.ShowWaitbar = 1;
else
Dat.ShowWaitbar = 0;
end
case 'precision'
if strcmpi(varargin{ii+1},'single')
Dat.precision = 'single';
end
case 'flipkspace'
flipkspace = varargin{ii+1};
if strcmpi(flipkspace,'off')
Dat.FlipKspace=0;
elseif strcmpi(flipkspace,'LR')
Dat.FlipKspace=1;
elseif strcmpi(flipkspace,'UD')
Dat.FlipKspace=2;
elseif strcmpi(flipkspace,'LRUD')
Dat.FlipKspace=3;
else
DATA=[];
if ~Dat.ShowErrors
msg_out = 'Bad value for property FlipKspace!';
else
error('Bad value for property FlipKspace!')
end
return
end
case 'flipind'
Dat.FlipInd = varargin{ii+1};
case 'outputfile'
Dat.OutputFile = varargin{ii+1};
case 'reorientepi'
if strcmpi(varargin{ii+1},'on')
Dat.ReorientEPI = true;
end
case 'removeepiphaseim'
if strcmpi(varargin{ii+1},'on')
Dat.RemoveEPIphaseIm = true;
end
case 'epiblocksize'
Dat.EPIBlockSize = round(varargin{ii+1});
case 'epiphasedarraydata'
if strcmpi(varargin{ii+1},'on')
Dat.EPIPhasedArrayData = true;
else
Dat.EPIPhasedArrayData = false;
end
case 'orientimages'
if strcmpi(varargin{ii+1},'off')
Dat.OrientImages = false;
end
otherwise
DATA=[];
if ~Dat.ShowErrors
msg_out = ['Unknown property "' varargin{ii} '"'];
else
error(['Unknown property "' varargin{ii} '"'])
end
return
end
end
end
% Parse filename
[fpath,fname,fext]=fileparts(filename);
if ~strcmpi(fname,'fid')
if isempty(fname)
fpath = [fpath,filesep];
else
if isempty(fpath)
fpath = [pwd,filesep,fname,fext,filesep];
else
fpath = [fpath,filesep,fname,fext,filesep];
end
end
fname = 'fid';
else
fpath = [fpath,filesep];
end
%% Read procpar if not given as input argument
if isempty(procpar) || nargin~=2
[procpar,proc_msg]=aedes_readprocpar([fpath,'procpar']);
if isempty(procpar)
DATA=[];
if ~Dat.ShowErrors
msg_out=proc_msg;
else
error(proc_msg)
end
return
end
end
%% Read phasetable if necessary
if isfield(procpar,'petable') && isempty(Dat.phasetable) && ...
~isempty(procpar.petable{1}) && ~strcmpi(procpar.petable{1},'n') && ...
Dat.Sorting
% Look in preferences for tablib-directory
try
tabpath = getpref('Aedes','TabLibDir');
catch
% If the table path was not found in preferences, try to use aedes
% directory
tmp_path = which('aedes');
if ~isempty(tmp_path)
aedes_path=fileparts(tmp_path);
tabpath = [aedes_path,filesep,'tablib',filesep];
else
% If all else fails, look in the current directory
tabpath = [pwd,filesep,'tablib',filesep];
end
end
[phasetable,msg] = aedes_readphasetable([tabpath,procpar.petable{1}]);
% If petable cannot be found, check if it is in procpar...
if isempty(phasetable) && isfield(procpar,'pe_table')
phasetable = {'t1',procpar.pe_table(:)};
elseif isempty(phasetable) && isfield(procpar,'pelist') && ...
~isempty(procpar.pelist) && isnumeric(procpar.pelist)
phasetable = {'t1',reshape(procpar.pelist,procpar.etl,[]).'};
end
% If the sequence is a fast spin echo, try to construct the phasetable
if isempty(phasetable) && isfield(procpar,'petable') && ...
strncmpi(procpar.petable{1},'fse',3)
err_msg = 'Could not find phasetable!';
if ~Dat.ShowErrors
msg_out=err_msg;
else
error(err_msg)
end
return
%phasetable = {'t1',l_CreateFsemsPhaseTable(procpar)};
end
%% Abort and throw error if phasetable cannot be read and the FID-file
% has not been sorted
if isempty(phasetable) && ~isfield(procpar,'flash_converted')
DATA=[];
if ~Dat.ShowErrors
msg_out={['Could not find the required phase table "' ...
procpar.petable{1} '" in the following folders'],...
['"' tabpath '".']};
else
error('Could not find the required phase table "%s" in %s',...
procpar.petable{1},tabpath)
end
return
elseif ( isempty(phasetable) && isfield(procpar,'flash_converted') ) || ...
~Dat.UsePhaseTable
%% If the FID-file has been sorted, the reading can continue but
% throw a warning.
fprintf(1,'Warning: aedes_readfid: Could not find phasetable "%s" in "%s"!\n',procpar.petable{1},tabpath)
Dat.phasetable=[];
else
Dat.phasetable = phasetable{1,2};
end
end
% Convert phasetable indices to MATLAB indices
if ~isempty(Dat.phasetable)
Dat.phasetable=Dat.phasetable+(-min(min(Dat.phasetable))+1);
else
Dat.UsePhaseTable = false;
end
%% Open FID-file
[file_fid,msg]=fopen([fpath,fname],'r','ieee-be');
if file_fid < 0
DATA=[];
if ~Dat.ShowErrors
msg_out=['Could not open file "' filename '" for reading.'];
else
error('Could not open file "%s" for reading.',filename)
end
return
end
% Read only header
if ReadHdr && ~ReadData
[hdr,msg_out]=l_ReadDataFileHeader(file_fid);
if ~isempty(msg_out)
DATA=[];
fclose(file_fid);
if Dat.ShowErrors
error(msg_out)
end
return
else
DATA.HDR.FileHeader = hdr.FileHeader;
DATA.FTDATA=[];
DATA.KSPACE=[];
DATA.PROCPAR=[];
DATA.PHASETABLE=[];
end
elseif ~ReadHdr && ReadData % Header structure given as input argument
% Read only data.
[hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar);
if ~isempty(msg_out)
DATA=[];
fclose(file_fid);
if Dat.ShowErrors
error(msg_out)
end
return
else
DATA.HDR.FileHeader=hdr.FileHeader;
DATA.HDR.BlockHeaders = hdr.BlockHeaders;
DATA.FTDATA=data;
DATA.KSPACE=kspace;
DATA.PROCPAR=procpar;
DATA.PHASETABLE=Dat.phasetable;
end
elseif ReadHdr && ReadData % Read headers and data
[hdr,msg_out]=l_ReadDataFileHeader(file_fid);
[hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar);
if ~isempty(msg_out)
DATA=[];
fclose(file_fid);
if Dat.ShowErrors
error(msg_out)
end
return
else
DATA.HDR.FileHeader=hdr.FileHeader;
DATA.HDR.BlockHeaders = hdr.BlockHeaders;
DATA.FTDATA=data;
DATA.KSPACE=kspace;
DATA.PROCPAR=procpar;
DATA.PHASETABLE=Dat.phasetable;
end
end
% Set file name and path to the HDR structure
DATA.HDR.fname=fname;
DATA.HDR.fpath=fpath;
% Set parameter values
DATA.HDR.Param.ReturnKSpace = Dat.ReturnKSpace;
DATA.HDR.Param.ReturnFTData = Dat.ReturnFTData;
if Dat.ZeroPadding==0
DATA.HDR.Param.ZeroPadding = 'off';
elseif Dat.ZeroPadding==1
DATA.HDR.Param.ZeroPadding = 'on';
else
DATA.HDR.Param.ZeroPadding = 'auto';
end
if ~Dat.DCcorrection
DATA.HDR.Param.DCcorrection = 'off';
else
DATA.HDR.Param.DCcorrection = 'on';
end
if Dat.Sorting==0
DATA.HDR.Param.Sorting = 'off';
else
DATA.HDR.Param.Sorting = 'on';
end
if Dat.FlipKspace==0
DATA.HDR.Param.FlipKspace = 'off';
elseif Dat.FlipKspace==1
DATA.HDR.Param.FlipKspace = 'LR';
elseif Dat.FlipKspace==2
DATA.HDR.Param.FlipKspace = 'UD';
elseif Dat.FlipKspace==3
DATA.HDR.Param.FlipKspace = 'LRUD';
end
DATA.HDR.Param.FlipInd = Dat.FlipInd;
% Close file
fclose(file_fid);
%% Write slices to NIfTI files
if ischar(Dat.OutputFile)
if isempty(Dat.OutputFile)
[fp,fn,fe]=fileparts(DATA.HDR.fpath(1:end-1));
fprefix=fn;
niipath = [pwd,filesep];
else
[fp,fn,fe]=fileparts(Dat.OutputFile);
if strcmpi(fe,'.nii')
fprefix=fn;
else
fprefix=[fn,fe];
end
if isempty(fp)
niipath = [pwd,filesep];
else
niipath = [fp,filesep];
end
end
% Create file names
filenames={};
for ii=1:size(DATA.FTDATA,3)
filenames{ii}=sprintf('%s%03d%s',[fprefix,'_'],ii,'.nii');
end
nFiles = length(filenames);
h=aedes_wbar(0,sprintf('Saving slice 1/%d in NIfTI format...',nFiles));
for ii=1:nFiles
aedes_wbar(ii/nFiles,h,sprintf('Saving slice %d/%d in NIfTI format...',ii, ...
nFiles))
[done,msg]=aedes_write_nifti(DATA.FTDATA(:,:,ii),...
[niipath,filenames{ii}],'DataType','single',...
'FileType',2);
end
delete(h)
if ~done
warning('Error occurred while writing NIfTI-files. Could not write file(s)!')
end
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read Data File (Main) Header
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hdr,msg_out]=l_ReadDataFileHeader(file_fid)
msg_out='';
%% Read Data File Header
try
FH.nblocks = fread(file_fid,1,'long'); % Number of blocks in file
FH.ntraces = fread(file_fid,1,'long'); % Number of traces per block
FH.np = fread(file_fid,1,'long'); % Number of elements per trace
FH.ebytes = fread(file_fid,1,'long'); % Number of bytes per element
FH.tbytes = fread(file_fid,1,'long'); % Number of bytes per trace
FH.bbytes = fread(file_fid,1,'long'); % Number of bytes per block
FH.vers_id = fread(file_fid,1,'short'); % Software version, file_id status bits
FH.status = fread(file_fid,1,'short'); % Status of whole file
FH.nbheaders = fread(file_fid,1,'long'); % Number of block headers per block
hdr.FileHeader = FH;
%% Parse status bits
hdr.FileHeader.status=[];
tmp_str = {'S_DATA',... % 0 = no data, 1 = data
'S_SPEC',... % 0 = FID, 1 = spectrum
'S_32',... % 0 = 16 bit, 1 = 32 bit integer
'S_FLOAT',... % 0 = integer, 1 = floating point
'S_COMPLEX',... % 0 = real, 1 = complex
'S_HYPERCOMPLEX',... % 1 = hypercomplex
'',... % Unused bit
'S_ACQPAR',... % 0 = not Acqpar, 1 = Acqpar
'S_SECND',... % 0 = first FT, 1 = second FT
'S_TRANSF',... % 0 = regular, 1 = transposed
'',... % Unused bit
'S_NP',... % 1 = np dimension is active
'S_NF',... % 1 = nf dimension is active
'S_NI',... % 1 = ni dimension is active
'S_NI2',... % 1 = ni2 dimension is active
''... % Unused bit
};
status_bits = fliplr(double(dec2bin(FH.status,16))==49);
for ii=1:length(tmp_str)
if ~isempty(tmp_str{ii})
hdr.FileHeader.status.(tmp_str{ii}) = status_bits(ii);
end
end
catch
hdr=[];
msg_out=['Error occurred while reading file header from "' filename '".'];
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read Block Headers and Data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar)
hdr.BlockHeaders = [];
%hdr.HypercmplxBHeaders = [];
data=[];
kspace=[];
count = [];
msg_out='';
if isfield(procpar,'seqcon')
procpar.seqcon=procpar.seqcon{1};
else
procpar.seqcon='nnnnn';
end
%% Determine acquisition type from seqcon parameter
if all(procpar.seqcon=='n') % 1D spectroscopy
Dat.AcqType=1;
elseif all(procpar.seqcon(4:5)=='n') % 2D imaging
Dat.AcqType=2;
elseif procpar.seqcon(5)=='n' % 3D imaging
Dat.AcqType=3;
else % 4D imaging
Dat.AcqType=4;
end
%% Double-check that AcqType is not 1D spectral
if isfield(procpar,'nv') && procpar.nv==0
Dat.AcqType=1;
end
%% Use some heuristical approach to "triple-check" that the data is not
% 1D spectral
if ~isempty(strfind(procpar.seqfil{1},'STEAM')) || ...
~isempty(strfind(procpar.seqfil{1},'steam')) || ...
~isempty(strfind(procpar.seqfil{1},'LASER')) || ...
~isempty(strfind(procpar.seqfil{1},'laser')) || ...
~isempty(strfind(procpar.seqfil{1},'PRESS')) || ...
~isempty(strfind(procpar.seqfil{1},'press')) || ...
~isempty(strfind(procpar.seqfil{1},'1PULSE')) || ...
~isempty(strfind(procpar.seqfil{1},'1pulse')) || ...
~isempty(strfind(procpar.seqfil{1},'2PULSE')) || ...
~isempty(strfind(procpar.seqfil{1},'2pulse'))
Dat.AcqType=1;
end
UsePhaseTable=Dat.UsePhaseTable;
%% If the data has been converted with flashc, the seqcon parameter needs
% to be changed
if isfield(procpar,'flash_converted')
if isfield(procpar,'ni') && procpar.ni>1
procpar.seqcon(3)='s';
elseif ((procpar.seqcon(2)=='c') && (procpar.seqcon(3)=='c'))
procpar.seqcon(2)='s';
end
UsePhaseTable=false; % Do not try to order data if flashc has been run
end
%% Set number of transients
if ~isfield(procpar,'nt')
nt=1;
else
nt=procpar.nt;
end
%% Determine if the arrayed acquisition was used
if isfield(procpar,'array') && ~isempty(procpar.array{1})
if length(procpar.array)==1 && ~iscell(procpar.array{1}) && ...
strcmp(procpar.array{1},'pad') && all(procpar.pad==0)
% Skip the array parsing if the array is a dummy "pad"-array...
Dat.isAcqArrayed = false;
Dat.ArrayLength = 1;
else
Dat.isAcqArrayed = true;
Dat.ArrayLength = [];
% Determine array length
for ii=1:length(procpar.array)
if iscell(procpar.array{ii})
Dat.ArrayLength(ii) = length(procpar.(procpar.array{ii}{1}));
else
Dat.ArrayLength(ii) = length(procpar.(procpar.array{ii}));
end
end
Dat.ArrayLength = prod(Dat.ArrayLength);
end
else
Dat.isAcqArrayed = false;
Dat.ArrayLength = 1;
end
%% Determine if the data is EPI data
if ( isfield(procpar,'readres') && isfield(procpar,'phaseres') ) || ...
( isfield(procpar,'apptype') && strcmp(procpar.apptype{1},'imEPI') )
Dat.isEPIdata = true;
else
Dat.isEPIdata = false;
end
BlockHeadStatusLabels = {'S_DATA',... % 0 = no data, 1 = data
'S_SPEC',... % 0 = FID, 1 = spectrum
'S_32',... % 0 = 16 bit, 1 = 32 bit integer
'S_FLOAT',... % 0 = integer, 1 = floating point
'S_COMPLEX',... % 0 = real, 1 = complex
'S_HYPERCOMPLEX',... % 1 = hypercomplex
'',... % Unused bit
'MORE_BLOCKS',... % 0 = absent, 1 = present
'NP_CMPLX',... % 0 = real, 1 = complex
'NF_CMPLX',... % 0 = real, 1 = complex
'NI_CMPLX',... % 0 = real, 1 = complex
'NI2_CMPLX',... % 0 = real, 1 = complex
'',... % Unused bit
'',... % Unused bit
'',... % Unuesd bit
''... % Unused bit
};
BlockHeadModeLabels = {'NP_PHMODE',... % 1 = ph mode
'NP_AVMODE',... % 1 = av mode
'NP_PWRMODE',... % 1 = pwr mode
'',... % Unused bit
'NF_PHMODE',... % 1 = ph mode
'NF_AVMODE',... % 1 = av mode
'NF_PWRMODE',... % 1 = pwr mode
'',... % Unused bit
'NI_PHMODE',... % 1 = ph mode
'NI_AVMODE',... % 1 = av mode
'NI_PWRMODE',... % 1 = pwr mode
'',... % Unused bit
'NI2_PHMODE',... % 1 = ph mode
'NI2_AVMODE',... % 1 = av mode
'NI2_PWRMODE',... % 1 = pwr mode
''... % Unused bit
};
% The nbheaders -field is not correct in some cases. Thus, this field
% cannot be trusted and the real nbheaders has to be calculated from the
% byte counts.
tmp_bytes=hdr.FileHeader.bbytes-hdr.FileHeader.tbytes*hdr.FileHeader.ntraces;
header_bytes=28;
if rem(tmp_bytes,header_bytes)~=0
msg_out = 'Block header byte count does not match with file header';
return
else
nbheaders = tmp_bytes/28;
end
%nbheaders = hdr.FileHeader.nbheaders;
%% Allocate space for k-space
% kspace=zeros(hdr.FileHeader.np/2,...
% hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);
if ~Dat.FastDataRead
if any(Dat.AcqType==[1 2])
switch procpar.seqcon(2:3)
case {'cc','sc'}
kspace = complex(zeros(hdr.FileHeader.np/2,...
hdr.FileHeader.ntraces,...
hdr.FileHeader.nblocks,Dat.precision));
otherwise
kspace = complex(zeros(hdr.FileHeader.np/2,...
hdr.FileHeader.nblocks,...
hdr.FileHeader.ntraces,Dat.precision));
end
else
kspace = complex(zeros(hdr.FileHeader.np/2,...
hdr.FileHeader.ntraces,...
hdr.FileHeader.nblocks,Dat.precision));
end
else
%kspace = [];
kspace = complex(zeros(hdr.FileHeader.np/2*hdr.FileHeader.ntraces,...
hdr.FileHeader.nblocks,Dat.precision));
end
%% - The older robust (but also slower) way for reading the data.
%% When the blocksize is relatively small, this is also quite fast.
if ~Dat.FastDataRead
% Initialize waitbar
if Dat.ShowWaitbar
wb_h = aedes_wbar(0/hdr.FileHeader.nblocks,...
{['Reading ',num2str(Dat.AcqType),'D VNMR data (seqcon: "' procpar.seqcon '")'],...
['(Processing data block ' ...
num2str(0) '/' num2str(hdr.FileHeader.nblocks) ')']});
end
%% Read data and block headers
for ii=1:hdr.FileHeader.nblocks
%% Update waitbar
if Dat.ShowWaitbar
aedes_wbar(ii/hdr.FileHeader.nblocks,...
wb_h,...
{['Reading ',num2str(Dat.AcqType),'D VNMR data (seqcon: "' procpar.seqcon '")'],...
['(Processing data block ' ...
num2str(ii) '/' num2str(hdr.FileHeader.nblocks) ')']})
end
%% Read block header and hypercomplex header
for kk=1:nbheaders
%% Read block header
if kk==1
hdr.BlockHeaders.scale = fread(file_fid,1,'short'); % Scaling factor
tmp_status = fread(file_fid,1,'short'); % Status of data in block
hdr.BlockHeaders.status = [];
hdr.BlockHeaders.index = fread(file_fid,1,'short'); % Block index
tmp_mode = fread(file_fid,1,'short'); % Mode of data in block
hdr.BlockHeaders.mode = [];
hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
hdr.BlockHeaders.lpval = fread(file_fid,1,'float'); % f2 (2D-f1) left phase in phasefile
hdr.BlockHeaders.rpval = fread(file_fid,1,'float'); % f2 (2D-f1) right phase in phasefile
hdr.BlockHeaders.lvl = fread(file_fid,1,'float'); % level drift correction
hdr.BlockHeaders.tlt = fread(file_fid,1,'float'); % tilt drift correction
%% Parse status and mode bits
status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
for tt=1:length(BlockHeadStatusLabels)
if ~isempty(BlockHeadStatusLabels{tt})
hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
end
if ~isempty(BlockHeadModeLabels{tt})
hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
end
end
else %% Read hypercomplex header
fread(file_fid,1,'short'); % Spare
hdr.BlockHeaders.HCHstatus = fread(file_fid,1,'short');
fread(file_fid,1,'short'); % Spare
fread(file_fid,1,'short'); % Spare
fread(file_fid,1,'long'); % Spare
hdr.BlockHeaders.HCHlpval1 = fread(file_fid,1,'float');
hdr.BlockHeaders.HCHrpval1 = fread(file_fid,1,'float');
fread(file_fid,1,'float'); % Spare
fread(file_fid,1,'float'); % Spare
end
end
%% Check block index to be sure about the data type
if hdr.BlockHeaders.index~=ii
fprintf(1,'Warning: Index mismatch in "%s" block %d\n',fopen(file_fid),ii);
% Use information from the file header instead of the BlockHeader if
% there is a mismatch in blockheader index...
useFileHeader = true;
else
useFileHeader = false;
end
%% Determine data precision
if useFileHeader
if hdr.FileHeader.status.S_FLOAT==1
prec_str = ['single=>',Dat.precision];
elseif hdr.FileHeader.status.S_32==1 ...
&& hdr.FileHeader.status.S_FLOAT==0
prec_str = ['int32=>',Dat.precision];
elseif hdr.FileHeader.status.S_32==0 ...
&& hdr.FileHeader.status.S_FLOAT==0
prec_str = ['int16=>',Dat.precision];
end
else
if hdr.BlockHeaders.status.S_FLOAT==1
prec_str = ['single=>',Dat.precision];
elseif hdr.BlockHeaders.status.S_32==1 ...
&& hdr.BlockHeaders.status.S_FLOAT==0
prec_str = ['int32=>',Dat.precision];
elseif hdr.BlockHeaders.status.S_32==0 ...
&& hdr.BlockHeaders.status.S_FLOAT==0
prec_str = ['int16=>',Dat.precision];
end
end
% Read k-space
tmp=fread(file_fid,...
[hdr.FileHeader.np,hdr.FileHeader.ntraces],...
prec_str);
% Store complex block and perform DC correction
if ~Dat.DCcorrection || ( nt(1)>1 )
complex_block = complex(tmp(1:2:end,:),tmp(2:2:end,:));
else
complex_block = complex(tmp(1:2:end,:)-hdr.BlockHeaders.lvl,...
tmp(2:2:end,:)-hdr.BlockHeaders.tlt);
end
%% Store and order k-space values
if any(Dat.AcqType==[1 2])
switch procpar.seqcon(2:3)
case {'cc','sc'}
kspace(:,:,ii) = complex_block;
otherwise
kspace(:,ii,:) = complex_block;
end
else
kspace(:,:,ii) = complex_block;
end
% Do not save blockheaders by default. When reading data files with a lot of
% blocks (e.g. over 1000) the processing of the DATA structure can be
% slowed down considerably. If you for some reason want to save also the
% block headers in the DATA structure comment out the code line below.
hdr.BlockHeaders = [];
end % for ii=1:hdr.
if Dat.ShowWaitbar
delete(wb_h)
end
else
%% -------------------------------------------------------------------
%% Fast Method for reading data. This may fail with some datas and can
%% also require a relatively large amount of memory. This
%% method should be used for EPI datas that contain a large number
%% of block headers...
% Check the size of the FID-file
d=dir(fopen(file_fid));
file_sz = d.bytes/1024/1024; % File size in MB
if file_sz<Dat.FileChunkSize
nBlocks = 1;
else
nBlocks = ceil(file_sz/Dat.FileChunkSize); % Read data in 500 MB blocks
end
% Initialize waitbar
if Dat.ShowWaitbar
if nBlocks==1
wb_h = aedes_calc_wait(['Reading ',num2str(Dat.AcqType),...
'D VNMR data (seqcon: "' procpar.seqcon '")']);
else
wb_h = aedes_wbar(1/nBlocks,{['Reading ',num2str(Dat.AcqType),...
'D VNMR data (seqcon: "' procpar.seqcon '")'],...
sprintf('Reading block 0/%d',nBlocks)});
end
end
% The first block header is read and it is assumed that the values in
% the other block headers don't change.
hdr.BlockHeaders.scale = fread(file_fid,1,'short'); % Scaling factor
tmp_status = fread(file_fid,1,'short'); % Status of data in block
hdr.BlockHeaders.status = [];
hdr.BlockHeaders.index = fread(file_fid,1,'short'); % Block index
tmp_mode = fread(file_fid,1,'short'); % Mode of data in block
hdr.BlockHeaders.mode = [];
hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
hdr.BlockHeaders.lpval = fread(file_fid,1,'float'); % f2 (2D-f1) left phase in phasefile
hdr.BlockHeaders.rpval = fread(file_fid,1,'float'); % f2 (2D-f1) right phase in phasefile
hdr.BlockHeaders.lvl = fread(file_fid,1,'float'); % level drift correction
hdr.BlockHeaders.tlt = fread(file_fid,1,'float'); % tilt drift correction
%% Parse status and mode bits
status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
for tt=1:length(BlockHeadStatusLabels)
if ~isempty(BlockHeadStatusLabels{tt})
hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
end
if ~isempty(BlockHeadModeLabels{tt})
hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
end
end
%% Determine data precision
if hdr.BlockHeaders.status.S_FLOAT==1
prec_str = ['single=>',Dat.precision];
prec = 4; % Precision in bytes
elseif hdr.BlockHeaders.status.S_32==1 ...
&& hdr.BlockHeaders.status.S_FLOAT==0
prec_str = ['int32=>',Dat.precision];
prec = 4;
elseif hdr.BlockHeaders.status.S_32==0 ...
&& hdr.BlockHeaders.status.S_FLOAT==0
prec_str = ['int16=>',Dat.precision];
prec = 2;
end
% Seek the file back to the beginning of the first block header
fseek(file_fid,-28,0);
% Determine the number of values that will result from block header(s)
nVals = (nbheaders*28)/prec;
nbh = floor(hdr.FileHeader.nblocks/nBlocks);
szh = nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces;
for ii=1:nBlocks
if nBlocks~=1
aedes_wbar(ii/nBlocks,wb_h,{['Reading ',num2str(Dat.AcqType),...
'D VNMR data (seqcon: "' procpar.seqcon '")'],...
sprintf('Reading block %d/%d',ii,nBlocks)});
end
% Read the whole data including block headers etc...
if ii==nBlocks
tmp = fread(file_fid,inf,prec_str);
else
tmp = fread(file_fid,nbh*szh,prec_str);
end
tmp=reshape(tmp,nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces,[]);
tmp(1:nVals,:)=[];
if ii==nBlocks
inds = ((ii-1)*nbh+1):size(kspace,2);
else
inds = ((ii-1)*nbh+1):ii*nbh;
end
% Do DC-correction if necessary
if ~Dat.DCcorrection || ( nt(1)>1 )
kspace(:,inds)=complex(tmp(1:2:end,:,:),tmp(2:2:end,:,:));
else
kspace(:,inds)=complex(tmp(1:2:end,:,:)-hdr.BlockHeaders.lvl,...
tmp(2:2:end,:,:)-hdr.BlockHeaders.tlt);
end
end
clear tmp
% Transform to 3D matrix
kspace = reshape(kspace,hdr.FileHeader.np/2,...
hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);
%% Store and order k-space values
if any(Dat.AcqType==[1 2])
switch procpar.seqcon(2:3)
case {'cc','sc'}
%kspace(:,:,ii) = complex_block;
otherwise
%kspace(:,ii,:) = complex_block;
kspace = permute(kspace,[1 3 2]);
end
else
%kspace(:,:,ii) = complex_block;
end
if Dat.ShowWaitbar
delete(wb_h)
pause(0.1)
end
end
% Remove singleton dimensions from kspace
kspace = squeeze(kspace);
% Check if raw kspace should be returned
if Dat.ReturnRawKspace
%% Delete waitbar
if Dat.ShowWaitbar && ishandle(wb_h)
delete(wb_h)
end
return
end
%% Support for RASER sequence ---------------------------
if isfield(procpar,'teType')
if strcmpi(procpar.sing_sh,'y') && strcmpi(procpar.teType,'c')
% Separate reference scans from the data
if isfield(procpar,'refscan') && strcmp(procpar.refscan,'y')
kspace=permute(reshape(kspace,procpar.np/2,procpar.ne/2,2,[]),[1 2 4 3]);
noise = kspace(:,:,:,1);
kspace = kspace(:,:,:,2);
% DC-correction
dc_level = squeeze(mean(mean(noise)));
for ii=1:size(kspace,3)
kspace(:,:,ii)=kspace(:,:,ii)-dc_level(ii);
end
end
if strcmpi(procpar.readType,'s')
sw = procpar.sw;
dwell = cumsum(ones(1,size(kspace,1))*(1/sw));
else
error('This feature has not been implemented yet!!!')
end
dwell = dwell*1000;
%kspace = permute(kspace,[2 1 3]);
% Phase correction
dims = size(kspace);
nI = size(kspace,3);
nv = size(kspace,4);
ne = size(kspace,2);
np = size(kspace,1);
g = 42.58;
G = procpar.gvox2;
Ds = procpar.vox2/100;
pos2 = procpar.pos2/10;
fudge = 0;
fudge2 = 1;
tt_dwell=repmat(dwell(:),1,ne);
tt_ne = repmat((1:ne),length(dwell),1);
tt_np = repmat((1:np)',1,ne);
i=sqrt(-1);
phi = (-2.*pi.*g.*G.*Ds.*tt_dwell.*(tt_ne./ne - 1/2 - pos2/Ds))+fudge.*tt_np;
phi = repmat(phi,[1,1,size(kspace,3)]);
kspace = abs(kspace).*exp(i*(angle(kspace)+phi));
% Flip every other column
sz(1)=size(kspace,1);
sz(2)=size(kspace,2);
sz(3)=size(kspace,3);
sz(4)=size(kspace,4);
kspace=reshape(kspace,sz(1),prod(sz(2:4)));
kspace(:,1:2:end) = flipud(kspace(:,1:2:end));
kspace = reshape(kspace,sz);
%kspace = reshape(permute(kspace,[2 1 3]),sz(2),[],1);
%kspace(:,1:2:end) = flipud(kspace(:,1:2:end));
%kspace=permute(reshape(kspace,sz(2),[],sz(3)),[2 1 3]);
else
data = [];
error('RASER sequence of this type has not been implemeted yet!')
end
% Reshape into 4D matrix
kspace = reshape(kspace,size(kspace,1),size(kspace,2),[],size(kspace,3));
% Reorient if requested
if Dat.OrientImages
kspace = flipdim(kspace,2);
end
data = abs(fftshift(fft(fftshift(fft(kspace,[],3),3),[],1),1));
% Delete kspace if not returned
if ~Dat.ReturnKSpace
kspace=[];
end
return
end
%% Support for fat/water measurements ----------------------
if isfield(procpar,'seqfil') && strcmpi(procpar.seqfil,'ge3d_csi2')
% Split kspace into fat and water (in 4th dimesion)
kspace=reshape(permute(reshape(kspace,256,2,[]),[1 3 4 2]),256,128,[],2);
% Fourier transform data
if any(Dat.ZeroPadding==[1 2])
data_sz = [procpar.np/2,procpar.nf,procpar.nv2,2];
data = zeros(data_sz,class(kspace));
data(:,:,:,1) = abs(fftshift(fftn(kspace(:,:,:,1),data_sz(1:3))));
data(:,:,:,2) = abs(fftshift(fftn(kspace(:,:,:,2),data_sz(1:3))));
else
data = zeros(size(kspace),class(kspace));
data(:,:,:,1) = abs(fftshift(fftn(kspace(:,:,:,1))));
data(:,:,:,2) = abs(fftshift(fftn(kspace(:,:,:,2))));
end
% Delete kspace if not returned
if ~Dat.ReturnKSpace
kspace=[];
end
return
end
% Check the number of receivers (PI stuff)
if isfield(procpar,'rcvrs') && ~isempty(procpar.rcvrs) && ...
length(procpar.rcvrs{1})>1
% Multiple receivers used
if isfield(procpar,'apptype') && strcmp(procpar.apptype{1},'imEPI')
% Store multiple receiver data in EPI measurements in 5th dimension
% and calculate sum-of-squares image
nRcv = length(find(procpar.rcvrs{1}=='y'));
if ~isfield(procpar,'epiref_type') || strcmpi(procpar.epiref_type,'single')
nRef = 1;
elseif strcmpi(procpar.epiref_type,'triple')
nRef = 3;
end
nVols = size(kspace,3)/nRcv-nRef;
if Dat.EPIPhasedArrayData
data = zeros(procpar.nv,procpar.np/2,procpar.ns,nVols+nRef,nRcv,'single');
else
data = zeros(procpar.nv,procpar.np/2,procpar.ns,nVols+nRef,'single');
end
kssz=size(kspace);
blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
nBlocks = ceil((size(kspace,3)/nRcv-nRef)/blksz);
lnum = length(num2str(nBlocks));
lnumstr = num2str(lnum);
bsl = lnum*2+1;
fprintf(1,'Processing data in blocks of %d volumes\n',blksz)
fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks);
for ii=1:nBlocks
fprintf(1,repmat('\b',1,bsl));
fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks);
tmp_data = [];
for kk=1:nRcv
inds_ref = kk:nRcv:nRcv*nRef;
inds_im = (nRcv*nRef+kk):nRcv:kssz(3);
inds = cat(2,inds_ref,inds_im(((ii-1)*blksz+1):min(ii*blksz,nVols)));
tmp_kspace = l_ReconstructKspace(kspace(:,:,inds),procpar,Dat);
tmp_data(:,:,:,:,kk) = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
end
if Dat.EPIPhasedArrayData
data_block = abs(tmp_data);
else
data_block = sqrt(sum(tmp_data.*conj(tmp_data),5));
end
if ii==1
data(:,:,:,1:size(tmp_data,4),:) = data_block;
elseif ii==nBlocks
data_block(:,:,:,1:nRef,:)=[];
data(:,:,:,((ii-1)*blksz+1+nRef):(nVols+nRef),:) = data_block;
else
data_block(:,:,:,1:nRef,:)=[];
data(:,:,:,((ii-1)*blksz+1+nRef):(ii*blksz+nRef),:) = data_block;
end
end
fprintf(1,'\n')
% Remove reference image if requested
if Dat.isEPIdata && Dat.RemoveEPIphaseIm
data(:,:,:,1:nRef,:)=[];
end
if Dat.OrientImages && ~isempty(procpar) && ...
isfield(procpar,'orient') && any(Dat.AcqType==[2 3 4])
orient = procpar.orient{1};
if any(strcmpi(orient,{'xyz','trans90','cor90','sag90'}))
data = flipdim(aedes_rot3d(data,1,3),2);
elseif strcmpi(orient,'oblique')
data = flipdim(flipdim(data,1),2);
else
data = flipdim(flipdim(data,1),2);
end
end
else
nRcv = length(find(procpar.rcvrs{1}=='y'));
data = [];
kspace2 = [];
for ii=1:nRcv
tmp_kspace = l_ReconstructKspace(kspace(:,:,ii),procpar,Dat);
kspace2(:,:,:,:,ii)=tmp_kspace;
if Dat.ReturnFTData
tmp_data = l_CalculateFFT(tmp_kspace,procpar,Dat);
data(:,:,:,:,ii) = tmp_data;
end
end
kspace = kspace2;
kspace2=[];
end
else
% Only one receiver
kspace = l_ReconstructKspace(kspace,procpar,Dat);
data = l_CalculateFFT(kspace,procpar,Dat);
end
% Delete kspace if not returned
if ~Dat.ReturnKSpace
kspace=[];
end
%===============================================
% Reconstruct kspace for different sequences
%===============================================
function kspace=l_ReconstructKspace(kspace,procpar,Dat)
%% Flip images for certain measurements
if Dat.FlipKspace~=0
if ischar(Dat.FlipInd)
if strcmpi(Dat.FlipInd,'all')
flipind = 1:size(kspace,3);
elseif strcmpi(Dat.FlipInd,'alt')
flipind = 2:2:size(kspace,3);
end
else
flipind = Dat.FlipInd;
end
for ii=flipind
if Dat.FlipKspace==1
kspace(:,:,ii)=fliplr(kspace(:,:,ii));
elseif Dat.FlipKspace==2
kspace(:,:,ii)=flipud(kspace(:,:,ii));
else
kspace(:,:,ii)=flipud(fliplr(kspace(:,:,ii)));
end
end
end
% Handle arrayed and RARE sequences
if Dat.isAcqArrayed && Dat.AcqType~=1 && Dat.Sorting && ~Dat.isEPIdata
%if ~strcmpi(procpar.seqcon(2:3),'cc') && ( isempty(Dat.phasetable) | ...
% isfield(procpar,'flash_converted') )
if (( isempty(Dat.phasetable) && ~strcmpi(procpar.seqcon(2:3),'sc') ) && ...
~(strcmpi(procpar.seqcon(2:3),'cc') && Dat.ArrayLength==size(kspace,3))) || ...
isfield(procpar,'flash_converted')
% Sort uncompressed arrayed data
ks_order = 1:procpar.nv*Dat.ArrayLength;
tmp = 0:procpar.nv:(procpar.ne-1)*procpar.nv;
tmp=repmat(tmp,procpar.nv*Dat.ArrayLength,1);
ks_order = reshape(ks_order,Dat.ArrayLength,procpar.nv).';
ks_order = ks_order(:);
ks_order = repmat(ks_order,1,procpar.ne);
ks_order = ks_order+tmp;
ks_order = ks_order.';
if procpar.ns==1
kspace=kspace(:,ks_order);
else
kspace = kspace(:,ks_order,:);
end
% Reshape into 3D matrix
kspace=reshape(kspace,[size(kspace,1) procpar.nv Dat.ArrayLength* ...
procpar.ns]);
else
%% Sort RARE type data
if Dat.UsePhaseTable && ~isempty(Dat.phasetable)
%Dat.phasetable = Dat.phasetable';
if Dat.isAcqArrayed && Dat.ArrayLength>1 && strcmpi(procpar.seqcon(2:3),'cs')
kspace = permute(reshape(reshape(kspace,size(kspace,1),[]),size(kspace,1),Dat.ArrayLength,[]),[1 3 2]);
else
Dat.phasetable = Dat.phasetable.';
end
kspace(:,Dat.phasetable(:),:) = kspace;
end
end
elseif Dat.isEPIdata
%% Support for EPI-measurements
% EPI data is measured with old VNMR
if isfield(procpar,'readres') && isfield(procpar,'phaseres')
% Number of slices
tmp_ns=length(procpar.pss);
if isfield(procpar,'navecho') && strcmpi(procpar.navecho{1},'y')
tmp_nv = procpar.nv-procpar.nseg;
else
tmp_nv = procpar.nv;
end
kspace = reshape(kspace,[size(kspace,1) ...
size(kspace,2)/tmp_nv tmp_nv]);
kspace = permute(kspace,[1 3 2]);
% Reshape to 4D matrix
kspace = reshape(kspace,[size(kspace,1) size(kspace,2) ...
tmp_ns size(kspace,3)/tmp_ns]);
elseif isfield(procpar,'apptype') && strcmp(procpar.apptype{1},'imEPI')
% If EPI data is measured with new VNMRj system
% Number of slices
ns = length(procpar.pss);
% Reshape kspace to 4D matrix
if isfield(procpar,'fract_ky') && procpar.fract_ky~=procpar.nv/2
kspace = reshape(kspace,procpar.np/2,procpar.nv/2+procpar.fract_ky,...
ns,[]);
if Dat.ZeroPadding~=0
kspace(procpar.np/2,procpar.nv,1)=0;
end
else
kspace = reshape(kspace,procpar.np/2,procpar.nv,...
ns,[]);
end
% Flip even kspace lines
if ~Dat.FastDataRead
for tt=2:2:size(kspace,2)
kspace(:,tt,:,:) = flipdim(kspace(:,tt,:,:),1);
end
else
kspace(:,2:2:end,:,:) = flipdim(kspace(:,2:2:end,:,:),1);
end
% Sort centric measurements
if isfield(procpar,'ky_order') && strcmpi(procpar.ky_order,'c')
kspace(:,Dat.phasetable(:),:,:)=kspace;
end
%kspace = flipdim(kspace,1);
% EPI phase correction -------------------------
if ~isfield(procpar,'epiref_type') || strcmpi(procpar.epiref_type,'single')
% Single reference pointwise phase correction
% Get the reference image
ref_im = kspace(:,:,:,1);
% Do phase correction.
phase_e = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
for kk=2:size(kspace,4)
kspace(:,:,:,kk) = ifft(fft(kspace(:,:,:,kk),[],1).*phase_e,[],1);
end
elseif strcmpi(procpar.epiref_type,'triple')
% Triple reference pointwise phase correction
% Get the reference images
ref1 = kspace(:,:,:,1);
phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
ref2 = flipdim(kspace(:,:,:,3),1);
phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
im1 = flipdim(kspace(:,:,:,2),1);
% Correct phase for reversed read gradient image
rev_phase = fft(im1,[],1).*phase2;
for kk=4:size(kspace,4)
kspace(:,:,:,kk)=ifft(rev_phase+fft(kspace(:,:,:,kk),[],1).*phase1);
end
end
end
else
% This section contains various "chewing gum and ironwire" type
% patches that hopefully work...
if strcmp(procpar.seqcon,'nncnn')
kspace = permute(kspace,[1 3 2]);
elseif strcmp(procpar.seqcon,'nccnn') && length(size(kspace))==2 && ...
procpar.ns>=1 && Dat.AcqType~=1
if ~isempty(Dat.phasetable)
kssz = size(kspace);
phsz = size(Dat.phasetable);
kspace=permute(reshape(kspace,procpar.np/2,...
phsz(2),...
kssz(2)/phsz(2)),[1 3 2]);
kspace = permute(reshape(kspace,...
procpar.np/2,procpar.ns,kssz(2)/procpar.ns),[1 3 2]);
if Dat.Sorting
kspace(:,Dat.phasetable(:),:)=kspace;
end
else
kspace = reshape(kspace,...
[procpar.np/2 procpar.ns ...
size(kspace,2)/procpar.ns]);
kspace=permute(kspace,[1 3 2]);
end
elseif strcmp(procpar.seqcon,'nscsn') && length(size(kspace))==3 && ...
Dat.AcqType~=1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Support for 3D fast spin-echo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(Dat.phasetable)
for ii=1:size(kspace,3)
tmp = kspace(:,:,ii);
kssz = size(tmp);
phsz = size(Dat.phasetable);
tmp=permute(reshape(tmp,procpar.np/2,...
phsz(2),...
kssz(2)/phsz(2)),[1 3 2]);
tmp = permute(reshape(tmp,...
procpar.np/2,procpar.ns,kssz(2)/procpar.ns),[1 3 2]);
if Dat.Sorting
tmp(:,Dat.phasetable(:),:)=tmp;
end
kspace(:,:,ii)=tmp;
end
end
elseif Dat.AcqType~=1 && isfield(procpar,'nv')
if length(size(kspace))==2 && ...
( size(kspace,1)~=(procpar.np/2) || size(kspace,2)~=procpar.nv )
%% Reshape the kspace to the appropriate size
kspace=reshape(kspace,[procpar.np/2 procpar.nv size(kspace,2)/ ...
procpar.nv]);
elseif length(size(kspace))==3
kspace = reshape(kspace,procpar.np/2,procpar.nv,[]);
if Dat.Sorting && ~isempty(Dat.phasetable)
kspace(:,Dat.phasetable(:),:) = kspace;
end
end
end
%%% Support for Teemu's ASE3D-data
%if strcmpi(procpar.seqcon,'ncccn')
% %% Reshape the kspace to the appropriate size
% kspace=reshape(kspace,[procpar.np/2 procpar.nv size(kspace,2)/procpar.nv]);
%end
end
%=========================================
% Fourier Transform Data
%=========================================
function data=l_CalculateFFT(kspace,procpar,Dat)
data=[];
% Return image/spectral data --------------------------------
if Dat.ReturnFTData
%% Fourier transform spectral data
if Dat.AcqType==1
wb_h = aedes_wbar(0,'Fourier transforming spectral data');
%data=zeros(size(kspace),'single');
data=kspace;
sz=size(kspace,2)*size(kspace,3);
count=1;
for kk=1:size(kspace,3)
for ii=1:size(kspace,2)
%% Update waitbar
if Dat.ShowWaitbar
aedes_wbar(count/sz,...
wb_h,...
{'Fourier transforming spectral data',...
['(Processing spectrum ' num2str(count) '/' num2str(sz) ')']})
end
data(:,ii,kk) = abs(fftshift(fft(data(:,ii,kk))));
count=count+1;
end
end
%% Fourier transform image data
else
%% Zeropadding auto
if Dat.ZeroPadding==2
if Dat.isEPIdata && isfield(procpar,'readres') && isfield(procpar,'phaseres')
%% Determine the image size for EPI data for procpar "readres"
%% and "phaseres" fields
data_sz = [procpar.readres procpar.phaseres];
if isequal(data_sz,[size(kspace,1),size(kspace,2)])
DoZeroPadding=false;
else
DoZeroPadding=true;
end
data_sz = [procpar.readres ...
procpar.phaseres size(kspace,3) size(kspace,4)];
else
%% Check if zeropadding is necessary.
ks_sz_sorted = sort([size(kspace,1),size(kspace,2)]);
fov_sz_sorted = sort([procpar.lro,procpar.lpe]);
sliceRelativeDim = ks_sz_sorted(2)/ks_sz_sorted(1);
FOVrelativeDim = fov_sz_sorted(2)/fov_sz_sorted(1);
if FOVrelativeDim~=sliceRelativeDim
ind=find([size(kspace,1),size(kspace,2)]==ks_sz_sorted(1));
data_sz = size(kspace);
data_sz(ind) = round((fov_sz_sorted(1)/fov_sz_sorted(2))*ks_sz_sorted(2));
DoZeroPadding=true;
else
data_sz = size(kspace);
DoZeroPadding=false;
end
end
%% Force zeropadding on. Force images to be square.
elseif Dat.ZeroPadding==1
ks_sz_max = max(size(kspace,1),size(kspace,2));
data_sz = size(kspace);
data_sz(1:2)=ks_sz_max;
DoZeroPadding=true;
%% Force zeropadding off
else
data_sz = size(kspace);
DoZeroPadding=false;
end
%% Fourier transform 2D and EPI image data
if Dat.AcqType==2 || Dat.isEPIdata
sz=size(kspace,3);
if DoZeroPadding
data=zeros(data_sz,class(kspace));
data(1:size(kspace,1),1:size(kspace,2),:,:)=kspace;
else
data=kspace;
end
if ~Dat.ReturnKSpace
kspace = [];
end
if Dat.ShowWaitbar
[wb_h,txh]=aedes_calc_wait('Fourier transforming image data...');
end
%data=permute(fftshift(fftshift(abs(fft(permute(fft(data),[2 1 3 4]))),1),2),[2 1 3 4]);
if ~Dat.FastDataRead
for tt=1:size(data,4)
data(:,:,:,tt)=fftshift(fftshift(abs(fft(fft(data(:,:,:,tt),[],1),[],2)),1),2);
end
else
data=fftshift(fftshift(abs(fft(fft(data,[],1),[],2)),1),2);
end
else
%% Fourier transform 3D image data
if Dat.ShowWaitbar
[wb_h,txh]=aedes_calc_wait('Fourier transforming 3D image data...');
end
if DoZeroPadding
% Check zeropadding in the 3rd dimension if zeropadding is set to
% "auto"
if Dat.ZeroPadding==2
if isfield(procpar,'lpe2') && isfield(procpar,'lpe')
data_sz(3) = max(1,round((procpar.lpe2/procpar.lpe)*data_sz(2)*Dat.ArrayLength));
end
end
data=zeros(data_sz,class(kspace));
data = abs(fftshift(fftn(kspace,data_sz)));
else
data=zeros(data_sz,class(kspace));
data = abs(fftshift(fftn(kspace)));
end
end
% Remove reference image if requested
if Dat.isEPIdata && Dat.RemoveEPIphaseIm
data = data(:,:,:,2:end);
end
% Reorient images and remove the reference image if requested
if Dat.OrientImages && ~isempty(procpar) && ...
isfield(procpar,'orient') && any(Dat.AcqType==[2 3 4])
orient = procpar.orient{1};
if any(strcmpi(orient,{'xyz','trans90','cor90','sag90'}))
data = flipdim(aedes_rot3d(data,1,3),2);
elseif strcmpi(orient,'oblique')
data = flipdim(flipdim(data,1),2);
else
data = flipdim(flipdim(data,1),2);
end
end
% $$$ %% Fourier transform 4D image data
% $$$ elseif Dat.AcqType==4
% $$$ data = abs(fftshift(fftshift(fft2(kspace),2),1));
end
%% Delete waitbar
if Dat.ShowWaitbar && ishandle(wb_h)
delete(wb_h)
end
end
%========================================================
% A function for creating phasetables for fse and fsems
%========================================================
function phasetable=l_CreateFsemsPhaseTable(procpar)
phasetable = [];
if ~isfield(procpar,'etl') || ~isfield(procpar,'kzero')
return
end
etl = procpar.etl;
kzero = procpar.kzero;
nv = procpar.nv;
t = (-nv/2+1):(nv/2);
t = t(:);
tt = flipud(t);
tt = tt(:);
phasetable = [reshape(t(1:nv/2),[],etl);...
flipud(reshape(tt(1:nv/2),[],etl))];
phasetable = circshift(fliplr(phasetable),[0 kzero-1]);
%% - EOF - %%
return