source: aedes_readvnmr.m @ 178

Last change on this file since 178 was 177, checked in by tjniskan, 8 years ago
  • Added support for reading complex data from bruker 2dseq files.
  • Fixed scaling of 3D data from 2dseq files.
  • Aedes now prompts for handling of complex data (absolute, real

or imag parts) if only complex data is read from a file.

  • Added an option for EPI phase correction in aedes_readvnmr and epi_recon

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

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

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