source: aedes_readfid.m @ 81

Last change on this file since 81 was 81, checked in by tjniskan, 10 years ago
  • Added (again) some "chewing gum and ironwire" into aedes_readfid.m.

The handling and sorting of VNMR data should be rewritten in the near
future.

  • Added a function in aedes_readfid.m for creating fse/fsems

phasetables. It is, however, not in use at the moment so phasetables
should still be in the tablib directory.

M aedes_readfid.m
M aedes_revision.m

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

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