source: aedes_readfid.m @ 103

Last change on this file since 103 was 103, checked in by tjniskan, 9 years ago
  • Added support for multiple receivers in aedes_readfid.m
  • Aedes can now read data stored in MAT-file from a variable that starts

with "data" (not case sensitive)

M aedes_data_read.m
M aedes_readfid.m
M aedes_revision.m

File size: 52.8 KB
RevLine 
[80]1function [DATA,msg_out] = aedes_readfid(filename,varargin)
2% AEDES_READFID - Read VNMR (Varian) FID-files
[2]3%   
4%
5% Synopsis:
[80]6%        [DATA,msg]=aedes_readfid(filename,'PropertyName1',value1,'PropertyName2',value2,...)
7%        [DATA,msg]=aedes_readfid(filename,'header')
8%        [DATA,msg]=aedes_readfid(filename,output_filename)
[2]9%
10% Description:
11%        The function reads VNMR (Varian) FID-file and return a structure
12%        DATA with fields DataFormat, HDR, FTDATA, KSPACE, PROCPAR, and
13%        PHASETABLE. The fields of the DATA structure are constructed as
14%        follows:
15%       
16%        DATA
17%          |-> DataFormat         (string identifier for data format 'vnmr')
18%          |-> HDR
19%                |-> FileHeader   (data file header)
20%                |-> BlockHeaders (data block headers, not stored by default)
21%                |-> fname        (file name)
22%                |-> fpath        (file path)
[80]23%                |-> Param        (parameter values used by AEDES_READFID to read the data)
[2]24%          |-> FTDATA             (real fourier transformed image data)
25%          |-> KSPACE             (complex k-space, empty by default)
26%          |-> PROCPAR            (parameters from procpar file)
27%          |-> PHASETABLE         (phasetable)
28%
29%        The DATA structure is returned as empty (without HDR, FTDATA,
30%        KSPACE, and PROCPAR fields) if an error has occurred during
31%        reading. The error message is returned in the second output
[80]32%        argument msg. If AEDES_READFID is called with only one output argument
[2]33%        (i.e. without MSG), the error message is echoed to the workspace.
34%       
35%        The first input argument is either a string containing full path to
36%        the FID-file or the header structure HDR. Only the data file header
37%        can be read by giving a string 'header' as the second input
38%        argument.
39%
40%        By default the k-space is not returned, i.e. the field KSPACE is
41%        empty. The returned data can be adjusted by using the 'return'
42%        property and values 1, 2, or 3 (see below for more information).
43%
[80]44%        The supported property-value pairs in AEDES_READFID (property strings
[2]45%        are not case sensitive):
46%
47%        Property:        Value:                Description:
48%        *********        ******                ************
49%        'Return'      :  [ {1} | 2 | 3 | 4 ]   % 1=return only ftdata (default)
50%                                               % 2=return only k-space
51%                                               % 3=return both ftdata & kspace
52%                                               % 4=return raw kspace
53%
[21]54%        'FastRead'    :  [{'off'} | 'on' ]     % Enables an experimental
55%                                               % method for very fast reading
56%                                               % of fid-files. This not as
57%                                               % robust as the older
58%                                               % method and can consume a lot
[22]59%                                               % of memory.
[21]60%
[2]61%        'OutputFile'  :  filename              % Writes the slices into
62%                                               % NIfTI-format files
63%                                               % using the given filename.
64%
65%        'DCcorrection':  [ {'off'} | 'on' ]    % 'on'=use DC correction
66%                                               % 'off'=don't use DC correction
67%                                               % (default)
68%
69%        'Procpar'     :  (procpar-structure)   % Input procpar
70%                                               % structure. If this
71%                                               % property is not set the
72%                                               % procpar structure
73%                                               % will be read using
[80]74%                                               % AEDES_READPROCPAR
[2]75%
76%        'ZeroPadding' :  ['off'|'on'|{'auto'}] % 'off' = force
77%                                               % zeropadding off
78%                                               % 'on' = force
79%                                               % zeropadding on (force
80%                                               % images to be square)
81%                                               % 'auto' = autoselect
82%                                               % using relative FOV
83%                                               % dimensions (default)
84%
85%        'PhaseTable'  : (custom phase table)   % User-specified
86%                                               % line order for
87%                                               % k-space. The phase table
88%                                               % is obtained from the file
89%                                               % specified in
90%                                               % procpar.petable if
91%                                               % necessary.
92%
93%        'sorting'      : [ 'off' | {'on'} ]    % 'off' =disable k-space
94%                                               % sorting, i.e. ignore the
95%                                               % phase table and/or arrays.
96%                                               % 'on' =sort k-space using
97%                                               % phase table and/or array
98%                                               % information if necessary
99%                                               % (default)
100%
101%        'wbar'        : [ {'on'} | 'off' ]     % 'on'  = show waitbar
102%                                               % 'off' = don't show waitbar
103%
104%        'FlipKspace'  : [ {'off'} | 'LR' |
105%                             'UD' | 'LRUD' ]   % 'off' = don't flip (default)
106%                                               % 'LR' = left/right
107%                                               % 'UD' = up/down
108%                                               % 'LRUD' = left/right and
109%                                               % up/down
110%
111%        'FlipInd'     : [ {'all'} | 'alt' |
112%                          (custom vector)  ]   % 'all' = flip all slices
113%                                               % (default)
114%                                               % 'alt' = flip every
115%                                               % second slice
116%                                               % (custom vector) = A
117%                                               % custom vector containing
118%                                               % indices to the flipped slices
119%
[36]120%        'Precision'   : ['single'|{'double'}]  % Precision of the
121%                                               % outputted data.
122%                                               % 'single' = 32-bit float
123%                                               % 'double' = 64-bit float
124%
[75]125%        'OrientImages': [ {'on'} | 'off' ]     % Orient FTDATA after
126%                                               % reading the data using
127%                                               % PROCPAR.orient property.
[67]128%
129%        'RemoveEPIphaseIm' : ['on'|{'off'}]    % Remove the phase image
130%                                               % from EPI data. This
131%                                               % option only affects EPI
132%                                               % images.
133%
[2]134% Examples:
[80]135%        [DATA,msg]=aedes_readfid(filename)             % Read image data from 'filename'
136%        [DATA,msg]=aedes_readfid(filename,'header')    % Read only data file header
[2]137%
[80]138%        [DATA,msg]=aedes_readfid(filename,'return',1)  % Return only image data (default)
139%        [DATA,msg]=aedes_readfid(filename,'return',2)  % Return only k-space
140%        [DATA,msg]=aedes_readfid(filename,'return',3)  % Return both image data and k-space
[2]141%       
142%        % Read first data file header and then image data and k-space
[80]143%        [DATA,msg]=aedes_readfid(filename,'header')
144%        [DATA,msg]=aedes_readfid(DATA.HDR,'return',3)
[2]145%
146%        % Read VNMR-data with default options and save slices in NIfTI
147%        % format
[80]148%        DATA=aedes_readfid('example.fid','saved_slices.nii');
[2]149%
150%        % If you want to use non-default options and also write
151%        % NIfTI-files, use the property/value formalism, for example
[80]152%        DATA=aedes_readfid('example.fid','ZeroPadding','off',...
[2]153%                     'OutputFile','saved_slices.nii');
154%
155% See also:
[80]156%        AEDES_READFIDPREFS, AEDES_READPROCPAR, AEDES_READPHASETABLE,
157%        AEDES_DATA_READ, AEDES_WRITE_NIFTI, AEDES
[2]158
[37]159% This function is a part of Aedes - A graphical tool for analyzing
[36]160% medical images
[2]161%
[36]162% Copyright (C) 2006 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uku.fi>
163%
[45]164% Department of Physics, Department of Neurobiology
[39]165% University of Kuopio, FINLAND
[36]166%
167% This program may be used under the terms of the GNU General Public
168% License version 2.0 as published by the Free Software Foundation
169% and appearing in the file LICENSE.TXT included in the packaging of
170% this program.
171%
172% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
173% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
[2]174
175
[36]176
177
[2]178if nargout==2
179  Dat.ShowErrors = false;
180  msg_out='';
181else
182  Dat.ShowErrors = true;
183end
184
[35]185%% ---------------------
186% Defaults
187Dat.ReturnKSpace  = false;
188Dat.ReturnFTData  = true;
189Dat.DCcorrection  = false;
190Dat.ZeroPadding = 2;
191Dat.Sorting = true;
192Dat.UsePhaseTable = true;
[77]193Dat.FastDataRead = true;
194Dat.precision = 'single';
[2]195
[67]196
[2]197%% Other defaults
198Dat.ShowWaitbar   = true;
199procpar=[];
200Dat.phasetable=[];
201Dat.FlipKspace = 0;
202Dat.FlipInd = 'all';
203Dat.OutputFile = false;
204Dat.ReturnRawKspace = false;
[67]205Dat.ReorientEPI = false;
206Dat.RemoveEPIphaseIm = false;
[75]207Dat.OrientImages = true;
[2]208
209%% -------------------------------------------------------------------
210
211
212%% Set data format label
213DATA.DataFormat = 'vnmr';
214
215%% Parse input arguments
216if nargin==0 || isempty(filename)
217 
218  %% Get default directory
219  try
[37]220    defaultdir = getpref('Aedes','GetDataFileDir');
[2]221  catch
222    defaultdir = [pwd,filesep];
223  end
224 
225  %% If no input arguments are given, ask for a file
226  [f_name, f_path, f_index] = uigetfile( ...
227       {'fid;FID','Varian FID-files (fid)'; ...
228        '*.*', 'All Files (*.*)'}, ...
229        'Select VNMR (Varian) FID file',defaultdir);
[22]230  if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed
[2]231    DATA=[];
232    msg_out='Canceled';
233    return
234  end
235  ReadHdr = true;
236  ReadData = true;
237  filename=[f_path,f_name];
238 
239  %% Set default directory to preferences
[37]240  setpref('Aedes','GetDataFileDir',f_path)
[2]241 
242end
243
244if nargin==1
245  if isstruct(filename)
246    hdr = filename;
247    filename = [hdr.fpath,hdr.fname];
248    ReadHdr = false;
249    ReadData = true;
250% $$$     ReturnKSpace = false;
251% $$$     ReturnFTData = true;
252   
253  elseif ischar(filename)
254    ReadHdr = true;
255    ReadData = true;
256% $$$     ReturnKSpace = false;
257% $$$     ReturnFTData = true;
258  end
259elseif nargin==2
260  if strcmpi(varargin{1},'header')
261    ReadHdr = true;
262    ReadData = false;
263% $$$     ReturnKSpace = false;
264  elseif ischar(varargin{1})
265    ReadHdr = true;
266    ReadData = true;
267    Dat.OutputFile = varargin{1};
268  else
269    if ~Dat.ShowErrors
270      DATA=[];
271      msg_out=sprintf('Invalid second input argument "%s".',varargin{1});
272      return
273    else
274      error('Invalid second input argument "%s".',varargin{1})
275    end
276  end
277else
278  if isstruct(filename)
279    hdr = filename;
280    filename = [hdr.fpath,hdr.fname];
281    ReadHdr = false;
282    ReadData = true;
283  elseif isempty(filename)
284    [f_name, f_path, f_index] = uigetfile( ...
285        {'fid;FID','Varian FID-files (fid)'; ...
286         '*.*', 'All Files (*.*)'}, ...
287        'Select VNMR (Varian) FID file');
[22]288    if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed
[2]289      DATA=[];
290      msg_out='Canceled';
291      return
292    end
293    ReadHdr = true;
294    ReadData = true;
295    filename=[f_path,f_name];
296  else
297    ReadHdr = true;
298    ReadData = true;
299  end
300 
301  for ii=1:2:length(varargin)
302    switch lower(varargin{ii})
303     case 'return'
304      if length(varargin)<2
305        DATA=[];
306        if ~Dat.ShowErrors
307          msg_out='"Return" value not specified!';
308        else
309          error('"Return" value not specified!')
310        end
311        return
312      else
313        if varargin{ii+1}==1
314          Dat.ReturnKSpace = false;
315          Dat.ReturnFTData = true;
316        elseif varargin{ii+1}==2
317          Dat.ReturnKSpace = true;
318          Dat.ReturnFTData = false;
319        elseif varargin{ii+1}==3
320          Dat.ReturnKSpace = true;
321          Dat.ReturnFTData = true;
322                elseif varargin{ii+1}==4
323                  Dat.ReturnRawKspace = true;
324                  Dat.ReturnKSpace = true;
325          Dat.ReturnFTData = false;
326        end
327      end
328     
329     case {'dc','dccorrection','dccorr'}
330      if strcmpi(varargin{ii+1},'on')
331        Dat.DCcorrection = true;
332      else
333        Dat.DCcorrection = false;
334      end
335     
336     case 'procpar'
337      procpar=varargin{ii+1};
338     
339     case 'zeropadding'
[38]340      if ischar(varargin{ii+1})
[2]341        if strcmpi(varargin{ii+1},'on')
342          Dat.ZeroPadding = 1; % on
343        elseif strcmpi(varargin{ii+1},'off')
344          Dat.ZeroPadding = 0; % off
345        else
346          Dat.ZeroPadding = 2; % auto
347        end
348      else
349        % Undocumented
350        Dat.ZeroPadding = 3; % Custom
351        Dat.CustomZeroPadding = varargin{ii+1};
352      end
353     
354     case 'phasetable'
355      Dat.phasetable = varargin{ii+1};
356     
357     case 'sorting'
[22]358          if strcmpi(varargin{ii+1},'on')
[2]359        Dat.UsePhaseTable = true;
360        Dat.Sorting = true;
361      else
362        Dat.UsePhaseTable = false;
363        Dat.Sorting = false;
[21]364          end
365     
366         case 'fastread'
367           if strcmpi(varargin{ii+1},'on')
[22]368                 Dat.FastDataRead = true;
369           else
370                 Dat.FastDataRead = false;
371           end
[21]372           
[2]373     case 'wbar'
374      if strcmpi(varargin{ii+1},'on')
375        Dat.ShowWaitbar = 1;
376      else
377        Dat.ShowWaitbar = 0;
[36]378          end
379         
380         case 'precision'
381           if strcmpi(varargin{ii+1},'single')
382                 Dat.precision = 'single';
383           end
384               
385          case 'flipkspace'
[2]386       
387        flipkspace = varargin{ii+1};
388        if strcmpi(flipkspace,'off')
389          Dat.FlipKspace=0;
390        elseif strcmpi(flipkspace,'LR')
391          Dat.FlipKspace=1;
392        elseif strcmpi(flipkspace,'UD')
393          Dat.FlipKspace=2;
394        elseif strcmpi(flipkspace,'LRUD')
395          Dat.FlipKspace=3;
396        else
397          DATA=[];
398          if ~Dat.ShowErrors
[22]399            msg_out = 'Bad value for property FlipKspace!';
[2]400          else
401            error('Bad value for property FlipKspace!')
402          end
403          return
404        end
405       
406      case 'flipind'
407        Dat.FlipInd = varargin{ii+1};
408       
409      case 'outputfile'
410        Dat.OutputFile = varargin{ii+1};
411       
[67]412      case 'reorientepi'
413        if strcmpi(varargin{ii+1},'on')
414          Dat.ReorientEPI = true;
415        end
416       
417      case 'removeepiphaseim'
418        if strcmpi(varargin{ii+1},'on')
419          Dat.RemoveEPIphaseIm = true;
420        end
[75]421      case 'orientimages'
422        if strcmpi(varargin{ii+1},'off')
423          Dat.OrientImages = false;
424        end
[2]425     otherwise
426      DATA=[];
427      if ~Dat.ShowErrors
428        msg_out = ['Unknown property "' varargin{ii} '"'];
429      else
430        error(['Unknown property "' varargin{ii} '"'])
431      end
432      return
433    end
434  end
435end
436
437% Parse filename
438[fpath,fname,fext]=fileparts(filename);
439if ~strcmpi(fname,'fid')
440  if isempty(fname)
441    fpath = [fpath,filesep];
442  else
[53]443        if isempty(fpath)
444          fpath = [pwd,filesep,fname,fext,filesep];
445        else
446          fpath = [fpath,filesep,fname,fext,filesep];
447        end
[2]448  end
449  fname = 'fid';
450else
451  fpath = [fpath,filesep];
452end
453
454%% Read procpar if not given as input argument
455if isempty(procpar) || nargin~=2
[80]456  [procpar,proc_msg]=aedes_readprocpar([fpath,'procpar']);
[2]457  if isempty(procpar)
458    DATA=[];
459    if ~Dat.ShowErrors
460      msg_out=proc_msg;
461    else
462      error(proc_msg)
463    end
464    return
465  end
466end
467
468%% Read phasetable if necessary
469if isfield(procpar,'petable') && isempty(Dat.phasetable) && ...
[71]470    ~isempty(procpar.petable{1}) && ~strcmpi(procpar.petable{1},'n') && ...
471    Dat.Sorting
[2]472  % Look in preferences for tablib-directory
473  try
[37]474    tabpath = getpref('Aedes','TabLibDir');
[2]475  catch
[37]476    % If the table path was not found in preferences, try to use aedes
[2]477    % directory
[37]478    tmp_path = which('aedes');
[2]479    if ~isempty(tmp_path)
[37]480      aedes_path=fileparts(tmp_path);
481      tabpath = [aedes_path,filesep,'tablib',filesep];
[2]482    else
483      % If all else fails, look in the current directory
484      tabpath = [pwd,filesep,'tablib',filesep];
485    end
486  end
[80]487  [phasetable,msg] = aedes_readphasetable([tabpath,procpar.petable{1}]);
[2]488 
[71]489  % If petable cannot be found, check if it is in procpar...
490  if isempty(phasetable) && isfield(procpar,'pe_table')
491    phasetable = {'t1',procpar.pe_table(:)};
492  end
493 
[82]494  % If the sequence is a fast spin echo, try to construct the phasetable
495  if isempty(phasetable) && isfield(procpar,'petable') && ...
496      strncmpi(procpar.petable{1},'fse',3)
497    phasetable = {'t1',l_CreateFsemsPhaseTable(procpar)};
498  end
499 
[2]500  %% Abort and throw error if phasetable cannot be read and the FID-file
501  % has not been sorted
502  if isempty(phasetable) && ~isfield(procpar,'flash_converted')
503    DATA=[];
504    if ~Dat.ShowErrors
505      msg_out={['Could not find the required phase table "' ...
506                procpar.petable{1} '" in the following folders'],...
507               ['"' tabpath '".']};
508    else
509      error('Could not find the required phase table "%s" in %s',...
510            procpar.petable{1},tabpath)
511    end
512    return
[22]513  elseif ( isempty(phasetable) && isfield(procpar,'flash_converted') ) || ...
[2]514      ~Dat.UsePhaseTable
515    %% If the FID-file has been sorted, the reading can continue but
516    % throw a warning.
[80]517    fprintf(1,'Warning: aedes_readfid: Could not find phasetable "%s" in "%s"!\n',procpar.petable{1},tabpath)
[2]518    Dat.phasetable=[];
519  else
520    Dat.phasetable = phasetable{1,2};
521  end
522end
523
524% Convert phasetable indices to MATLAB indices
525if ~isempty(Dat.phasetable)
526  Dat.phasetable=Dat.phasetable+(-min(min(Dat.phasetable))+1);
527else
528  Dat.UsePhaseTable = false;
529end
530
531%% Open FID-file
532[file_fid,msg]=fopen([fpath,fname],'r','ieee-be');
533if file_fid < 0
534  DATA=[];
535  if ~Dat.ShowErrors
536    msg_out=['Could not open file "' filename '" for reading.'];
537  else
538    error('Could not open file "%s" for reading.',filename)
539  end
540  return
541end
542
543% Read only header
544if ReadHdr && ~ReadData
545  [hdr,msg_out]=l_ReadDataFileHeader(file_fid);
546  if ~isempty(msg_out)
547    DATA=[];
548    fclose(file_fid);
549    if Dat.ShowErrors
550      error(msg_out)
551    end
552    return
553  else
554    DATA.HDR.FileHeader = hdr.FileHeader;
555    DATA.FTDATA=[];
556    DATA.KSPACE=[];
557    DATA.PROCPAR=[];
558    DATA.PHASETABLE=[];
559  end
560elseif ~ReadHdr && ReadData % Header structure given as input argument
561                            % Read only data.
562  [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar);
563  if ~isempty(msg_out)
564    DATA=[];
565    fclose(file_fid);
566    if Dat.ShowErrors
567      error(msg_out)
568    end
569    return
570  else
571    DATA.HDR.FileHeader=hdr.FileHeader;
572    DATA.HDR.BlockHeaders = hdr.BlockHeaders;
573    DATA.FTDATA=data;
574    DATA.KSPACE=kspace;
575    DATA.PROCPAR=procpar;
576    DATA.PHASETABLE=Dat.phasetable;
577  end
578elseif ReadHdr && ReadData  % Read headers and data
579  [hdr,msg_out]=l_ReadDataFileHeader(file_fid);
580  [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar);
581  if ~isempty(msg_out)
582    DATA=[];
583    fclose(file_fid);
584    if Dat.ShowErrors
585      error(msg_out)
586    end
587    return
588  else
589    DATA.HDR.FileHeader=hdr.FileHeader;
590    DATA.HDR.BlockHeaders = hdr.BlockHeaders;
591    DATA.FTDATA=data;
592    DATA.KSPACE=kspace;
593    DATA.PROCPAR=procpar;
594    DATA.PHASETABLE=Dat.phasetable;
595  end
596end
597
598% Set file name and path to the HDR structure
599DATA.HDR.fname=fname;
600DATA.HDR.fpath=fpath;
601
602% Set parameter values
603DATA.HDR.Param.ReturnKSpace = Dat.ReturnKSpace;
604DATA.HDR.Param.ReturnFTData = Dat.ReturnFTData;
605if Dat.ZeroPadding==0
606  DATA.HDR.Param.ZeroPadding = 'off';
607elseif Dat.ZeroPadding==1
608  DATA.HDR.Param.ZeroPadding = 'on';
609else
610  DATA.HDR.Param.ZeroPadding = 'auto';
611end
612if ~Dat.DCcorrection
613  DATA.HDR.Param.DCcorrection = 'off';
614else
615  DATA.HDR.Param.DCcorrection = 'on';
616end
617if Dat.Sorting==0
618  DATA.HDR.Param.Sorting = 'off';
619else
620  DATA.HDR.Param.Sorting = 'on';
621end
622if Dat.FlipKspace==0
623  DATA.HDR.Param.FlipKspace = 'off';
624elseif Dat.FlipKspace==1
625  DATA.HDR.Param.FlipKspace = 'LR';
626elseif Dat.FlipKspace==2
627  DATA.HDR.Param.FlipKspace = 'UD';
628elseif Dat.FlipKspace==3
629  DATA.HDR.Param.FlipKspace = 'LRUD';
630end
631DATA.HDR.Param.FlipInd = Dat.FlipInd;
632
633
634% Close file
635fclose(file_fid);
636
637%% Write slices to NIfTI files
638if ischar(Dat.OutputFile)
639  if isempty(Dat.OutputFile)
640    [fp,fn,fe]=fileparts(DATA.HDR.fpath(1:end-1));
641    fprefix=fn;
642    niipath = [pwd,filesep];
643  else
644    [fp,fn,fe]=fileparts(Dat.OutputFile);
645    if strcmpi(fe,'.nii')
646      fprefix=fn;
647    else
648      fprefix=[fn,fe];
649    end
650    if isempty(fp)
651      niipath = [pwd,filesep];
652    else
653      niipath = [fp,filesep];
654    end
655  end
656 
657  % Create file names
658  filenames={};
659  for ii=1:size(DATA.FTDATA,3)
660    filenames{ii}=sprintf('%s%03d%s',[fprefix,'_'],ii,'.nii');
661  end
662  nFiles = length(filenames);
663   
[80]664  h=aedes_wbar(0,sprintf('Saving slice 1/%d in NIfTI format...',nFiles));
[2]665  for ii=1:nFiles
[80]666    aedes_wbar(ii/nFiles,h,sprintf('Saving slice %d/%d in NIfTI format...',ii, ...
[2]667                             nFiles))
[80]668    [done,msg]=aedes_write_nifti(DATA.FTDATA(:,:,ii),...
[2]669                           [niipath,filenames{ii}],'DataType','single',...
670                           'FileType',2);
671  end
672  delete(h)
673 
674  if ~done
675    warning('Error occurred while writing NIfTI-files. Could not write file(s)!')
676  end
677 
678end
679return
680
681%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
682% Read Data File (Main) Header
683%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
684function [hdr,msg_out]=l_ReadDataFileHeader(file_fid)
685
686msg_out='';
687%% Read Data File Header
688try
689  FH.nblocks = fread(file_fid,1,'long');   % Number of blocks in file
690  FH.ntraces = fread(file_fid,1,'long');   % Number of traces per block
691  FH.np = fread(file_fid,1,'long');        % Number of elements per trace
692  FH.ebytes = fread(file_fid,1,'long');    % Number of bytes per element
693  FH.tbytes = fread(file_fid,1,'long');    % Number of bytes per trace
694  FH.bbytes = fread(file_fid,1,'long');    % Number of bytes per block
695  FH.vers_id = fread(file_fid,1,'short');  % Software version, file_id status bits
696  FH.status = fread(file_fid,1,'short');   % Status of whole file
697  FH.nbheaders = fread(file_fid,1,'long'); % Number of block headers per block
698 
699  hdr.FileHeader = FH;
700 
701  %% Parse status bits
702  hdr.FileHeader.status=[];
703  tmp_str = {'S_DATA',...          % 0 = no data, 1 = data
704             'S_SPEC',...          % 0 = FID, 1 = spectrum
705             'S_32',...            % 0 = 16 bit, 1 = 32 bit integer
706             'S_FLOAT',...         % 0 = integer, 1 = floating point
707             'S_COMPLEX',...       % 0 = real, 1 = complex
708             'S_HYPERCOMPLEX',...  % 1 = hypercomplex
709             '',...                % Unused bit
710             'S_ACQPAR',...        % 0 = not Acqpar, 1 = Acqpar
711             'S_SECND',...         % 0 = first FT, 1 = second FT
712             'S_TRANSF',...        % 0 = regular, 1 = transposed
713             '',...                % Unused bit
714             'S_NP',...            % 1 = np dimension is active
715             'S_NF',...            % 1 = nf dimension is active
716             'S_NI',...            % 1 = ni dimension is active
717             'S_NI2',...           % 1 = ni2 dimension is active
718             ''...                 % Unused bit
719            };
720  status_bits = fliplr(double(dec2bin(FH.status,16))==49);
721  for ii=1:length(tmp_str)
722    if ~isempty(tmp_str{ii})
723      hdr.FileHeader.status.(tmp_str{ii}) = status_bits(ii);
724    end
725  end
726catch
727  hdr=[];
728  msg_out=['Error occurred while reading file header from "' filename '".'];
729  return
730end
731
732
733%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
734% Read Block Headers and Data
735%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
736function [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar)
737
738hdr.BlockHeaders = [];
739%hdr.HypercmplxBHeaders = [];
740data=[];
741kspace=[];
742count = [];
743msg_out='';
744
745if isfield(procpar,'seqcon')
[103]746  procpar.seqcon=procpar.seqcon{1};
[2]747else
[103]748  procpar.seqcon='nnnnn';
[2]749end
750
751%% Determine acquisition type from seqcon parameter
[103]752if all(procpar.seqcon=='n')          % 1D spectroscopy
753  Dat.AcqType=1;
754elseif all(procpar.seqcon(4:5)=='n') % 2D imaging
755  Dat.AcqType=2;
756elseif procpar.seqcon(5)=='n'        % 3D imaging
757  Dat.AcqType=3;
[2]758else                         % 4D imaging
[103]759  Dat.AcqType=4;
[2]760end
761
762%% Double-check that AcqType is not 1D spectral
763if isfield(procpar,'nv') && procpar.nv==0
[103]764  Dat.AcqType=1;
[2]765end
766
767%% Use some heuristical approach to "triple-check" that the data is not
768%  1D spectral
769if ~isempty(strfind(procpar.seqfil{1},'STEAM')) || ...
770      ~isempty(strfind(procpar.seqfil{1},'steam')) || ...
771      ~isempty(strfind(procpar.seqfil{1},'LASER')) || ...
772      ~isempty(strfind(procpar.seqfil{1},'laser')) || ...
773      ~isempty(strfind(procpar.seqfil{1},'PRESS')) || ...
774      ~isempty(strfind(procpar.seqfil{1},'press')) || ...
775      ~isempty(strfind(procpar.seqfil{1},'1PULSE')) || ...
776      ~isempty(strfind(procpar.seqfil{1},'1pulse')) || ...
777      ~isempty(strfind(procpar.seqfil{1},'2PULSE')) || ...
778      ~isempty(strfind(procpar.seqfil{1},'2pulse'))
[103]779  Dat.AcqType=1;
[2]780end
781
782UsePhaseTable=Dat.UsePhaseTable;
783
784
785%% If the data has been converted with flashc, the seqcon parameter needs
786% to be changed
787if isfield(procpar,'flash_converted')
788  if isfield(procpar,'ni') && procpar.ni>1
[103]789    procpar.seqcon(3)='s';
790  elseif ((procpar.seqcon(2)=='c') && (procpar.seqcon(3)=='c'))
791    procpar.seqcon(2)='s';
[2]792  end
793  UsePhaseTable=false; % Do not try to order data if flashc has been run
794end
795
796%% Set number of transients
797if ~isfield(procpar,'nt')
798  nt=1;
799else
800  nt=procpar.nt;
801end
802
803%% Determine if the arrayed acquisition was used
[53]804if isfield(procpar,'array') && ~isempty(procpar.array{1})
[2]805 
[60]806  if length(procpar.array)==1 && ~iscell(procpar.array{1}) && ...
807      strcmp(procpar.array{1},'pad') && all(procpar.pad==0)
808    % Skip the array parsing if the array is a dummy "pad"-array...
[103]809    Dat.isAcqArrayed = false;
810    Dat.ArrayLength = 1;
[58]811  else
[103]812    Dat.isAcqArrayed = true;
813    Dat.ArrayLength = [];
[60]814   
815    % Determine array length
816    for ii=1:length(procpar.array)
817      if iscell(procpar.array{ii})
[103]818        Dat.ArrayLength(ii) = length(procpar.(procpar.array{ii}{1}));
[60]819      else
[103]820        Dat.ArrayLength(ii) = length(procpar.(procpar.array{ii}));
[60]821      end
822    end
[103]823    Dat.ArrayLength = prod(Dat.ArrayLength);
[2]824  end
825else
[103]826  Dat.isAcqArrayed = false;
827  Dat.ArrayLength = 1;
[2]828end
829
830%% Determine if the data is EPI data
[71]831if ( isfield(procpar,'readres') && isfield(procpar,'phaseres') ) || ...
832    ( isfield(procpar,'apptype') && strcmp(procpar.apptype{1},'imEPI') )
[2]833  Dat.isEPIdata = true;
834else
835  Dat.isEPIdata = false;
836end
837
838BlockHeadStatusLabels = {'S_DATA',...       % 0 = no data, 1 = data
839                    'S_SPEC',...          % 0 = FID, 1 = spectrum
840                    'S_32',...            % 0 = 16 bit, 1 = 32 bit integer
841                    'S_FLOAT',...         % 0 = integer, 1 = floating point
842                    'S_COMPLEX',...       % 0 = real, 1 = complex
843                    'S_HYPERCOMPLEX',...  % 1 = hypercomplex
844                    '',...                % Unused bit
845                    'MORE_BLOCKS',...     % 0 = absent, 1 = present
846                    'NP_CMPLX',...        % 0 = real, 1 = complex
847                    'NF_CMPLX',...        % 0 = real, 1 = complex
848                    'NI_CMPLX',...        % 0 = real, 1 = complex
849                    'NI2_CMPLX',...       % 0 = real, 1 = complex
850                    '',...                % Unused bit
851                    '',...                % Unused bit
852                    '',...                % Unuesd bit
853                    ''...                 % Unused bit
854                   };
855
856BlockHeadModeLabels = {'NP_PHMODE',...   % 1 = ph mode
857                    'NP_AVMODE',...    % 1 = av mode
858                    'NP_PWRMODE',...   % 1 = pwr mode
859                    '',...             % Unused bit
860                    'NF_PHMODE',...    % 1 = ph mode
861                    'NF_AVMODE',...    % 1 = av mode
862                    'NF_PWRMODE',...   % 1 = pwr mode
863                    '',...             % Unused bit
864                    'NI_PHMODE',...    % 1 = ph mode
865                    'NI_AVMODE',...    % 1 = av mode
866                    'NI_PWRMODE',...   % 1 = pwr mode
867                    '',...             % Unused bit
868                    'NI2_PHMODE',...   % 1 = ph mode
869                    'NI2_AVMODE',...   % 1 = av mode
870                    'NI2_PWRMODE',...  % 1 = pwr mode
871                    ''...              % Unused bit
872                   };
873
874
875% The nbheaders -field is not correct in some cases. Thus, this field
876% cannot be trusted and the real nbheaders has to be calculated from the
877% byte counts.
878tmp_bytes=hdr.FileHeader.bbytes-hdr.FileHeader.tbytes*hdr.FileHeader.ntraces;
879header_bytes=28;
880if rem(tmp_bytes,header_bytes)~=0
881  msg_out = 'Block header byte count does not match with file header';
882  return
883else
884  nbheaders = tmp_bytes/28;
885end
[37]886%nbheaders = hdr.FileHeader.nbheaders;
[2]887
888%% Allocate space for k-space
[84]889% kspace=zeros(hdr.FileHeader.np/2,...
890%              hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);
[2]891
[84]892if ~Dat.FastDataRead
[103]893  if any(Dat.AcqType==[1 2])
894    switch procpar.seqcon(2:3)
[84]895      case {'cc','sc'}
896        kspace = zeros(hdr.FileHeader.np/2,...
897          hdr.FileHeader.ntraces,...
898          hdr.FileHeader.nblocks,Dat.precision);
899      otherwise
900        kspace = zeros(hdr.FileHeader.np/2,...
901          hdr.FileHeader.nblocks,...
902          hdr.FileHeader.ntraces,Dat.precision);
903    end
904  else
[2]905    kspace = zeros(hdr.FileHeader.np/2,...
[84]906      hdr.FileHeader.ntraces,...
907      hdr.FileHeader.nblocks,Dat.precision);
[2]908  end
909else
[84]910  kspace = [];
[2]911end
912
[21]913%% - The older robust (but also slower) way for reading the data.
914%% When the blocksize is relatively small, this is also quite fast.
915if ~Dat.FastDataRead
[2]916
[21]917  % Initialize waitbar
918  if Dat.ShowWaitbar
[80]919        wb_h = aedes_wbar(0/hdr.FileHeader.nblocks,...
[103]920          {['Reading ',num2str(Dat.AcqType),'D VNMR data (seqcon: "' procpar.seqcon '")'],...
[21]921          ['(Processing data block ' ...
922          num2str(0) '/' num2str(hdr.FileHeader.nblocks) ')']});
923  end
[2]924
[21]925  %% Read data and block headers
926  for ii=1:hdr.FileHeader.nblocks
927        %% Update waitbar
928        if Dat.ShowWaitbar
[80]929          aedes_wbar(ii/hdr.FileHeader.nblocks,...
[21]930                wb_h,...
[103]931                {['Reading ',num2str(Dat.AcqType),'D VNMR data (seqcon: "' procpar.seqcon '")'],...
[21]932                ['(Processing data block ' ...
933                num2str(ii) '/' num2str(hdr.FileHeader.nblocks) ')']})
934        end
[2]935
[21]936        %% Read block header and hypercomplex header
937        for kk=1:nbheaders
938          %% Read block header
939          if kk==1
940                hdr.BlockHeaders.scale = fread(file_fid,1,'short');  % Scaling factor
941                tmp_status = fread(file_fid,1,'short'); % Status of data in block
942                hdr.BlockHeaders.status = [];
943                hdr.BlockHeaders.index = fread(file_fid,1,'short');  % Block index
944                tmp_mode = fread(file_fid,1,'short');   % Mode of data in block
945                hdr.BlockHeaders.mode = [];
946                hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
947                hdr.BlockHeaders.lpval = fread(file_fid,1,'float');  % f2 (2D-f1) left phase in phasefile
948                hdr.BlockHeaders.rpval = fread(file_fid,1,'float');  % f2 (2D-f1) right phase in phasefile
949                hdr.BlockHeaders.lvl = fread(file_fid,1,'float');    % level drift correction
950                hdr.BlockHeaders.tlt = fread(file_fid,1,'float');    % tilt drift correction
951
952                %% Parse status and mode bits
953                status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
954                mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
955                for tt=1:length(BlockHeadStatusLabels)
956                  if ~isempty(BlockHeadStatusLabels{tt})
957                        hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
958                  end
959                  if ~isempty(BlockHeadModeLabels{tt})
960                        hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
961                  end
962                end
963
964
965          else %% Read hypercomplex header
966                fread(file_fid,1,'short'); % Spare
967                hdr.BlockHeaders.HCHstatus = fread(file_fid,1,'short');
968                fread(file_fid,1,'short'); % Spare
969                fread(file_fid,1,'short'); % Spare
970                fread(file_fid,1,'long'); % Spare
971                hdr.BlockHeaders.HCHlpval1 = fread(file_fid,1,'float');
972                hdr.BlockHeaders.HCHrpval1 = fread(file_fid,1,'float');
973                fread(file_fid,1,'float'); % Spare
974                fread(file_fid,1,'float'); % Spare
975          end
976        end
[37]977       
978        %% Check block index to be sure about the data type
979        if hdr.BlockHeaders.index~=ii
980          fprintf(1,'Warning: Index mismatch in "%s" block %d\n',fopen(file_fid),ii);
981         
982          % Use information from the file header instead of the BlockHeader if
983          % there is a mismatch in blockheader index...
984          useFileHeader = true;
985        else
986          useFileHeader = false;
987        end
988       
[21]989        %% Determine data precision
[37]990        if useFileHeader
991          if hdr.FileHeader.status.S_FLOAT==1
992                prec_str = ['single=>',Dat.precision];
993          elseif hdr.FileHeader.status.S_32==1 ...
994                  && hdr.FileHeader.status.S_FLOAT==0
995                prec_str = ['int32=>',Dat.precision];
996          elseif hdr.FileHeader.status.S_32==0 ...
997                  && hdr.FileHeader.status.S_FLOAT==0
998                prec_str = ['int16=>',Dat.precision];
999          end
1000         
1001        else
1002          if hdr.BlockHeaders.status.S_FLOAT==1
1003                prec_str = ['single=>',Dat.precision];
1004          elseif hdr.BlockHeaders.status.S_32==1 ...
1005                  && hdr.BlockHeaders.status.S_FLOAT==0
1006                prec_str = ['int32=>',Dat.precision];
1007          elseif hdr.BlockHeaders.status.S_32==0 ...
1008                  && hdr.BlockHeaders.status.S_FLOAT==0
1009                prec_str = ['int16=>',Dat.precision];
1010          end
[21]1011        end
1012
1013        % Read k-space
1014        tmp=fread(file_fid,...
1015          [hdr.FileHeader.np,hdr.FileHeader.ntraces],...
1016          prec_str);
1017
1018
1019        % Store complex block and perform DC correction
1020        if ~Dat.DCcorrection || ( nt(1)>1 )
1021          complex_block = complex(tmp(1:2:end,:),tmp(2:2:end,:));
1022        else
1023          complex_block = complex(tmp(1:2:end,:)-hdr.BlockHeaders.lvl,...
1024                tmp(2:2:end,:)-hdr.BlockHeaders.tlt);
1025        end
1026
1027
1028        %% Store and order k-space values
[103]1029        if any(Dat.AcqType==[1 2])
1030          switch procpar.seqcon(2:3)
[21]1031                case {'cc','sc'}
1032                  kspace(:,:,ii) = complex_block;
1033                otherwise
1034                  kspace(:,ii,:) = complex_block;
1035          end
1036        else
1037          kspace(:,:,ii) = complex_block;
1038        end
1039
1040        % Do not save blockheaders by default. When reading data files with a lot of
1041        % blocks (e.g. over 1000) the processing of the DATA structure can be
1042        % slowed down considerably. If you for some reason want to save also the
1043        % block headers in the DATA structure comment out the code line below.
[84]1044        hdr.BlockHeaders = [];
[21]1045  end % for ii=1:hdr.
1046
[2]1047  if Dat.ShowWaitbar
[21]1048        delete(wb_h)
1049  end
[2]1050 
[21]1051else
[22]1052  %% -------------------------------------------------------------------
[21]1053  %% Fast Method for reading data. This may fail with some datas and can
1054  %% also require a relatively large amount of memory. This
1055  %% method should be used for EPI datas that contain a large number
1056  %% of block headers...
1057 
1058  % Initialize waitbar
1059  if Dat.ShowWaitbar
[103]1060        wb_h = aedes_calc_wait(['Reading ',num2str(Dat.AcqType),...
1061          'D VNMR data (seqcon: "' procpar.seqcon '")']);
[2]1062  end
1063 
[22]1064  % The first block header is read and it is assumed that the values in
1065  % the other block headers don't change. 
1066  hdr.BlockHeaders.scale = fread(file_fid,1,'short');  % Scaling factor
1067  tmp_status = fread(file_fid,1,'short'); % Status of data in block
1068  hdr.BlockHeaders.status = [];
1069  hdr.BlockHeaders.index = fread(file_fid,1,'short');  % Block index
1070  tmp_mode = fread(file_fid,1,'short');   % Mode of data in block
1071  hdr.BlockHeaders.mode = [];
1072  hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
1073  hdr.BlockHeaders.lpval = fread(file_fid,1,'float');  % f2 (2D-f1) left phase in phasefile
1074  hdr.BlockHeaders.rpval = fread(file_fid,1,'float');  % f2 (2D-f1) right phase in phasefile
1075  hdr.BlockHeaders.lvl = fread(file_fid,1,'float');    % level drift correction
1076  hdr.BlockHeaders.tlt = fread(file_fid,1,'float');    % tilt drift correction
1077 
1078  %% Parse status and mode bits
1079  status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
1080  mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
1081  for tt=1:length(BlockHeadStatusLabels)
[78]1082    if ~isempty(BlockHeadStatusLabels{tt})
1083      hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
1084    end
1085    if ~isempty(BlockHeadModeLabels{tt})
1086      hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
1087    end
[22]1088  end
1089 
1090  %% Determine data precision
1091  if hdr.BlockHeaders.status.S_FLOAT==1
[78]1092    prec_str = ['single=>',Dat.precision];
1093    prec = 4; % Precision in bytes
[22]1094  elseif hdr.BlockHeaders.status.S_32==1 ...
[78]1095      && hdr.BlockHeaders.status.S_FLOAT==0
1096    prec_str = ['int32=>',Dat.precision];
1097    prec = 4;
[22]1098  elseif hdr.BlockHeaders.status.S_32==0 ...
[78]1099      && hdr.BlockHeaders.status.S_FLOAT==0
1100    prec_str = ['int16=>',Dat.precision];
1101    prec = 2;
[2]1102  end
1103 
[22]1104  % Seek the file back to the beginning of the first block header
1105  fseek(file_fid,-28,0);
1106 
[21]1107  % Determine the number of values that will result from block header(s)
1108  nVals = (nbheaders*28)/prec;
[2]1109 
[21]1110  % Read the whole data including block headers etc...
[78]1111  kspace = fread(file_fid,inf,prec_str);
1112  kspace=reshape(kspace,nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces,[]);
[2]1113 
[21]1114  % Remove block headers from the data
[78]1115  kspace=kspace(nVals+1:end,:);
[2]1116 
[21]1117  % Transform to complex values
[78]1118  kspace = reshape(kspace,hdr.FileHeader.np,...
[21]1119        hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);
[2]1120 
[22]1121  % Do DC-correction if necessary
1122  if ~Dat.DCcorrection || ( nt(1)>1 )
[78]1123    kspace=complex(kspace(1:2:end,:,:),kspace(2:2:end,:,:));
[22]1124  else
[78]1125    kspace=complex(kspace(1:2:end,:,:)-hdr.BlockHeaders.lvl,...
1126      kspace(2:2:end,:,:)-hdr.BlockHeaders.tlt);
[22]1127  end
[2]1128 
1129  %% Store and order k-space values
[103]1130  if any(Dat.AcqType==[1 2])
1131    switch procpar.seqcon(2:3)
[78]1132      case {'cc','sc'}
1133        %kspace(:,:,ii) = complex_block;
1134      otherwise
1135        %kspace(:,ii,:) = complex_block;
1136        kspace = permute(kspace,[1 3 2]);
1137    end
[2]1138  else
[78]1139    %kspace(:,:,ii) = complex_block;
[2]1140  end
1141 
[21]1142  if Dat.ShowWaitbar
1143        delete(wb_h)
1144  end
1145 
1146end
[2]1147
1148% Remove singleton dimensions from kspace
1149kspace = squeeze(kspace);
1150
1151% Check if raw kspace should be returned
1152if Dat.ReturnRawKspace
1153  %% Delete waitbar
[36]1154  if Dat.ShowWaitbar && ishandle(wb_h)
1155        delete(wb_h)
1156  end
[2]1157  return
1158end
1159
[35]1160%% Support for RASER sequence ---------------------------
1161if isfield(procpar,'teType')
1162 
1163  if strcmpi(procpar.sing_sh,'y') && strcmpi(procpar.teType,'c')
1164       
[71]1165    % Separate reference scans from the data
1166    if isfield(procpar,'refscan') && strcmp(procpar.refscan,'y')
1167      kspace=permute(reshape(kspace,procpar.np/2,procpar.ne/2,2,[]),[1 2 4 3]);
1168      noise = kspace(:,:,:,1);
1169      kspace = kspace(:,:,:,2);
1170      % DC-correction
1171      dc_level = squeeze(mean(mean(noise)));
1172      for ii=1:size(kspace,3)
1173        kspace(:,:,ii)=kspace(:,:,ii)-dc_level(ii);
1174      end
1175    end
1176   
1177    if strcmpi(procpar.readType,'s')
1178      sw = procpar.sw;
1179      dwell = cumsum(ones(1,size(kspace,1))*(1/sw));
1180    else
1181      error('This feature has not been implemented yet!!!')
1182    end
1183    dwell = dwell*1000;
1184    %kspace = permute(kspace,[2 1 3]);
1185   
1186   
1187   
1188    % Phase correction
1189    dims = size(kspace);
1190    nI = size(kspace,3);
1191    nv = size(kspace,4);
1192    ne = size(kspace,2);
1193    np = size(kspace,1);
1194    g = 42.58;
1195    G = procpar.gvox2;
1196    Ds = procpar.vox2/100;
1197    pos2 = procpar.pos2/10;
1198    fudge = 0;
1199    fudge2 = 1;
[83]1200   
[71]1201    tt_dwell=repmat(dwell(:),1,ne);
1202    tt_ne = repmat((1:ne),length(dwell),1);
1203    tt_np = repmat((1:np)',1,ne);
1204   
1205    i=sqrt(-1);
1206    phi = (-2.*pi.*g.*G.*Ds.*tt_dwell.*(tt_ne./ne - 1/2 - pos2/Ds))+fudge.*tt_np;
1207    phi = repmat(phi,[1,1,size(kspace,3)]);
1208    kspace = abs(kspace).*exp(i*(angle(kspace)+phi));
1209   
1210    % Flip every other column
1211    sz(1)=size(kspace,1);
1212    sz(2)=size(kspace,2);
1213    sz(3)=size(kspace,3);
1214    sz(4)=size(kspace,4);
1215    kspace=reshape(kspace,sz(1),prod(sz(2:4)));
1216    kspace(:,1:2:end) = flipud(kspace(:,1:2:end));
1217    kspace = reshape(kspace,sz);
1218   
1219    %kspace = reshape(permute(kspace,[2 1 3]),sz(2),[],1);
1220    %kspace(:,1:2:end) = flipud(kspace(:,1:2:end));
1221    %kspace=permute(reshape(kspace,sz(2),[],sz(3)),[2 1 3]);
1222   
[35]1223  else
[71]1224    data = [];
1225    error('RASER sequence of this type has not been implemeted yet!')
[35]1226  end
1227 
1228  % Reshape into 4D matrix
1229  kspace = reshape(kspace,size(kspace,1),size(kspace,2),[],size(kspace,3));
1230 
[67]1231  % Reorient if requested
[75]1232  if Dat.OrientImages
[67]1233    kspace = flipdim(kspace,2);
1234  end
1235 
[36]1236  data = abs(fftshift(fft(fftshift(fft(kspace,[],3),3),[],1),1));
[35]1237 
1238  % Delete kspace if not returned
1239  if ~Dat.ReturnKSpace
[71]1240    kspace=[];
[35]1241  end
1242 
1243  return
1244 
1245end
1246
[103]1247%% Support for fat/water measurements ----------------------
1248if isfield(procpar,'seqfil') && strcmpi(procpar.seqfil,'ge3d_csi2')
1249 
1250  % Split kspace into fat and water (in 4th dimesion)
1251  kspace=reshape(permute(reshape(kspace,256,2,[]),[1 3 4 2]),256,128,[],2);
1252 
1253  % Fourier transform data
1254  if any(Dat.ZeroPadding==[1 2])
1255    data_sz = [procpar.np/2,procpar.nf,procpar.nv2,2];
1256    data = zeros(data_sz,class(kspace));
1257    data(:,:,:,1) = abs(fftshift(fftn(kspace(:,:,:,1),data_sz(1:3))));
1258    data(:,:,:,2) = abs(fftshift(fftn(kspace(:,:,:,2),data_sz(1:3))));
1259  else
1260    data = zeros(size(kspace),class(kspace));
1261    data(:,:,:,1) = abs(fftshift(fftn(kspace(:,:,:,1))));
1262    data(:,:,:,2) = abs(fftshift(fftn(kspace(:,:,:,2))));
1263  end
1264 
1265  % Delete kspace if not returned
1266  if ~Dat.ReturnKSpace
1267    kspace=[];
1268  end
1269 
1270  return
1271end
[35]1272
[103]1273% Check the number of receivers (PI stuff)
1274if isfield(procpar,'rcvrs') && ~isempty(procpar.rcvrs)
1275  % Multiple receivers used
1276  nRcv = length(find(procpar.rcvrs{1}=='y'));
1277  data = [];
1278  kspace2 = [];
1279  for ii=1:nRcv
1280    tmp_kspace = l_ReconstructKspace(kspace(:,:,ii),procpar,Dat);
1281    kspace2(:,:,:,ii)=tmp_kspace;
1282    if Dat.ReturnFTData
1283      tmp_data = l_CalculateFFT(tmp_kspace,procpar,Dat);
1284      data(:,:,:,ii) = tmp_data;
1285    end
1286  end
1287  kspace = kspace2;
1288  kspace2=[];
1289else
1290  % Only one receiver
1291  kspace = l_ReconstructKspace(kspace,procpar,Dat);
1292  data = l_CalculateFFT(kspace,procpar,Dat);
1293end
1294
1295% Delete kspace if not returned
1296if ~Dat.ReturnKSpace
1297  kspace=[];
1298end
1299
1300%===============================================
1301% Reconstruct kspace for different sequences
1302%===============================================
1303function kspace=l_ReconstructKspace(kspace,procpar,Dat)
1304
1305
[2]1306%% Flip images for certain measurements
1307if Dat.FlipKspace~=0
1308  if ischar(Dat.FlipInd)
1309    if strcmpi(Dat.FlipInd,'all')
1310      flipind = 1:size(kspace,3);
1311    elseif strcmpi(Dat.FlipInd,'alt')
1312      flipind = 2:2:size(kspace,3);
1313    end
1314  else
1315    flipind = Dat.FlipInd;
1316  end
1317 
1318  for ii=flipind
1319    if Dat.FlipKspace==1
1320      kspace(:,:,ii)=fliplr(kspace(:,:,ii));
1321    elseif Dat.FlipKspace==2
1322      kspace(:,:,ii)=flipud(kspace(:,:,ii));
1323    else
1324      kspace(:,:,ii)=flipud(fliplr(kspace(:,:,ii)));
1325    end
1326  end
1327end
1328
1329% Handle arrayed and RARE sequences
[103]1330if Dat.isAcqArrayed && Dat.AcqType~=1 && Dat.Sorting && ~Dat.isEPIdata
1331  %if ~strcmpi(procpar.seqcon(2:3),'cc') && ( isempty(Dat.phasetable) | ...
[2]1332  %        isfield(procpar,'flash_converted') )
[103]1333  if (( isempty(Dat.phasetable) && ~strcmpi(procpar.seqcon(2:3),'sc') ) && ...
1334      ~(strcmpi(procpar.seqcon(2:3),'cc') && Dat.ArrayLength==size(kspace,3))) || ...
[81]1335      isfield(procpar,'flash_converted')
[2]1336                                 
1337    % Sort uncompressed arrayed data
[103]1338    ks_order = 1:procpar.nv*Dat.ArrayLength;
[2]1339    tmp = 0:procpar.nv:(procpar.ne-1)*procpar.nv;
[103]1340    tmp=repmat(tmp,procpar.nv*Dat.ArrayLength,1);
1341    ks_order = reshape(ks_order,Dat.ArrayLength,procpar.nv).';
[2]1342    ks_order = ks_order(:);
1343    ks_order = repmat(ks_order,1,procpar.ne);
1344    ks_order = ks_order+tmp;
1345    ks_order = ks_order.';
1346   
1347    if procpar.ns==1
1348      kspace=kspace(:,ks_order);
1349    else
1350      kspace = kspace(:,ks_order,:);
1351    end
1352   
1353    % Reshape into 3D matrix
[103]1354    kspace=reshape(kspace,[size(kspace,1) procpar.nv Dat.ArrayLength* ...
[2]1355                        procpar.ns]);
1356  else
1357    %% Sort RARE type data
[81]1358    if UsePhaseTable && ~isempty(Dat.phasetable)
[93]1359      %Dat.phasetable = Dat.phasetable';
[103]1360      if Dat.isAcqArrayed && Dat.ArrayLength>1 && strcmpi(procpar.seqcon(2:3),'cs')
1361        kspace = permute(reshape(reshape(kspace,size(kspace,1),[]),size(kspace,1),Dat.ArrayLength,[]),[1 3 2]);
[95]1362      else
1363        Dat.phasetable = Dat.phasetable.';
[93]1364      end
[2]1365      kspace(:,Dat.phasetable(:),:) = kspace;
1366    end
1367  end
1368elseif Dat.isEPIdata
[71]1369  %% Support for EPI-measurements
[2]1370 
[28]1371 
[71]1372  % EPI data is measured with old VNMR
1373  if isfield(procpar,'readres') && isfield(procpar,'phaseres')
1374   
1375    % Number of slices
1376    tmp_ns=length(procpar.pss);
1377   
1378    if isfield(procpar,'navecho') && strcmpi(procpar.navecho{1},'y')
1379      tmp_nv = procpar.nv-procpar.nseg;
1380    else
1381      tmp_nv = procpar.nv;
1382    end
1383    kspace = reshape(kspace,[size(kspace,1) ...
1384      size(kspace,2)/tmp_nv tmp_nv]);
1385    kspace = permute(kspace,[1 3 2]);
1386   
1387    % Reshape to 4D matrix
1388    kspace = reshape(kspace,[size(kspace,1) size(kspace,2) ...
1389      tmp_ns size(kspace,3)/tmp_ns]);
1390   
1391   
1392  elseif isfield(procpar,'apptype') && strcmp(procpar.apptype{1},'imEPI')
1393    % If EPI data is measured with new VNMRj system
1394   
1395    % Number of slices
1396    ns = length(procpar.pss);
1397   
1398    % Reshape kspace to 4D matrix
[99]1399    if isfield(procpar,'fract_ky') && procpar.fract_ky~=procpar.nv/2
1400      kspace = reshape(kspace,procpar.np/2,procpar.nv/2+procpar.fract_ky,...
1401        ns,[]);
1402      kspace(procpar.np/2,procpar.nv,1)=0;
1403    else
1404      kspace = reshape(kspace,procpar.np/2,procpar.nv,...
1405        ns,[]);
1406    end
[71]1407   
1408    % Flip even kspace lines
[84]1409    if ~Dat.FastDataRead
1410      for tt=2:2:size(kspace,2)
1411        kspace(:,tt,:,:) = flipdim(kspace(:,tt,:,:),1);
1412      end
1413    else
1414      kspace(:,2:2:end,:,:) = flipdim(kspace(:,2:2:end,:,:),1);
1415    end
[71]1416   
[84]1417   
[71]1418    % Sort centric measurements
1419    if isfield(procpar,'ky_order') && strcmpi(procpar.ky_order,'c')
1420      kspace(:,Dat.phasetable(:),:,:)=kspace;
1421    end
1422    %kspace = flipdim(kspace,1);
1423   
1424    % Get the reference image(s)
1425    ref_im = kspace(:,:,:,1);
1426   
1427    % Do phase correction.
1428    phase_e = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
1429    for kk=2:size(kspace,4)
1430      kspace(:,:,:,kk) = ifft(fft(kspace(:,:,:,kk),[],1).*phase_e,[],1);
1431    end
[38]1432  end
[2]1433 
1434else
1435  % This section contains various "chewing gum and ironwire" type
1436  % patches that hopefully work...
1437 
[103]1438  if strcmp(procpar.seqcon,'nncnn')
[82]1439    kspace = permute(kspace,[1 3 2]);
[103]1440  elseif strcmp(procpar.seqcon,'nccnn') && length(size(kspace))==2 && ...
1441      procpar.ns>=1 && Dat.AcqType~=1
[2]1442    if ~isempty(Dat.phasetable)
[82]1443      kssz = size(kspace);
1444      phsz = size(Dat.phasetable);
1445      kspace=permute(reshape(kspace,procpar.np/2,...
1446        phsz(2),...
1447        kssz(2)/phsz(2)),[1 3 2]);
1448      kspace = permute(reshape(kspace,...
1449        procpar.np/2,procpar.ns,kssz(2)/procpar.ns),[1 3 2]);
1450      if Dat.Sorting
1451        kspace(:,Dat.phasetable(:),:)=kspace;
1452      end
1453    else
1454      kspace = reshape(kspace,...
1455        [procpar.np/2 procpar.ns ...
1456        size(kspace,2)/procpar.ns]);
1457      kspace=permute(kspace,[1 3 2]);
1458    end
1459   
[103]1460  elseif strcmp(procpar.seqcon,'nscsn') && length(size(kspace))==3 && ...
1461          Dat.AcqType~=1
[23]1462        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1463        % Support for 3D fast spin-echo
1464        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1465        if ~isempty(Dat.phasetable)
1466          for ii=1:size(kspace,3)
1467                tmp = kspace(:,:,ii);
1468                kssz = size(tmp);
1469                phsz = size(Dat.phasetable);
1470                tmp=permute(reshape(tmp,procpar.np/2,...
1471                  phsz(2),...
1472                  kssz(2)/phsz(2)),[1 3 2]);
1473                tmp = permute(reshape(tmp,...
1474                  procpar.np/2,procpar.ns,kssz(2)/procpar.ns),[1 3 2]);
1475                if Dat.Sorting
1476                  tmp(:,Dat.phasetable(:),:)=tmp;
1477                end
1478                kspace(:,:,ii)=tmp;
1479          end
1480        end
1481       
1482       
[103]1483  elseif Dat.AcqType~=1 && isfield(procpar,'nv')
[2]1484    if length(size(kspace))==2 && ...
[21]1485              ( size(kspace,1)~=(procpar.np/2) || size(kspace,2)~=procpar.nv )
[2]1486      %% Reshape the kspace to the appropriate size
1487      kspace=reshape(kspace,[procpar.np/2 procpar.nv size(kspace,2)/ ...
1488                          procpar.nv]);
[89]1489    elseif length(size(kspace))==3
1490      kspace = reshape(kspace,procpar.np/2,procpar.nv,[]);
[94]1491      if Dat.Sorting && ~isempty(Dat.phasetable)
[89]1492        kspace(:,Dat.phasetable(:),:) = kspace;
1493      end
[2]1494    end
1495  end
1496 
1497 
1498  %%% Support for Teemu's ASE3D-data
[103]1499  %if strcmpi(procpar.seqcon,'ncccn')
[2]1500  %  %% Reshape the kspace to the appropriate size
1501  %  kspace=reshape(kspace,[procpar.np/2 procpar.nv size(kspace,2)/procpar.nv]);
1502  %end
1503end
1504
[71]1505
[103]1506
1507%=========================================
1508% Fourier Transform Data
1509%=========================================
1510function data=l_CalculateFFT(kspace,procpar,Dat)
1511data=[];
1512
[2]1513% Return image/spectral data --------------------------------
1514if Dat.ReturnFTData
1515  %% Fourier transform spectral data
[103]1516  if Dat.AcqType==1
[80]1517        wb_h = aedes_wbar(0,'Fourier transforming spectral data');
[22]1518    %data=zeros(size(kspace),'single');
1519        data=kspace;
[2]1520    sz=size(kspace,2)*size(kspace,3);
1521    count=1;
1522    for kk=1:size(kspace,3)
1523      for ii=1:size(kspace,2)
1524        %% Update waitbar
1525        if Dat.ShowWaitbar
[80]1526          aedes_wbar(count/sz,...
[2]1527               wb_h,...
1528               {'Fourier transforming spectral data',...
1529                ['(Processing spectrum ' num2str(count) '/' num2str(sz) ')']})
[22]1530                end
1531        data(:,ii,kk) = abs(fftshift(fft(data(:,ii,kk))));
[2]1532        count=count+1;
1533      end
1534    end
1535 
1536  %% Fourier transform image data
1537  else
1538   
1539    %% Zeropadding auto
1540    if Dat.ZeroPadding==2
1541     
[71]1542      if Dat.isEPIdata && isfield(procpar,'readres') && isfield(procpar,'phaseres')
[2]1543       
1544        %% Determine the image size for EPI data for procpar "readres"
1545        %% and "phaseres" fields
1546        data_sz = [procpar.readres procpar.phaseres];
[71]1547        if isequal(data_sz,[size(kspace,1),size(kspace,2)])
1548          DoZeroPadding=false;
1549        else
1550          DoZeroPadding=true;
1551        end
1552        data_sz = [procpar.readres ...
1553          procpar.phaseres size(kspace,3) size(kspace,4)];
[2]1554      else
1555       
[71]1556        %% Check if zeropadding is necessary.
[2]1557        ks_sz_sorted = sort([size(kspace,1),size(kspace,2)]);
1558        fov_sz_sorted = sort([procpar.lro,procpar.lpe]);
1559        sliceRelativeDim = ks_sz_sorted(2)/ks_sz_sorted(1);
1560        FOVrelativeDim = fov_sz_sorted(2)/fov_sz_sorted(1);
1561        if FOVrelativeDim~=sliceRelativeDim
1562          ind=find([size(kspace,1),size(kspace,2)]==ks_sz_sorted(1));
1563          data_sz = size(kspace);
1564          data_sz(ind) = round((fov_sz_sorted(1)/fov_sz_sorted(2))*ks_sz_sorted(2));
1565          DoZeroPadding=true;
1566        else
1567          data_sz = size(kspace);
1568          DoZeroPadding=false;
1569        end
1570      end
1571     
[71]1572      %% Force zeropadding on. Force images to be square.
[2]1573    elseif Dat.ZeroPadding==1
1574      ks_sz_max = max(size(kspace,1),size(kspace,2));
1575      data_sz = size(kspace);
1576      data_sz(1:2)=ks_sz_max;
1577      DoZeroPadding=true;
1578     
[71]1579      %% Force zeropadding off
[2]1580    else
1581      data_sz = size(kspace);
1582      DoZeroPadding=false;
1583    end
1584   
1585    %% Fourier transform 2D and EPI image data
[103]1586    if Dat.AcqType==2 || Dat.isEPIdata
[2]1587      sz=size(kspace,3);
[71]1588     
1589      if DoZeroPadding
1590        data=zeros(data_sz,class(kspace));
1591        data(1:size(kspace,1),1:size(kspace,2),:,:)=kspace;
1592      else
1593        data=kspace;
1594      end
[84]1595      if ~Dat.ReturnKSpace
1596        kspace = [];
1597      end
[71]1598      if Dat.ShowWaitbar
[80]1599        [wb_h,txh]=aedes_calc_wait('Fourier transforming image data...');
[71]1600      end
[84]1601     
[71]1602      %data=permute(fftshift(fftshift(abs(fft(permute(fft(data),[2 1 3 4]))),1),2),[2 1 3 4]);
[84]1603      if ~Dat.FastDataRead
1604        for tt=1:size(data,4)
1605          data(:,:,:,tt)=fftshift(fftshift(abs(fft(fft(data(:,:,:,tt),[],1),[],2)),1),2);
1606        end
1607      else
1608        data=fftshift(fftshift(abs(fft(fft(data,[],1),[],2)),1),2);
1609      end
[2]1610    else
1611      %% Fourier transform 3D image data
1612      if Dat.ShowWaitbar
[80]1613        [wb_h,txh]=aedes_calc_wait('Fourier transforming 3D image data...');
[2]1614      end
1615     
1616      if DoZeroPadding
[23]1617        % Check zeropadding in the 3rd dimension if zeropadding is set to
1618        % "auto"
[71]1619        if Dat.ZeroPadding==2
1620          if isfield(procpar,'lpe2') && isfield(procpar,'lpe')
[103]1621            data_sz(3) = max(1,round((procpar.lpe2/procpar.lpe)*data_sz(2)*Dat.ArrayLength));
[71]1622          end
1623        end
[38]1624        data=zeros(data_sz,class(kspace));
[2]1625        data = abs(fftshift(fftn(kspace,data_sz)));
1626      else
[38]1627        data=zeros(data_sz,class(kspace));
[71]1628        data = abs(fftshift(fftn(kspace)));
1629      end
[2]1630    end
1631   
[75]1632    % Remove reference image if requested
1633    if Dat.isEPIdata && Dat.RemoveEPIphaseIm
1634      data = data(:,:,:,2:end);
1635    end
1636   
[71]1637    % Reorient images and remove the reference image if requested
[75]1638    if Dat.OrientImages && ~isempty(procpar) && ...
[103]1639        isfield(procpar,'orient') && any(Dat.AcqType==[2 3 4])
[75]1640      orient = procpar.orient{1};
1641      if any(strcmpi(orient,{'xyz','trans90','cor90','sag90'}))
[80]1642        data = flipdim(aedes_rot3d(data,1,3),2);
[75]1643      elseif strcmpi(orient,'oblique')
1644        data = flipdim(flipdim(data,1),2);
[71]1645      else
[75]1646        data = flipdim(flipdim(data,1),2);
[71]1647      end
1648    end
1649   
[75]1650   
[71]1651   
[2]1652% $$$   %% Fourier transform 4D image data
[103]1653% $$$   elseif Dat.AcqType==4
[2]1654% $$$     data = abs(fftshift(fftshift(fft2(kspace),2),1));
1655  end
[103]1656 
1657  %% Delete waitbar
1658  if Dat.ShowWaitbar && ishandle(wb_h)
1659    delete(wb_h)
1660  end
[2]1661end
1662
1663
1664
[103]1665
1666
[81]1667%========================================================
1668% A function for creating phasetables for fse and fsems
1669%========================================================
1670function phasetable=l_CreateFsemsPhaseTable(procpar)
[2]1671
[81]1672phasetable = [];
1673if ~isfield(procpar,'etl') || ~isfield(procpar,'kzero')
1674  return
1675end
1676etl = procpar.etl;
1677kzero = procpar.kzero;
1678nv = procpar.nv;
[2]1679
[81]1680t = (-nv/2+1):(nv/2);
1681t = t(:);
1682tt = flipud(t);
1683tt = tt(:);
1684
[82]1685phasetable = [reshape(t(1:nv/2),[],etl);...
1686  flipud(reshape(tt(1:nv/2),[],etl))];
[81]1687phasetable = circshift(fliplr(phasetable),[0 kzero-1]);
1688
1689
[2]1690%% - EOF - %%
[82]1691return
Note: See TracBrowser for help on using the repository browser.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008