source: aedes_readvnmr.m @ 153

Last change on this file since 153 was 153, checked in by tjniskan, 8 years ago
  • Added support for Various EPI sequences in vnmr_recon/epi_recon.m.

M aedes_readvnmr.m
M aedes_revision.m
M vnmr_recon/epi_recon.m

File size: 42.2 KB
Line 
1function [DATA,msg_out] = aedes_readvnmr(filename,varargin)
2% AEDES_READVNMR - Read VNMR (Varian) FID-files
3%   
4%
5% Synopsis:
6%        [DATA,msg]=aedes_readvnmr(filename,'PropertyName1',value1,'PropertyName2',value2,...)
7%        [DATA,msg]=aedes_readvnmr(filename,'header')
8%        [DATA,msg]=aedes_readvnmr(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_READVNMR 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_READVNMR (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 is not as
57%                                               % robust as the older
58%                                               % method and can consume a lot
59%                                               % of memory.
60%
61%        'SumOfSquares' :  [ {1} | 2 | 3 ]      % 1=Return only the
62%                                               % sum-of-squares image
63%                                               % for multireceiver
64%                                               % data (default).
65%                                               % 2=Return both SoQ and
66%                                               % individual receiver data
67%                                               % 3=Return only individual
68%                                               % receiver data
69%                                               % NOTE: This property has
70%                                               % no effect on single
71%                                               % receiver data.
72%
73%        'FileChunkSize' : [ megabytes ]        % Chunk size in megabytes for
74%                                               % reading fid-files
75%                                               % (default=500 MB).
76%                                               % Lowering the chunk size
77%                                               % helps to conserve memory
78%                                               % when reading large files,
79%                                               % but is slower. This
80%                                               % property has no effect if
81%                                               % FastRead='off'.
82%                                               
83%
84%        'OutputFile'  :  filename              % Writes the slices into
85%                                               % NIfTI-format files
86%                                               % using the given filename.
87%
88%        'DCcorrection':  [ {'off'} | 'on' ]    % 'on'=use DC correction
89%                                               % 'off'=don't use DC correction
90%                                               % (default)
91%
92%        'Procpar'     :  (procpar-structure)   % Input procpar
93%                                               % structure. If this
94%                                               % property is not set the
95%                                               % procpar structure
96%                                               % will be read using
97%                                               % AEDES_READPROCPAR
98%
99%        'ZeroPadding' :  ['off'|'on'|{'auto'}] % 'off' = force
100%                                               % zeropadding off
101%                                               % 'on' = force
102%                                               % zeropadding on (force
103%                                               % images to be square)
104%                                               % 'auto' = autoselect
105%                                               % using relative FOV
106%                                               % dimensions (default)
107%
108%        'PhaseTable'  : (custom phase table)   % User-specified
109%                                               % line order for
110%                                               % k-space. The phase table
111%                                               % is obtained from the file
112%                                               % specified in
113%                                               % procpar.petable if
114%                                               % necessary.
115%
116%        'sorting'      : [ 'off' | {'on'} ]    % 'off' =disable k-space
117%                                               % sorting, i.e. ignore the
118%                                               % phase table and/or arrays.
119%                                               % 'on' =sort k-space using
120%                                               % phase table and/or array
121%                                               % information if necessary
122%                                               % (default)
123%
124%        'wbar'        : [ {'on'} | 'off' ]     % 'on'  = show waitbar
125%                                               % 'off' = don't show waitbar
126%
127%        'Precision'   : [{'single'}|'double']  % Precision of the
128%                                               % outputted data.
129%                                               % 'single' = 32-bit float
130%                                               % 'double' = 64-bit float
131%
132%        'OrientImages': [ {'on'} | 'off' ]     % Orient FTDATA after
133%                                               % reading the data using
134%                                               % PROCPAR.orient property.
135%
136%        'RemoveEPIphaseIm' : [{'on'}|'off']    % Remove the phase image
137%                                               % from EPI data. This
138%                                               % option only affects EPI
139%                                               % images.
140%
141%        'EPIBlockSize' : [integer value]       % Block size (number of
142%                                               % volumes) used for
143%                                               % processing multireceiver
144%                                               % EPI data. Default=300
145%
146%        'EPIPhasedArrayData' : ['on'|{'off'}]  % Return data from
147%                                               % individual receivers from
148%                                               % phased array EPI files.
149%
150% Examples:
151%        [DATA,msg]=aedes_readvnmr(filename)             % Read image data from 'filename'
152%        [DATA,msg]=aedes_readvnmr(filename,'header')    % Read only data file header
153%
154%        [DATA,msg]=aedes_readvnmr(filename,'return',1)  % Return only image data (default)
155%        [DATA,msg]=aedes_readvnmr(filename,'return',2)  % Return only k-space
156%        [DATA,msg]=aedes_readvnmr(filename,'return',3)  % Return both image data and k-space
157%       
158%        % Read first data file header and then image data and k-space
159%        [DATA,msg]=aedes_readvnmr(filename,'header')
160%        [DATA,msg]=aedes_readvnmr(DATA.HDR,'return',3)
161%
162%        % Read VNMR-data with default options and save slices in NIfTI
163%        % format
164%        DATA=aedes_readvnmr('example.fid','saved_slices.nii');
165%
166%        % If you want to use non-default options and also write
167%        % NIfTI-files, use the property/value formalism, for example
168%        DATA=aedes_readvnmr('example.fid','ZeroPadding','off',...
169%                     'OutputFile','saved_slices.nii');
170%
171% See also:
172%        AEDES_READFIDPREFS, AEDES_READPROCPAR, AEDES_READPHASETABLE,
173%        AEDES_DATA_READ, AEDES_WRITE_NIFTI, AEDES
174
175% This function is a part of Aedes - A graphical tool for analyzing
176% medical images
177%
178% Copyright (C) 2010 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uef.fi>
179%
180% Department of Physics, Department of Neurobiology
181% University of Kuopio, FINLAND
182%
183% This program may be used under the terms of the GNU General Public
184% License version 2.0 as published by the Free Software Foundation
185% and appearing in the file LICENSE.TXT included in the packaging of
186% this program.
187%
188% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
189% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
190
191
192if nargout==2
193  Dat.ShowErrors = false;
194  msg_out='';
195else
196  Dat.ShowErrors = true;
197end
198
199%% ---------------------
200% Defaults
201Dat.ReturnKSpace  = false;
202Dat.ReturnFTData  = true;
203Dat.ReturnRawKspace = false;
204Dat.DCcorrection  = false;
205Dat.ZeroPadding = 2;
206Dat.Sorting = true;
207Dat.UsePhaseTable = true;
208Dat.FastDataRead = true;
209Dat.precision = 'single';
210Dat.FileChunkSize = 500; % Chunk size (in MB) for FastRead
211
212%% Other defaults
213Dat.ShowWaitbar   = true;
214procpar=[];
215Dat.phasetable=[];
216Dat.OutputFile = false;
217Dat.ReorientEPI = false;
218Dat.RemoveEPIphaseIm = true;
219Dat.EPIBlockSize = 100;
220Dat.EPIPhasedArrayData = false;
221Dat.OrientImages = true;
222Dat.SumOfSquares = 1;
223Dat.UseCustomRecon = true;
224
225% Default custom reconstruction file directory
226tmp = which('aedes');
227if ~isempty(tmp)
228  tmp_path=fileparts(tmp);
229  recon_dir = [tmp_path,filesep,'vnmr_recon',filesep];
230else
231  % If all else fails, look in the current directory
232  recon_dir = [pwd,filesep,'vnmr_recon',filesep];
233end
234
235%% Set data format label
236DATA.DataFormat = 'vnmr';
237
238%% Parse input arguments
239if nargin==0 || isempty(filename)
240 
241  %% Get default directory
242  try
243    defaultdir = getpref('Aedes','GetDataFileDir');
244  catch
245    defaultdir = [pwd,filesep];
246  end
247 
248  %% If no input arguments are given, ask for a file
249  [f_name, f_path, f_index] = uigetfile( ...
250       {'fid;FID','Varian FID-files (fid)'; ...
251        '*.*', 'All Files (*.*)'}, ...
252        'Select VNMR (Varian) FID file',defaultdir);
253  if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed
254    DATA=[];
255    msg_out='Canceled';
256    return
257  end
258  ReadHdr = true;
259  ReadData = true;
260  filename=[f_path,f_name];
261 
262  %% Set default directory to preferences
263  setpref('Aedes','GetDataFileDir',f_path)
264 
265end
266
267if nargin==1
268  if isstruct(filename)
269    hdr = filename;
270    filename = [hdr.fpath,hdr.fname];
271    ReadHdr = false;
272    ReadData = true;
273  elseif ischar(filename)
274    ReadHdr = true;
275    ReadData = true;
276  end
277elseif nargin==2
278  if strcmpi(varargin{1},'header')
279    ReadHdr = true;
280    ReadData = false;
281  elseif ischar(varargin{1})
282    ReadHdr = true;
283    ReadData = true;
284    Dat.OutputFile = varargin{1};
285  else
286    if ~Dat.ShowErrors
287      DATA=[];
288      msg_out=sprintf('Invalid second input argument "%s".',varargin{1});
289      return
290    else
291      error('Invalid second input argument "%s".',varargin{1})
292    end
293  end
294else
295  if isstruct(filename)
296    hdr = filename;
297    filename = [hdr.fpath,hdr.fname];
298    ReadHdr = false;
299    ReadData = true;
300  elseif isempty(filename)
301    [f_name, f_path, f_index] = uigetfile( ...
302        {'fid;FID','Varian FID-files (fid)'; ...
303         '*.*', 'All Files (*.*)'}, ...
304        'Select VNMR (Varian) FID file');
305    if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed
306      DATA=[];
307      msg_out='Canceled';
308      return
309    end
310    ReadHdr = true;
311    ReadData = true;
312    filename=[f_path,f_name];
313  else
314    ReadHdr = true;
315    ReadData = true;
316  end
317 
318  for ii=1:2:length(varargin)
319    switch lower(varargin{ii})
320      case 'return'
321        if length(varargin)<2
322          DATA=[];
323          if ~Dat.ShowErrors
324            msg_out='"Return" value not specified!';
325          else
326            error('"Return" value not specified!')
327          end
328          return
329        else
330          if varargin{ii+1}==1
331            Dat.ReturnKSpace = false;
332            Dat.ReturnFTData = true;
333          elseif varargin{ii+1}==2
334            Dat.ReturnKSpace = true;
335            Dat.ReturnFTData = false;
336          elseif varargin{ii+1}==3
337            Dat.ReturnKSpace = true;
338            Dat.ReturnFTData = true;
339          elseif varargin{ii+1}==4
340            Dat.ReturnRawKspace = true;
341            Dat.ReturnKSpace = true;
342            Dat.ReturnFTData = false;
343          end
344        end
345       
346      case {'dc','dccorrection','dccorr'}
347        if strcmpi(varargin{ii+1},'on')
348          Dat.DCcorrection = true;
349        else
350          Dat.DCcorrection = false;
351        end
352     
353      case 'procpar'
354        procpar=varargin{ii+1};
355       
356      case 'zeropadding'
357        if ischar(varargin{ii+1})
358          if strcmpi(varargin{ii+1},'on')
359            Dat.ZeroPadding = 1; % on
360          elseif strcmpi(varargin{ii+1},'off')
361            Dat.ZeroPadding = 0; % off
362        else
363          Dat.ZeroPadding = 2; % auto
364          end
365        else
366          % Undocumented
367          Dat.ZeroPadding = 3; % Custom
368          Dat.CustomZeroPadding = varargin{ii+1};
369        end
370       
371      case {'sumofsquares','soq'}
372        Dat.SumOfSquares = varargin{ii+1};
373         
374      case 'phasetable'
375        Dat.phasetable = varargin{ii+1};
376       
377      case 'sorting'
378        if strcmpi(varargin{ii+1},'on')
379          Dat.UsePhaseTable = true;
380          Dat.Sorting = true;
381        else
382          Dat.UsePhaseTable = false;
383          Dat.Sorting = false;
384        end
385       
386      case 'fastread'
387        if strcmpi(varargin{ii+1},'on')
388          Dat.FastDataRead = true;
389        else
390          Dat.FastDataRead = false;
391        end
392       
393      case 'filechunksize'
394        Dat.FileChunkSize = round(varargin{ii+1});
395       
396      case 'wbar'
397        if strcmpi(varargin{ii+1},'on')
398          Dat.ShowWaitbar = 1;
399        else
400          Dat.ShowWaitbar = 0;
401        end
402       
403      case 'precision'
404        if strcmpi(varargin{ii+1},'single')
405          Dat.precision = 'single';
406        elseif strcmpi(varargin{ii+1},'double')
407          Dat.precision = 'double';
408        else
409          warning('Unsupported precision "%s". Using single...',varargin{ii+1});
410          Dat.precision = 'single';
411        end
412       
413      case 'outputfile'
414        Dat.OutputFile = varargin{ii+1};
415       
416      case 'reorientepi'
417        if strcmpi(varargin{ii+1},'on')
418          Dat.ReorientEPI = true;
419        end
420       
421      case 'removeepiphaseim'
422        if strcmpi(varargin{ii+1},'on')
423          Dat.RemoveEPIphaseIm = true;
424        end
425      case 'epiblocksize'
426        Dat.EPIBlockSize = round(varargin{ii+1});
427       
428      case 'epiphasedarraydata'
429        if strcmpi(varargin{ii+1},'on')
430          Dat.EPIPhasedArrayData = true;
431        else
432          Dat.EPIPhasedArrayData = false;
433        end
434      case 'orientimages'
435        if strcmpi(varargin{ii+1},'off')
436          Dat.OrientImages = false;
437        end
438      otherwise
439        DATA=[];
440        if ~Dat.ShowErrors
441          msg_out = ['Unknown property "' varargin{ii} '"'];
442        else
443          error(['Unknown property "' varargin{ii} '"'])
444        end
445        return
446    end
447  end
448end
449
450
451
452% Parse filename
453[fpath,fname,fext]=fileparts(filename);
454if ~strcmpi(fname,'fid')
455  if isempty(fname)
456    fpath = [fpath,filesep];
457  else
458    if isempty(fpath)
459      fpath = [pwd,filesep,fname,fext,filesep];
460    else
461      fpath = [fpath,filesep,fname,fext,filesep];
462    end
463  end
464  fname = 'fid';
465else
466  fpath = [fpath,filesep];
467end
468
469%% Read procpar if not given as input argument
470if isempty(procpar) || nargin~=2
471  [procpar,proc_msg]=aedes_readprocpar([fpath,'procpar']);
472  if isempty(procpar)
473    DATA=[];
474    if ~Dat.ShowErrors
475      msg_out=proc_msg;
476    else
477      error(proc_msg)
478    end
479    return
480  end
481end
482
483% Don't use DC correction if number of transients is greater that 1
484if procpar.nt>1
485  Dat.DCcorrection = false;
486end
487
488%% Read phasetable if necessary
489if isfield(procpar,'petable') && isempty(Dat.phasetable) && ...
490    ~isempty(procpar.petable{1}) && ~strcmpi(procpar.petable{1},'n') && ...
491    Dat.Sorting
492  % Look in preferences for tablib-directory
493  try
494    tabpath = getpref('Aedes','TabLibDir');
495  catch
496    % If the table path was not found in preferences, try to use aedes
497    % directory
498    tmp_path = which('aedes');
499    if ~isempty(tmp_path)
500      aedes_path=fileparts(tmp_path);
501      tabpath = [aedes_path,filesep,'tablib',filesep];
502    else
503      % If all else fails, look in the current directory
504      tabpath = [pwd,filesep,'tablib',filesep];
505    end
506  end
507  [phasetable,msg] = aedes_readphasetable([tabpath,procpar.petable{1}]);
508 
509  % If petable cannot be found, check if it is in procpar...
510  if isempty(phasetable) && isfield(procpar,'pe_table')
511    phasetable = {'t1',procpar.pe_table(:)};
512  elseif isempty(phasetable) && isfield(procpar,'pelist') && ...
513      ~isempty(procpar.pelist) && isnumeric(procpar.pelist) && ...
514      length(procpar.pelist)>1
515    phasetable = {'t1',reshape(procpar.pelist,procpar.etl,[]).'};
516  end
517 
518  %% Abort and throw error if phasetable cannot be read and the FID-file
519  % has not been sorted
520  if isempty(phasetable) && ~isfield(procpar,'flash_converted')
521    DATA=[];
522    if ~Dat.ShowErrors
523      msg_out={['Could not find the required phase table "' ...
524                procpar.petable{1} '" in the following folders'],...
525               ['"' tabpath '".']};
526    else
527      error('Could not find the required phase table "%s" in %s',...
528            procpar.petable{1},tabpath)
529    end
530    return
531  elseif ( isempty(phasetable) && isfield(procpar,'flash_converted') ) || ...
532      ~Dat.UsePhaseTable
533    %% If the FID-file has been sorted, the reading can continue but
534    % throw a warning.
535    fprintf(1,'Warning: aedes_readfid: Could not find phasetable "%s" in "%s"!\n',procpar.petable{1},tabpath)
536    Dat.phasetable=[];
537  else
538    Dat.phasetable = phasetable{1,2};
539  end
540end
541
542% Convert phasetable indices to MATLAB indices
543if ~isempty(Dat.phasetable)
544  Dat.phasetable=Dat.phasetable+(-min(min(Dat.phasetable))+1);
545else
546  Dat.UsePhaseTable = false;
547end
548
549%% Open FID-file
550[file_fid,msg]=fopen([fpath,fname],'r','ieee-be');
551if file_fid < 0
552  DATA=[];
553  if ~Dat.ShowErrors
554    msg_out=['Could not open file "' filename '" for reading.'];
555  else
556    error('Could not open file "%s" for reading.',filename)
557  end
558  return
559end
560
561% Read only header
562if ReadHdr && ~ReadData
563  [hdr,msg_out]=l_ReadDataFileHeader(file_fid);
564  if ~isempty(msg_out)
565    DATA=[];
566    fclose(file_fid);
567    if Dat.ShowErrors
568      error(msg_out)
569    end
570    return
571  else
572    DATA.HDR.FileHeader = hdr.FileHeader;
573    DATA.FTDATA=[];
574    DATA.KSPACE=[];
575    DATA.PROCPAR=[];
576    DATA.PHASETABLE=[];
577  end
578elseif ~ReadHdr && ReadData % Header structure given as input argument
579                            % Read only data.
580  [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar);
581  if ~isempty(msg_out)
582    DATA=[];
583    fclose(file_fid);
584    if Dat.ShowErrors
585      error(msg_out)
586    end
587    return
588  else
589    DATA.HDR.FileHeader=hdr.FileHeader;
590    DATA.HDR.BlockHeaders = hdr.BlockHeaders;
591    DATA.FTDATA=data;
592    DATA.KSPACE=kspace;
593    DATA.PROCPAR=procpar;
594    DATA.PHASETABLE=Dat.phasetable;
595  end
596elseif ReadHdr && ReadData  % Read headers and data
597  [hdr,msg_out]=l_ReadDataFileHeader(file_fid);
598  [hdr,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat);
599  if ~isempty(msg_out)
600    DATA=[];
601    fclose(file_fid);
602    if Dat.ShowErrors
603      error(msg_out)
604    end
605    return
606  else
607    DATA.HDR.FileHeader=hdr.FileHeader;
608    DATA.HDR.BlockHeaders = hdr.BlockHeaders;
609    DATA.FTDATA=[];
610    DATA.KSPACE=[];
611    DATA.PROCPAR=procpar;
612    DATA.PHASETABLE=Dat.phasetable;
613  end
614end
615
616% Close file
617fclose(file_fid);
618
619
620% Set file name and path to the HDR structure
621DATA.HDR.fname=fname;
622DATA.HDR.fpath=fpath;
623
624% Set parameter values
625if ~Dat.DCcorrection
626  DATA.HDR.Param.DCcorrection = 'off';
627else
628  DATA.HDR.Param.DCcorrection = 'on';
629end
630
631% Reconstruct kspace =======================================
632if Dat.ReturnRawKspace
633  % If asked to return raw unsorted kspace, return immediately
634  DATA.KSPACE=kspace;
635  return
636else
637 
638  % Check if the sequence used has a custom reconstruct code
639  recon_func=l_GetReconFunc(recon_dir);
640  recon_func_ind = 0;
641  if ~isempty(recon_func)
642    % Get sequence name
643    seqfil = procpar.seqfil;
644   
645    % Check if any of the custom reconstruction files would
646    % like to do the reconstruction...
647    for ii=1:length(recon_func)
648      try
649        list=recon_func{ii}();
650      catch
651        warning('Error in custom reconsruction function "%s", skipping...',func2str(recon_func{ii}));
652        continue
653      end
654      if any(strcmp(seqfil,list))
655        recon_func_ind = ii;
656        break
657      end
658    end
659  end
660 
661  if recon_func_ind==0
662    Dat.UseCustomRecon = false;
663  end
664 
665  if Dat.UseCustomRecon
666    if Dat.ShowWaitbar
667      wbh=aedes_calc_wait(sprintf('%s\n%s',...
668        ['Using custom function ',upper(func2str(recon_func{recon_func_ind}))],...
669        ['to reconstruct sequence ',procpar.seqfil{1}]));
670      drawnow
671    end
672    [kspace,data,msg_out]=recon_func{recon_func_ind}(kspace,Dat,procpar);
673    if Dat.ShowWaitbar
674      close(wbh)
675    end
676    if isempty(data) && Dat.ReturnFTData
677      % Fourier transform data if not done in custom reconstruction code
678      [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar);
679    end
680  else
681    % Use the default reconstruction code
682    [kspace,msg_out]=l_ReconstructKspace(kspace,Dat,procpar);
683   
684    % Fourier transform data
685    if Dat.ReturnFTData
686      [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar);
687    end
688  end
689end
690if Dat.ReturnKSpace
691  DATA.KSPACE = kspace;
692else
693  DATA.KSPACE = [];
694end
695if Dat.ReturnFTData
696  DATA.FTDATA = data;
697else
698  DATA.FTDATA = [];
699end
700clear kspace data
701
702
703return
704
705%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
706% Read Data File (Main) Header
707%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
708function [hdr,msg_out]=l_ReadDataFileHeader(file_fid)
709
710msg_out='';
711%% Read Data File Header
712try
713  FH.nblocks = fread(file_fid,1,'long');   % Number of blocks in file
714  FH.ntraces = fread(file_fid,1,'long');   % Number of traces per block
715  FH.np = fread(file_fid,1,'long');        % Number of elements per trace
716  FH.ebytes = fread(file_fid,1,'long');    % Number of bytes per element
717  FH.tbytes = fread(file_fid,1,'long');    % Number of bytes per trace
718  FH.bbytes = fread(file_fid,1,'long');    % Number of bytes per block
719  FH.vers_id = fread(file_fid,1,'short');  % Software version, file_id status bits
720  FH.status = fread(file_fid,1,'short');   % Status of whole file
721  FH.nbheaders = fread(file_fid,1,'long'); % Number of block headers per block
722 
723  hdr.FileHeader = FH;
724 
725  %% Parse status bits
726  hdr.FileHeader.status=[];
727  tmp_str = {'S_DATA',...          % 0 = no data, 1 = data
728             'S_SPEC',...          % 0 = FID, 1 = spectrum
729             'S_32',...            % 0 = 16 bit, 1 = 32 bit integer
730             'S_FLOAT',...         % 0 = integer, 1 = floating point
731             'S_COMPLEX',...       % 0 = real, 1 = complex
732             'S_HYPERCOMPLEX',...  % 1 = hypercomplex
733             '',...                % Unused bit
734             'S_ACQPAR',...        % 0 = not Acqpar, 1 = Acqpar
735             'S_SECND',...         % 0 = first FT, 1 = second FT
736             'S_TRANSF',...        % 0 = regular, 1 = transposed
737             '',...                % Unused bit
738             'S_NP',...            % 1 = np dimension is active
739             'S_NF',...            % 1 = nf dimension is active
740             'S_NI',...            % 1 = ni dimension is active
741             'S_NI2',...           % 1 = ni2 dimension is active
742             ''...                 % Unused bit
743            };
744  status_bits = fliplr(double(dec2bin(FH.status,16))==49);
745  for ii=1:length(tmp_str)
746    if ~isempty(tmp_str{ii})
747      hdr.FileHeader.status.(tmp_str{ii}) = status_bits(ii);
748    end
749  end
750catch
751  hdr=[];
752  msg_out=['Error occurred while reading file header from "' filename '".'];
753  return
754end
755
756
757
758%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
759% Read Block Headers and Data
760%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
761function [hdr,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat)
762
763hdr.BlockHeaders = [];
764%hdr.HypercmplxBHeaders = [];
765kspace=[];
766count = [];
767msg_out='';
768
769
770BlockHeadStatusLabels = {...
771  'S_DATA',...       % 0 = no data, 1 = data
772  'S_SPEC',...          % 0 = FID, 1 = spectrum
773  'S_32',...            % 0 = 16 bit, 1 = 32 bit integer
774  'S_FLOAT',...         % 0 = integer, 1 = floating point
775  'S_COMPLEX',...       % 0 = real, 1 = complex
776  'S_HYPERCOMPLEX',...  % 1 = hypercomplex
777  '',...                % Unused bit
778  'MORE_BLOCKS',...     % 0 = absent, 1 = present
779  'NP_CMPLX',...        % 0 = real, 1 = complex
780  'NF_CMPLX',...        % 0 = real, 1 = complex
781  'NI_CMPLX',...        % 0 = real, 1 = complex
782  'NI2_CMPLX',...       % 0 = real, 1 = complex
783  '',...                % Unused bit
784  '',...                % Unused bit
785  '',...                % Unuesd bit
786  ''...                 % Unused bit
787  };
788
789BlockHeadModeLabels = {...
790  'NP_PHMODE',...   % 1 = ph mode
791  'NP_AVMODE',...    % 1 = av mode
792  'NP_PWRMODE',...   % 1 = pwr mode
793  '',...             % Unused bit
794  'NF_PHMODE',...    % 1 = ph mode
795  'NF_AVMODE',...    % 1 = av mode
796  'NF_PWRMODE',...   % 1 = pwr mode
797  '',...             % Unused bit
798  'NI_PHMODE',...    % 1 = ph mode
799  'NI_AVMODE',...    % 1 = av mode
800  'NI_PWRMODE',...   % 1 = pwr mode
801  '',...             % Unused bit
802  'NI2_PHMODE',...   % 1 = ph mode
803  'NI2_AVMODE',...   % 1 = av mode
804  'NI2_PWRMODE',...  % 1 = pwr mode
805  ''...              % Unused bit
806  };
807
808
809% The nbheaders -field is not correct in some cases. Thus, this field
810% cannot be trusted and the real nbheaders has to be calculated from the
811% byte counts.
812tmp_bytes=hdr.FileHeader.bbytes-hdr.FileHeader.tbytes*hdr.FileHeader.ntraces;
813header_bytes=28;
814if rem(tmp_bytes,header_bytes)~=0
815  msg_out = 'Block header byte count does not match with file header';
816  return
817else
818  nbheaders = tmp_bytes/28;
819end
820
821% Allocate space for kspace
822% kspace = complex(zeros(hdr.FileHeader.np/2,...
823%           hdr.FileHeader.ntraces*hdr.FileHeader.nblocks,...
824%           Dat.precision));
825
826kspace = complex(zeros(hdr.FileHeader.np/2,...
827        hdr.FileHeader.ntraces*hdr.FileHeader.nblocks,Dat.precision));
828
829
830%% - The older robust (but also slower) way for reading the data.
831%% When the blocksize is relatively small, this is also quite fast.
832if ~Dat.FastDataRead
833 
834  % Initialize waitbar
835  if Dat.ShowWaitbar
836    wb_h = aedes_wbar(0/hdr.FileHeader.nblocks,...
837      {['Reading VNMR data.'],...
838      ['(Processing data block ' ...
839      num2str(0) '/' num2str(hdr.FileHeader.nblocks) ')']});
840  end
841
842  %% Read data and block headers
843  for ii=1:hdr.FileHeader.nblocks
844    %% Update waitbar
845    if Dat.ShowWaitbar
846      aedes_wbar(ii/hdr.FileHeader.nblocks,...
847        wb_h,...
848        {['Reading VNMR data.'],...
849        ['(Processing data block ' ...
850        num2str(ii) '/' num2str(hdr.FileHeader.nblocks) ')']})
851    end
852
853    %% Read block header and hypercomplex header
854    for kk=1:nbheaders
855      %% Read block header
856      if kk==1
857        hdr.BlockHeaders.scale = fread(file_fid,1,'short');  % Scaling factor
858        tmp_status = fread(file_fid,1,'short'); % Status of data in block
859        hdr.BlockHeaders.status = [];
860        hdr.BlockHeaders.index = fread(file_fid,1,'short');  % Block index
861        tmp_mode = fread(file_fid,1,'short');   % Mode of data in block
862        hdr.BlockHeaders.mode = [];
863        hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
864        hdr.BlockHeaders.lpval = fread(file_fid,1,'float');  % f2 (2D-f1) left phase in phasefile
865        hdr.BlockHeaders.rpval = fread(file_fid,1,'float');  % f2 (2D-f1) right phase in phasefile
866        hdr.BlockHeaders.lvl = fread(file_fid,1,'float');    % level drift correction
867        hdr.BlockHeaders.tlt = fread(file_fid,1,'float');    % tilt drift correction
868
869        %% Parse status and mode bits
870        status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
871        mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
872        for tt=1:length(BlockHeadStatusLabels)
873          if ~isempty(BlockHeadStatusLabels{tt})
874            hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
875          end
876          if ~isempty(BlockHeadModeLabels{tt})
877            hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
878          end
879        end
880
881
882      else %% Read hypercomplex header
883        fread(file_fid,1,'short'); % Spare
884        hdr.BlockHeaders.HCHstatus = fread(file_fid,1,'short');
885        fread(file_fid,1,'short'); % Spare
886        fread(file_fid,1,'short'); % Spare
887        fread(file_fid,1,'long'); % Spare
888        hdr.BlockHeaders.HCHlpval1 = fread(file_fid,1,'float');
889        hdr.BlockHeaders.HCHrpval1 = fread(file_fid,1,'float');
890        fread(file_fid,1,'float'); % Spare
891        fread(file_fid,1,'float'); % Spare
892      end
893    end
894
895    %% Check block index to be sure about the data type
896    if hdr.BlockHeaders.index~=ii
897      fprintf(1,'Warning: Index mismatch in "%s" block %d\n',fopen(file_fid),ii);
898
899      % Use information from the file header instead of the BlockHeader if
900      % there is a mismatch in blockheader index...
901      useFileHeader = true;
902    else
903      useFileHeader = false;
904    end
905
906    %% Determine data precision
907    if useFileHeader
908      if hdr.FileHeader.status.S_FLOAT==1
909        prec_str = ['single=>',Dat.precision];
910      elseif hdr.FileHeader.status.S_32==1 ...
911          && hdr.FileHeader.status.S_FLOAT==0
912        prec_str = ['int32=>',Dat.precision];
913      elseif hdr.FileHeader.status.S_32==0 ...
914          && hdr.FileHeader.status.S_FLOAT==0
915        prec_str = ['int16=>',Dat.precision];
916      end
917
918    else
919      if hdr.BlockHeaders.status.S_FLOAT==1
920        prec_str = ['single=>',Dat.precision];
921      elseif hdr.BlockHeaders.status.S_32==1 ...
922          && hdr.BlockHeaders.status.S_FLOAT==0
923        prec_str = ['int32=>',Dat.precision];
924      elseif hdr.BlockHeaders.status.S_32==0 ...
925          && hdr.BlockHeaders.status.S_FLOAT==0
926        prec_str = ['int16=>',Dat.precision];
927      end
928    end
929
930    % Read k-space
931    tmp=fread(file_fid,...
932      [hdr.FileHeader.np,hdr.FileHeader.ntraces],...
933      prec_str);
934   
935
936    % Store complex block and perform DC correction
937    if ~Dat.DCcorrection
938      complex_block = complex(tmp(1:2:end,:),tmp(2:2:end,:));
939    else
940      complex_block = complex(tmp(1:2:end,:)-hdr.BlockHeaders.lvl,...
941        tmp(2:2:end,:)-hdr.BlockHeaders.tlt);
942    end
943   
944    inds = ((ii-1)*hdr.FileHeader.ntraces+1):(ii*hdr.FileHeader.ntraces);
945    kspace(:,inds) = complex_block;
946
947    % Do not save blockheaders by default. When reading data files with a lot of
948    % blocks (e.g. over 1000) the processing of the DATA structure can be
949    % slowed down considerably. If you for some reason want to save also the
950    % block headers in the DATA structure comment out the code line below.
951    hdr.BlockHeaders = [];
952  end % for ii=1:hdr.
953
954  if Dat.ShowWaitbar
955    delete(wb_h)
956  end
957
958else
959  %% -------------------------------------------------------------------
960  %% Fast Method for reading data. This may fail with some datas and can
961  %% also require a relatively large amount of memory. This
962  %% method should be used for EPI datas that contain a large number
963  %% of block headers...
964
965  % Check the size of the FID-file
966  d=dir(fopen(file_fid));
967  file_sz = d.bytes/1024/1024; % File size in MB
968  if file_sz<Dat.FileChunkSize
969    nBlocks = 1;
970  else
971    nBlocks = ceil(file_sz/Dat.FileChunkSize); % Read data in 500 MB blocks
972  end
973
974  % Initialize waitbar
975  if Dat.ShowWaitbar
976    if nBlocks==1
977      wb_h = aedes_calc_wait(['Reading VNMR data.']);
978    else
979      wb_h = aedes_wbar(1/nBlocks,{['Reading VNMR data.'],...
980        sprintf('Reading block 0/%d',nBlocks)});
981    end
982  end
983
984  % The first block header is read and it is assumed that the values in
985  % the other block headers don't change.
986  hdr.BlockHeaders.scale = fread(file_fid,1,'short');  % Scaling factor
987  tmp_status = fread(file_fid,1,'short'); % Status of data in block
988  hdr.BlockHeaders.status = [];
989  hdr.BlockHeaders.index = fread(file_fid,1,'short');  % Block index
990  tmp_mode = fread(file_fid,1,'short');   % Mode of data in block
991  hdr.BlockHeaders.mode = [];
992  hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
993  hdr.BlockHeaders.lpval = fread(file_fid,1,'float');  % f2 (2D-f1) left phase in phasefile
994  hdr.BlockHeaders.rpval = fread(file_fid,1,'float');  % f2 (2D-f1) right phase in phasefile
995  hdr.BlockHeaders.lvl = fread(file_fid,1,'float');    % level drift correction
996  hdr.BlockHeaders.tlt = fread(file_fid,1,'float');    % tilt drift correction
997
998  %% Parse status and mode bits
999  status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
1000  mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
1001  for tt=1:length(BlockHeadStatusLabels)
1002    if ~isempty(BlockHeadStatusLabels{tt})
1003      hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
1004    end
1005    if ~isempty(BlockHeadModeLabels{tt})
1006      hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
1007    end
1008  end
1009
1010  %% Determine data precision
1011  if hdr.BlockHeaders.status.S_FLOAT==1
1012    prec_str = ['single=>',Dat.precision];
1013    prec = 4; % Precision in bytes
1014  elseif hdr.BlockHeaders.status.S_32==1 ...
1015      && hdr.BlockHeaders.status.S_FLOAT==0
1016    prec_str = ['int32=>',Dat.precision];
1017    prec = 4;
1018  elseif hdr.BlockHeaders.status.S_32==0 ...
1019      && hdr.BlockHeaders.status.S_FLOAT==0
1020    prec_str = ['int16=>',Dat.precision];
1021    prec = 2;
1022  end
1023
1024  % Seek the file back to the beginning of the first block header
1025  fseek(file_fid,-28,0);
1026
1027  % Determine the number of values that will result from block header(s)
1028  nVals = (nbheaders*28)/prec;
1029
1030  nbh = floor(hdr.FileHeader.nblocks/nBlocks);
1031  szh = nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces;
1032  for ii=1:nBlocks
1033    if nBlocks~=1
1034      aedes_wbar(ii/nBlocks,wb_h,{['Reading VNMR data.'],...
1035        sprintf('Reading block %d/%d',ii,nBlocks)});
1036    end
1037
1038    % Read the whole data including block headers etc...
1039    if ii==nBlocks
1040      tmp = fread(file_fid,inf,prec_str);
1041    else
1042      tmp = fread(file_fid,nbh*szh,prec_str);
1043    end
1044    tmp=reshape(tmp,nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces,[]);
1045    tmp(1:nVals,:)=[];
1046    tmp=reshape(tmp,hdr.FileHeader.np,[]);
1047
1048    if ii==nBlocks
1049      inds = ((ii-1)*nbh*hdr.FileHeader.ntraces+1):size(kspace,2);
1050    else
1051      inds = ((ii-1)*nbh*hdr.FileHeader.ntraces+1):ii*nbh*hdr.FileHeader.ntraces;
1052    end
1053
1054    % Do DC-correction if necessary
1055    if ~Dat.DCcorrection
1056      kspace(:,inds)=complex(tmp(1:2:end,:,:),tmp(2:2:end,:,:));
1057    else
1058      kspace(:,inds)=complex(tmp(1:2:end,:,:)-hdr.BlockHeaders.lvl,...
1059        tmp(2:2:end,:,:)-hdr.BlockHeaders.tlt);
1060    end
1061  end
1062  clear tmp
1063
1064  if Dat.ShowWaitbar
1065    delete(wb_h)
1066    pause(0.1)
1067  end
1068
1069end
1070
1071% ==================================================================
1072% Reconstruct k-space
1073% ==================================================================
1074function [kspace,msg_out]=l_ReconstructKspace(kspace,Dat,procpar)
1075
1076msg_out = '';
1077
1078% Get the strategic parameters
1079seqcon = procpar.seqcon{1};
1080seqfil = procpar.seqfil{1};
1081np = procpar.np;
1082nv = procpar.nv;
1083nv2 = procpar.nv2;
1084nv3 = procpar.nv3;
1085ns = procpar.ns;
1086if isfield(procpar,'etl')
1087  etl = procpar.etl;
1088else
1089  etl = 1;
1090end
1091pss = procpar.pss;
1092nt = procpar.nt;
1093nf = procpar.nf;
1094ne = procpar.ne;
1095if isfield(procpar,'flash_converted')
1096  % Don't try to sort already sorted data...
1097  Dat.Sorting = false;
1098  seqcon(2:3) = 'sc';
1099end
1100
1101% Check dimensions
1102if nv==0
1103  % 1D-image (specrum)
1104  AcqType = 1;
1105elseif nv2==0
1106  % 2D-image
1107  AcqType = 2;
1108elseif nv3==0
1109  % 3D-image
1110  AcqType = 3;
1111else
1112  AcqType = 4;
1113end
1114
1115% Check number of receivers
1116if isfield(procpar,'rcvrs') && length(procpar.rcvrs{1})>1
1117  nRcvrs = length(find(procpar.rcvrs{1}=='y'));
1118else
1119  nRcvrs = 1;
1120end
1121
1122
1123% Check for arrayed acquisition
1124if isfield(procpar,'array') && ~isempty(procpar.array{1})
1125  %if length(procpar.array)==1 && ~iscell(procpar.array{1}) && ...
1126  %    strcmp(procpar.array{1},'pad') && all(procpar.pad==0)
1127  %  % Skip the array parsing if the array is a dummy "pad"-array...
1128  %  isAcqArrayed = false;
1129  %  ArrayLength = 1;
1130  %else
1131    isAcqArrayed = true;
1132    ArrayLength = [];
1133   
1134    % Determine array length
1135    for ii=1:length(procpar.array)
1136      if iscell(procpar.array{ii})
1137        ArrayLength(ii) = length(procpar.(procpar.array{ii}{1}));
1138      else
1139        ArrayLength(ii) = length(procpar.(procpar.array{ii}));
1140      end
1141    end
1142    ArrayLength = prod(ArrayLength);
1143  %end
1144else
1145  isAcqArrayed = false;
1146  ArrayLength = 1;
1147end
1148
1149
1150% Reconstruct k-space ----------------------
1151if AcqType==1
1152  % Reconstruct 1D data ...
1153elseif AcqType==2
1154 
1155  % Reconstruct 2D data
1156  if strcmpi(seqcon,'ncsnn')
1157    kspace = reshape(kspace,[np/2,etl,ns,ArrayLength,nv/etl,nRcvrs]);
1158    kspace = permute(kspace,[1 2 5 3 4 6]);
1159    kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1160  elseif strcmpi(seqcon,'nscnn')
1161    if isfield(procpar,'flash_converted')
1162      kspace = reshape(kspace,[np/2,nRcvrs,ArrayLength,nv,ns]);
1163      kspace = permute(kspace,[1 4 5 3 2]);
1164      kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1165    else
1166      kspace = reshape(kspace,[np/2,nRcvrs,nv,ArrayLength,ns]);
1167      kspace = permute(kspace,[1 3 5 4 2]);
1168      kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1169    end
1170  elseif strcmpi(seqcon,'nccnn')
1171    kspace = reshape(kspace,[np/2,etl,nRcvrs,ns,nv/etl,ArrayLength]);
1172    kspace = permute(kspace,[1 2 5 4 6 3]);
1173    kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1174  end
1175elseif AcqType==3
1176  % Reconstruct 3D data
1177  if strcmpi(seqcon,'nscsn')
1178    kspace = reshape(kspace,[np/2,etl,nv/etl,nRcvrs,ArrayLength*nv2]);
1179    kspace = permute(kspace,[1 2 3 5 4]);
1180    kspace = reshape(kspace,[np/2,nv,nv2,ArrayLength,nRcvrs]);
1181  end
1182else
1183  % Reconstruct nD data ...
1184end
1185
1186
1187% Sort data using phasetable ------------------
1188if Dat.Sorting && ~isempty(Dat.phasetable)
1189  Dat.phasetable = Dat.phasetable.';
1190  kspace(:,Dat.phasetable(:),:,:,:)=kspace;
1191end
1192
1193
1194
1195% ==================================================================
1196% Fourier transform k-space
1197% ==================================================================
1198function [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar)
1199
1200msg_out = '';
1201data = [];
1202
1203
1204% Check dimensions
1205if procpar.nv==0
1206  % 1D-image
1207  AcqType = 1;
1208elseif procpar.nv2==0
1209  % 2D-image
1210  AcqType = 2;
1211elseif procpar.nv3==0
1212  % 3D-image
1213  AcqType = 3;
1214else
1215  AcqType = 4;
1216end
1217
1218% Check number of receivers
1219if isfield(procpar,'rcvrs') && length(procpar.rcvrs{1})>1
1220  nRcvrs = length(find(procpar.rcvrs{1}=='y'));
1221else
1222  nRcvrs = 1;
1223end
1224
1225% If zeropadding is requested, calculate the padded size
1226if Dat.ZeroPadding~=0
1227  if Dat.ZeroPadding==1
1228    % Zeropad to square
1229    if AcqType==1
1230      padSize = procpar.np/2;
1231    elseif AcqType==2
1232      padSize = ones(1,2)*procpar.np/2;
1233      padSize(3) = size(kspace,3);
1234    else
1235      padSize = ones(1,3)*procpar.np/2;
1236    end
1237  else
1238    % Zeropadding is on "auto", i.e. zeropad to FOV
1239    lpe = procpar.lpe;
1240    lpe2 = procpar.lpe2;
1241    lro = procpar.lro;
1242    if AcqType==2
1243      % 2D data
1244      padSize = [procpar.np/2 ...
1245        procpar.np/2*(lpe/lro) ...
1246        size(kspace,3)];
1247    elseif AcqType==3 && lpe2~=0
1248      % 3D data
1249      padSize = [procpar.np/2 ...
1250        procpar.np/2*(lpe/lro) ...
1251        procpar.np/2*(lpe2/lro)];
1252    end
1253  end
1254  ks_sz = [size(kspace,1) ...
1255    size(kspace,2) ...
1256    size(kspace,3)];
1257  if any(padSize>ks_sz)
1258    kspace(padSize) = 0;
1259  end
1260else
1261  padSize = [size(kspace,1) ...
1262    size(kspace,2) ...
1263    size(kspace,3)];
1264end
1265
1266% Allocate space for Fourier transformed data
1267if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1268  data_sz = [padSize,size(kspace,4),size(kspace,5)+1];
1269  data = zeros(data_sz,Dat.precision);
1270else
1271  data = zeros(size(kspace),Dat.precision);
1272end
1273%data = [];
1274%if strcmpi(Dat.precision,'single')
1275%  data = single(data);
1276%end
1277
1278% Fourier transform data
1279if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1280  ind = [2:size(data,5)];
1281else
1282  ind = [1:size(data,5)];
1283end
1284if AcqType==1
1285  data(:,:,:,:,ind) = abs(fftshift(ifft(kspace,[],1),1));
1286elseif AcqType==2
1287  data(:,:,:,:,ind) = abs(fftshift(fftshift(ifft(ifft(kspace,[],1),[],2),1),2));
1288elseif AcqType==3
1289  data(:,:,:,:,ind) = abs(fftshift(fftshift(fftshift(ifft(ifft(ifft(kspace,[],1),[],2),[],3),1),2),3));
1290end
1291
1292% Calculate sum-of-squares image
1293if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1294  % Calculate sum-of-squares
1295  data(:,:,:,:,1) = sqrt(sum(data(:,:,:,:,ind).^2,5));
1296  data=abs(data);
1297  if Dat.SumOfSquares==1
1298    % Remove individual receiver data
1299    data=data(:,:,:,:,1);
1300  end
1301end
1302
1303% ==================================================================
1304% Find custom functions for VNMR k-space reconstruction
1305% ==================================================================
1306function recon_func=l_GetReconFunc(recon_dir)
1307
1308recon_func = {};
1309
1310if ~isdir(recon_dir)
1311  return
1312end
1313
1314dir_struct=dir(recon_dir);
1315recon_files = {dir_struct(~[dir_struct(:).isdir]).name};
1316if isempty(recon_files)
1317  return
1318end
1319
1320% Remove files that don't have .m extension
1321ind = regexpi(recon_files,'\.m$');
1322recon_files = {recon_files{~cellfun('isempty',ind)}};
1323
1324currentDir = pwd;
1325try
1326  cd(recon_dir);
1327  for ii=1:length(recon_files)
1328    recon_func{ii}=str2func(recon_files{ii}(1:end-2));
1329  end
1330  cd(currentDir);
1331catch
1332  cd(currentDir);
1333end
1334
1335
1336
Note: See TracBrowser for help on using the repository browser.

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