source: an2_readfid.m @ 73

Last change on this file since 73 was 73, checked in by tjniskan, 11 years ago
  • Fixed a bug in aedes.m when viewing "voxel time-series" from generic

matrix input

  • Fixed a small issue with some arrayed experiments in an2_readfid.m

and hopefully did not break anything this time

  • Few minor corrections and changes here and there...

M misclib/fmri_filter.m
M misclib/fmri_corr.m
M an2_revision.m
M an2_fitmaps.m
M aedes.m
M plugins/resting_state_fc.m
M an2_readfid.m

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

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