source: an2_readfid.m @ 44

Last change on this file since 44 was 39, checked in by tjniskan, 11 years ago
  • Added "Recent Files" quick access menu in File-menu
  • Fixed a bug in plugins/copy_data_to_workspace.m that made writing a

new name for a variable difficult

  • Fixed a bug in the file prefix editbox in an2_export_gui.m
  • Changed references to Analyze 2.0 in all license notices to Aedes.

The name change should now be complete...

M an2_export_gui.m
M an2_cellsprintf.m
M an2_calc_wait.m
M an2_check_file_exist.m
M an2_iseven.m
M an2_cellwrite.m
M an2_wbar.m
M an2_rot3d.m
M an2_readfdf.m
M an2_revision.m
M an2_viewprocpar.m
M an2_checkcthdr.m
M an2_readprocpar.m
M an2_fitmaps.m
M an2_read_nifti.m
M an2_data_read.m
M an2_resviewer.m
M an2_maptool.m
M aedes.m
M an2_res2table.m
M an2_copy_roi.m
M plugins/save_roi_as_mask.m
M plugins/write_difference_images.m
M plugins/plot_profile.m
M plugins/calculate_t2_map.m
M plugins/calculate_t1r_map.m
M plugins/view_kspace.m
M plugins/copy_data_to_workspace.m
M plugins/take_snapshot.m
M an2_inputdlg.m
M an2_roi_copy_gui.m
M an2_readctdata.m
M an2_readfid.m
M an2_readfidprefs.m
M an2_readtab.m
M an2_check_updates.m
M an2_killfigs.m
M an2_roi_stats.m
M an2_saveres.m
M an2_rotateflip.m
M an2_juigetfiles.m
M an2_gui_defaults.m
M an2_editstack.m
M an2_update.m
M an2_write_nifti.m

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

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