source: aedes_readvnmr.m @ 202

Last change on this file since 202 was 202, checked in by tjniskan, 7 years ago
  • Fixed a bug which prevented pelist from being used be aedes_readvnmr

M aedes_readvnmr.m
M aedes_revision.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 Dat.Sorting
523       
524  % Look in preferences for tablib-directory
525  if isfield(procpar,'petable') && ~strcmpi(procpar.petable{1},'n')
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        else
542                phasetable = [];
543        end
544               
545  % If petable cannot be found, check if it is in procpar...
546  if isempty(phasetable) && isfield(procpar,'pe_table')
547    phasetable = {'t1',procpar.pe_table(:)};
548  elseif isempty(phasetable) && isfield(procpar,'pelist') && ...
549      ~isempty(procpar.pelist) && isnumeric(procpar.pelist) && ...
550      length(procpar.pelist)>1
551    phasetable = {'t1',reshape(procpar.pelist,procpar.etl,[]).'};
552  end
553 
554  %% Abort and throw error if phasetable cannot be read and the FID-file
555  % has not been sorted
556  if isempty(phasetable) && ~isfield(procpar,'flash_converted')
557    DATA=[];
558    if ~Dat.ShowErrors
559      msg_out={['Could not find the required phase table "' ...
560                procpar.petable{1} '" in the following folders'],...
561               ['"' tabpath '".']};
562    else
563      error('Could not find the required phase table "%s" in %s',...
564            procpar.petable{1},tabpath)
565    end
566    return
567  elseif ( isempty(phasetable) && isfield(procpar,'flash_converted') ) || ...
568      ~Dat.UsePhaseTable
569    %% If the FID-file has been sorted, the reading can continue but
570    % throw a warning.
571    fprintf(1,'Warning: aedes_readfid: Could not find phasetable "%s" in "%s"!\n',procpar.petable{1},tabpath)
572    Dat.phasetable=[];
573  else
574    Dat.phasetable = phasetable{1,2};
575  end
576end
577
578% Convert phasetable indices to MATLAB indices
579if ~isempty(Dat.phasetable)
580  Dat.phasetable=Dat.phasetable+(-min(min(Dat.phasetable))+1);
581else
582  Dat.UsePhaseTable = false;
583end
584
585%% Open FID-file
586[file_fid,msg]=fopen([fpath,fname],'r','ieee-be');
587if file_fid < 0
588  DATA=[];
589  if ~Dat.ShowErrors
590    msg_out=['Could not open file "' filename '" for reading.'];
591  else
592    error('Could not open file "%s" for reading.',filename)
593  end
594  return
595end
596
597% Read only header
598if ReadHdr && ~ReadData
599  [hdr,msg_out]=l_ReadDataFileHeader(file_fid);
600  if ~isempty(msg_out)
601    DATA=[];
602    fclose(file_fid);
603    if Dat.ShowErrors
604      error(msg_out)
605    end
606    return
607  else
608    DATA.HDR.FileHeader = hdr.FileHeader;
609    DATA.FTDATA=[];
610    DATA.KSPACE=[];
611    DATA.PROCPAR=[];
612    DATA.PHASETABLE=[];
613  end
614elseif ~ReadHdr && ReadData % Header structure given as input argument
615                            % Read only data.
616  [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar);
617  if ~isempty(msg_out)
618    DATA=[];
619    fclose(file_fid);
620    if Dat.ShowErrors
621      error(msg_out)
622    end
623    return
624  else
625    DATA.HDR.FileHeader=hdr.FileHeader;
626    DATA.HDR.BlockHeaders = hdr.BlockHeaders;
627    DATA.FTDATA=data;
628    DATA.KSPACE=kspace;
629    DATA.PROCPAR=procpar;
630    DATA.PHASETABLE=Dat.phasetable;
631  end
632elseif ReadHdr && ReadData  % Read headers and data
633  [hdr,msg_out]=l_ReadDataFileHeader(file_fid);
634  [hdr,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat);
635  if ~isempty(msg_out)
636    DATA=[];
637    fclose(file_fid);
638    if Dat.ShowErrors
639      error(msg_out)
640    end
641    return
642  else
643    DATA.HDR.FileHeader=hdr.FileHeader;
644    DATA.HDR.BlockHeaders = hdr.BlockHeaders;
645    DATA.FTDATA=[];
646    DATA.KSPACE=[];
647    DATA.PROCPAR=procpar;
648    DATA.PHASETABLE=Dat.phasetable;
649  end
650end
651
652% Close file
653fclose(file_fid);
654
655
656% Set file name and path to the HDR structure
657DATA.HDR.fname=fname;
658DATA.HDR.fpath=fpath;
659
660% Set parameter values
661if ~Dat.DCcorrection
662  DATA.HDR.Param.DCcorrection = 'off';
663else
664  DATA.HDR.Param.DCcorrection = 'on';
665end
666
667% Reconstruct kspace =======================================
668if Dat.ReturnRawKspace
669  % If asked to return raw unsorted kspace, return immediately
670  DATA.KSPACE=kspace;
671  return
672else
673 
674  % Check if the sequence used has a custom reconstruct code
675  recon_func=l_GetReconFunc(recon_dir);
676  recon_func_ind = 0;
677  if ~isempty(recon_func)
678    % Get sequence name
679    seqfil = procpar.seqfil;
680   
681    % Check if any of the custom reconstruction files would
682    % like to do the reconstruction...
683    for ii=1:length(recon_func)
684      try
685        list=recon_func{ii}();
686      catch
687        warning('Error in custom reconsruction function "%s", skipping...',func2str(recon_func{ii}));
688        continue
689      end
690      if any(strcmp(seqfil,list))
691        recon_func_ind = ii;
692        break
693      end
694    end
695  end
696 
697  if recon_func_ind==0
698    Dat.UseCustomRecon = false;
699  end
700 
701  if Dat.UseCustomRecon
702    if Dat.ShowWaitbar
703      wbh=aedes_calc_wait(sprintf('%s\n%s',...
704        ['Using custom function ',upper(func2str(recon_func{recon_func_ind}))],...
705        ['to reconstruct sequence ',procpar.seqfil{1}]));
706      drawnow
707    end
708    [kspace,data,msg_out]=recon_func{recon_func_ind}(kspace,Dat,procpar);
709    if Dat.ShowWaitbar
710      close(wbh)
711    end
712    if isempty(data) && Dat.ReturnFTData
713      % Fourier transform data if not done in custom reconstruction code
714      [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar);
715    end
716  else
717    % Use the default reconstruction code
718    [kspace,msg_out,procpar]=l_ReconstructKspace(kspace,Dat,procpar);
719    DATA.PROCPAR = procpar;
720               
721    % Fourier transform data
722    if Dat.ReturnFTData
723      [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar);
724    end
725  end
726end
727if Dat.ReturnKSpace
728  DATA.KSPACE = kspace;
729else
730  DATA.KSPACE = [];
731end
732if Dat.ReturnFTData
733  DATA.FTDATA = data;
734else
735  DATA.FTDATA = [];
736end
737clear kspace data
738
739
740return
741
742%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
743% Read Data File (Main) Header
744%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
745function [hdr,msg_out]=l_ReadDataFileHeader(file_fid)
746
747msg_out='';
748%% Read Data File Header
749try
750  FH.nblocks = fread(file_fid,1,'long');   % Number of blocks in file
751  FH.ntraces = fread(file_fid,1,'long');   % Number of traces per block
752  FH.np = fread(file_fid,1,'long');        % Number of elements per trace
753  FH.ebytes = fread(file_fid,1,'long');    % Number of bytes per element
754  FH.tbytes = fread(file_fid,1,'long');    % Number of bytes per trace
755  FH.bbytes = fread(file_fid,1,'long');    % Number of bytes per block
756  FH.vers_id = fread(file_fid,1,'short');  % Software version, file_id status bits
757  FH.status = fread(file_fid,1,'short');   % Status of whole file
758  FH.nbheaders = fread(file_fid,1,'long'); % Number of block headers per block
759 
760  hdr.FileHeader = FH;
761 
762  %% Parse status bits
763  hdr.FileHeader.status=[];
764  tmp_str = {'S_DATA',...          % 0 = no data, 1 = data
765             'S_SPEC',...          % 0 = FID, 1 = spectrum
766             'S_32',...            % 0 = 16 bit, 1 = 32 bit integer
767             'S_FLOAT',...         % 0 = integer, 1 = floating point
768             'S_COMPLEX',...       % 0 = real, 1 = complex
769             'S_HYPERCOMPLEX',...  % 1 = hypercomplex
770             '',...                % Unused bit
771             'S_ACQPAR',...        % 0 = not Acqpar, 1 = Acqpar
772             'S_SECND',...         % 0 = first FT, 1 = second FT
773             'S_TRANSF',...        % 0 = regular, 1 = transposed
774             '',...                % Unused bit
775             'S_NP',...            % 1 = np dimension is active
776             'S_NF',...            % 1 = nf dimension is active
777             'S_NI',...            % 1 = ni dimension is active
778             'S_NI2',...           % 1 = ni2 dimension is active
779             ''...                 % Unused bit
780            };
781  status_bits = fliplr(double(dec2bin(FH.status,16))==49);
782  for ii=1:length(tmp_str)
783    if ~isempty(tmp_str{ii})
784      hdr.FileHeader.status.(tmp_str{ii}) = status_bits(ii);
785    end
786  end
787catch
788  hdr=[];
789  msg_out=['Error occurred while reading file header from "' filename '".'];
790  return
791end
792
793
794
795%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
796% Read Block Headers and Data
797%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
798function [hdr,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat)
799
800hdr.BlockHeaders = [];
801%hdr.HypercmplxBHeaders = [];
802kspace=[];
803count = [];
804msg_out='';
805
806
807BlockHeadStatusLabels = {...
808  'S_DATA',...       % 0 = no data, 1 = data
809  'S_SPEC',...          % 0 = FID, 1 = spectrum
810  'S_32',...            % 0 = 16 bit, 1 = 32 bit integer
811  'S_FLOAT',...         % 0 = integer, 1 = floating point
812  'S_COMPLEX',...       % 0 = real, 1 = complex
813  'S_HYPERCOMPLEX',...  % 1 = hypercomplex
814  '',...                % Unused bit
815  'MORE_BLOCKS',...     % 0 = absent, 1 = present
816  'NP_CMPLX',...        % 0 = real, 1 = complex
817  'NF_CMPLX',...        % 0 = real, 1 = complex
818  'NI_CMPLX',...        % 0 = real, 1 = complex
819  'NI2_CMPLX',...       % 0 = real, 1 = complex
820  '',...                % Unused bit
821  '',...                % Unused bit
822  '',...                % Unuesd bit
823  ''...                 % Unused bit
824  };
825
826BlockHeadModeLabels = {...
827  'NP_PHMODE',...   % 1 = ph mode
828  'NP_AVMODE',...    % 1 = av mode
829  'NP_PWRMODE',...   % 1 = pwr mode
830  '',...             % Unused bit
831  'NF_PHMODE',...    % 1 = ph mode
832  'NF_AVMODE',...    % 1 = av mode
833  'NF_PWRMODE',...   % 1 = pwr mode
834  '',...             % Unused bit
835  'NI_PHMODE',...    % 1 = ph mode
836  'NI_AVMODE',...    % 1 = av mode
837  'NI_PWRMODE',...   % 1 = pwr mode
838  '',...             % Unused bit
839  'NI2_PHMODE',...   % 1 = ph mode
840  'NI2_AVMODE',...   % 1 = av mode
841  'NI2_PWRMODE',...  % 1 = pwr mode
842  ''...              % Unused bit
843  };
844
845
846% The nbheaders -field is not correct in some cases. Thus, this field
847% cannot be trusted and the real nbheaders has to be calculated from the
848% byte counts.
849tmp_bytes=hdr.FileHeader.bbytes-hdr.FileHeader.tbytes*hdr.FileHeader.ntraces;
850header_bytes=28;
851if rem(tmp_bytes,header_bytes)~=0
852  msg_out = 'Block header byte count does not match with file header';
853  return
854else
855  nbheaders = tmp_bytes/28;
856end
857
858% Allocate space for kspace
859% kspace = complex(zeros(hdr.FileHeader.np/2,...
860%           hdr.FileHeader.ntraces*hdr.FileHeader.nblocks,...
861%           Dat.precision));
862
863kspace = complex(zeros(hdr.FileHeader.np/2,...
864        hdr.FileHeader.ntraces*hdr.FileHeader.nblocks,Dat.precision));
865
866
867%% - The older robust (but also slower) way for reading the data.
868%% When the blocksize is relatively small, this is also quite fast.
869if ~Dat.FastDataRead
870 
871  % Initialize waitbar
872  if Dat.ShowWaitbar
873    wb_h = aedes_wbar(0/hdr.FileHeader.nblocks,...
874      {['Reading VNMR data.'],...
875      ['(Processing data block ' ...
876      num2str(0) '/' num2str(hdr.FileHeader.nblocks) ')']});
877  end
878
879  %% Read data and block headers
880  for ii=1:hdr.FileHeader.nblocks
881    %% Update waitbar
882    if Dat.ShowWaitbar
883      aedes_wbar(ii/hdr.FileHeader.nblocks,...
884        wb_h,...
885        {['Reading VNMR data.'],...
886        ['(Processing data block ' ...
887        num2str(ii) '/' num2str(hdr.FileHeader.nblocks) ')']})
888    end
889
890    %% Read block header and hypercomplex header
891    for kk=1:nbheaders
892      %% Read block header
893      if kk==1
894        hdr.BlockHeaders.scale = fread(file_fid,1,'short');  % Scaling factor
895        tmp_status = fread(file_fid,1,'short'); % Status of data in block
896        hdr.BlockHeaders.status = [];
897        hdr.BlockHeaders.index = fread(file_fid,1,'short');  % Block index
898        tmp_mode = fread(file_fid,1,'short');   % Mode of data in block
899        hdr.BlockHeaders.mode = [];
900        hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
901        hdr.BlockHeaders.lpval = fread(file_fid,1,'float');  % f2 (2D-f1) left phase in phasefile
902        hdr.BlockHeaders.rpval = fread(file_fid,1,'float');  % f2 (2D-f1) right phase in phasefile
903        hdr.BlockHeaders.lvl = fread(file_fid,1,'float');    % level drift correction
904        hdr.BlockHeaders.tlt = fread(file_fid,1,'float');    % tilt drift correction
905
906        %% Parse status and mode bits
907        status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
908        mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
909        for tt=1:length(BlockHeadStatusLabels)
910          if ~isempty(BlockHeadStatusLabels{tt})
911            hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
912          end
913          if ~isempty(BlockHeadModeLabels{tt})
914            hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
915          end
916        end
917
918
919      else %% Read hypercomplex header
920        fread(file_fid,1,'short'); % Spare
921        hdr.BlockHeaders.HCHstatus = fread(file_fid,1,'short');
922        fread(file_fid,1,'short'); % Spare
923        fread(file_fid,1,'short'); % Spare
924        fread(file_fid,1,'long'); % Spare
925        hdr.BlockHeaders.HCHlpval1 = fread(file_fid,1,'float');
926        hdr.BlockHeaders.HCHrpval1 = fread(file_fid,1,'float');
927        fread(file_fid,1,'float'); % Spare
928        fread(file_fid,1,'float'); % Spare
929      end
930    end
931
932    %% Check block index to be sure about the data type
933    if hdr.BlockHeaders.index~=ii
934      fprintf(1,'Warning: Index mismatch in "%s" block %d\n',fopen(file_fid),ii);
935
936      % Use information from the file header instead of the BlockHeader if
937      % there is a mismatch in blockheader index...
938      useFileHeader = true;
939    else
940      useFileHeader = false;
941    end
942
943    %% Determine data precision
944    if useFileHeader
945      if hdr.FileHeader.status.S_FLOAT==1
946        prec_str = ['single=>',Dat.precision];
947      elseif hdr.FileHeader.status.S_32==1 ...
948          && hdr.FileHeader.status.S_FLOAT==0
949        prec_str = ['int32=>',Dat.precision];
950      elseif hdr.FileHeader.status.S_32==0 ...
951          && hdr.FileHeader.status.S_FLOAT==0
952        prec_str = ['int16=>',Dat.precision];
953      end
954
955    else
956      if hdr.BlockHeaders.status.S_FLOAT==1
957        prec_str = ['single=>',Dat.precision];
958      elseif hdr.BlockHeaders.status.S_32==1 ...
959          && hdr.BlockHeaders.status.S_FLOAT==0
960        prec_str = ['int32=>',Dat.precision];
961      elseif hdr.BlockHeaders.status.S_32==0 ...
962          && hdr.BlockHeaders.status.S_FLOAT==0
963        prec_str = ['int16=>',Dat.precision];
964      end
965    end
966
967    % Read k-space
968    tmp=fread(file_fid,...
969      [hdr.FileHeader.np,hdr.FileHeader.ntraces],...
970      prec_str);
971   
972
973    % Store complex block and perform DC correction
974    if ~Dat.DCcorrection
975      complex_block = complex(tmp(1:2:end,:),tmp(2:2:end,:));
976    else
977      complex_block = complex(tmp(1:2:end,:)-hdr.BlockHeaders.lvl,...
978        tmp(2:2:end,:)-hdr.BlockHeaders.tlt);
979    end
980   
981    inds = ((ii-1)*hdr.FileHeader.ntraces+1):(ii*hdr.FileHeader.ntraces);
982    kspace(:,inds) = complex_block;
983
984    % Do not save blockheaders by default. When reading data files with a lot of
985    % blocks (e.g. over 1000) the processing of the DATA structure can be
986    % slowed down considerably. If you for some reason want to save also the
987    % block headers in the DATA structure comment out the code line below.
988    hdr.BlockHeaders = [];
989  end % for ii=1:hdr.
990
991  if Dat.ShowWaitbar
992    delete(wb_h)
993  end
994
995else
996  %% -------------------------------------------------------------------
997  %% Fast Method for reading data. This may fail with some datas and can
998  %% also require a relatively large amount of memory. This
999  %% method should be used for EPI datas that contain a large number
1000  %% of block headers...
1001
1002  % Check the size of the FID-file
1003  d=dir(fopen(file_fid));
1004  file_sz = d.bytes/1024/1024; % File size in MB
1005  if file_sz<Dat.FileChunkSize
1006    nBlocks = 1;
1007  else
1008    nBlocks = ceil(file_sz/Dat.FileChunkSize); % Read data in 500 MB blocks
1009  end
1010
1011  % Initialize waitbar
1012  if Dat.ShowWaitbar
1013    if nBlocks==1
1014      wb_h = aedes_calc_wait(['Reading VNMR data.']);
1015    else
1016      wb_h = aedes_wbar(1/nBlocks,{['Reading VNMR data.'],...
1017        sprintf('Reading block 0/%d',nBlocks)});
1018    end
1019  end
1020
1021  % The first block header is read and it is assumed that the values in
1022  % the other block headers don't change.
1023  hdr.BlockHeaders.scale = fread(file_fid,1,'short');  % Scaling factor
1024  tmp_status = fread(file_fid,1,'short'); % Status of data in block
1025  hdr.BlockHeaders.status = [];
1026  hdr.BlockHeaders.index = fread(file_fid,1,'short');  % Block index
1027  tmp_mode = fread(file_fid,1,'short');   % Mode of data in block
1028  hdr.BlockHeaders.mode = [];
1029  hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
1030  hdr.BlockHeaders.lpval = fread(file_fid,1,'float');  % f2 (2D-f1) left phase in phasefile
1031  hdr.BlockHeaders.rpval = fread(file_fid,1,'float');  % f2 (2D-f1) right phase in phasefile
1032  hdr.BlockHeaders.lvl = fread(file_fid,1,'float');    % level drift correction
1033  hdr.BlockHeaders.tlt = fread(file_fid,1,'float');    % tilt drift correction
1034
1035  %% Parse status and mode bits
1036  status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
1037  mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
1038  for tt=1:length(BlockHeadStatusLabels)
1039    if ~isempty(BlockHeadStatusLabels{tt})
1040      hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
1041    end
1042    if ~isempty(BlockHeadModeLabels{tt})
1043      hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
1044    end
1045  end
1046
1047  %% Determine data precision
1048  if hdr.BlockHeaders.status.S_FLOAT==1
1049    prec_str = ['single=>',Dat.precision];
1050    prec = 4; % Precision in bytes
1051  elseif hdr.BlockHeaders.status.S_32==1 ...
1052      && hdr.BlockHeaders.status.S_FLOAT==0
1053    prec_str = ['int32=>',Dat.precision];
1054    prec = 4;
1055  elseif hdr.BlockHeaders.status.S_32==0 ...
1056      && hdr.BlockHeaders.status.S_FLOAT==0
1057    prec_str = ['int16=>',Dat.precision];
1058    prec = 2;
1059  end
1060
1061  % Seek the file back to the beginning of the first block header
1062  fseek(file_fid,-28,0);
1063
1064  % Determine the number of values that will result from block header(s)
1065  nVals = (nbheaders*28)/prec;
1066
1067  nbh = floor(hdr.FileHeader.nblocks/nBlocks);
1068  szh = nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces;
1069  for ii=1:nBlocks
1070    if nBlocks~=1
1071      aedes_wbar(ii/nBlocks,wb_h,{['Reading VNMR data.'],...
1072        sprintf('Reading block %d/%d',ii,nBlocks)});
1073    end
1074
1075    % Read the whole data including block headers etc...
1076    if ii==nBlocks
1077      tmp = fread(file_fid,inf,prec_str);
1078    else
1079      tmp = fread(file_fid,nbh*szh,prec_str);
1080    end
1081    tmp=reshape(tmp,nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces,[]);
1082    tmp(1:nVals,:)=[];
1083    tmp=reshape(tmp,hdr.FileHeader.np,[]);
1084
1085    if ii==nBlocks
1086      inds = ((ii-1)*nbh*hdr.FileHeader.ntraces+1):size(kspace,2);
1087    else
1088      inds = ((ii-1)*nbh*hdr.FileHeader.ntraces+1):ii*nbh*hdr.FileHeader.ntraces;
1089    end
1090
1091    % Do DC-correction if necessary
1092    if ~Dat.DCcorrection
1093      kspace(:,inds)=complex(tmp(1:2:end,:,:),tmp(2:2:end,:,:));
1094    else
1095      kspace(:,inds)=complex(tmp(1:2:end,:,:)-hdr.BlockHeaders.lvl,...
1096        tmp(2:2:end,:,:)-hdr.BlockHeaders.tlt);
1097    end
1098  end
1099  clear tmp
1100
1101  if Dat.ShowWaitbar
1102    delete(wb_h)
1103    pause(0.1)
1104  end
1105
1106end
1107
1108% ==================================================================
1109% Reconstruct k-space
1110% ==================================================================
1111function [kspace,msg_out,procpar]=l_ReconstructKspace(kspace,Dat,procpar)
1112
1113msg_out = '';
1114
1115% Get the strategic parameters
1116seqcon = procpar.seqcon{1};
1117seqfil = procpar.seqfil{1};
1118np = procpar.np;
1119nv = procpar.nv;
1120nv2 = procpar.nv2;
1121nv3 = procpar.nv3;
1122ns = procpar.ns;
1123if isfield(procpar,'etl')
1124  etl = procpar.etl;
1125else
1126  etl = 1;
1127end
1128pss = procpar.pss;
1129nt = procpar.nt;
1130nf = procpar.nf;
1131ne = procpar.ne;
1132if isfield(procpar,'flash_converted')
1133  % Don't try to sort already sorted data...
1134  Dat.Sorting = false;
1135  seqcon(2:3) = 'sc';
1136end
1137if isfield(procpar,'nseg')
1138        nseg = procpar.nseg;
1139else
1140        nseg = 1;
1141end
1142
1143
1144% Check dimensions
1145if nv==0
1146  % 1D-image (specrum)
1147  AcqType = 1;
1148elseif nv2==0
1149  % 2D-image
1150  AcqType = 2;
1151elseif nv3==0
1152  % 3D-image
1153  AcqType = 3;
1154else
1155  AcqType = 4;
1156end
1157
1158% Check number of receivers
1159if isfield(procpar,'rcvrs') && length(procpar.rcvrs{1})>1
1160  nRcvrs = length(find(procpar.rcvrs{1}=='y'));
1161else
1162  nRcvrs = 1;
1163end
1164
1165
1166% Check for arrayed acquisition
1167if isfield(procpar,'array') && ~isempty(procpar.array{1})
1168  %if length(procpar.array)==1 && ~iscell(procpar.array{1}) && ...
1169  %    strcmp(procpar.array{1},'pad') && all(procpar.pad==0)
1170  %  % Skip the array parsing if the array is a dummy "pad"-array...
1171  %  isAcqArrayed = false;
1172  %  ArrayLength = 1;
1173  %else
1174    isAcqArrayed = true;
1175    ArrayLength = [];
1176   
1177    % Determine array length
1178    for ii=1:length(procpar.array)
1179      if iscell(procpar.array{ii})
1180        ArrayLength(ii) = length(procpar.(procpar.array{ii}{1}));
1181      else
1182        ArrayLength(ii) = length(procpar.(procpar.array{ii}));
1183      end
1184    end
1185    ArrayLength = prod(ArrayLength);
1186  %end
1187else
1188  isAcqArrayed = false;
1189  ArrayLength = 1;
1190end
1191
1192
1193% Reconstruct k-space ----------------------
1194if AcqType==1
1195  % Reconstruct 1D data ...
1196elseif AcqType==2
1197 
1198  % Reconstruct 2D data
1199  if strcmpi(seqcon,'ncsnn')
1200    kspace = reshape(kspace,[np/2,etl,ns,nRcvrs,ArrayLength,nv/etl]);
1201    kspace = permute(kspace,[1 2 6 3 5 4]);
1202    kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1203  elseif strcmpi(seqcon,'nscnn')
1204    if isfield(procpar,'flash_converted')
1205      kspace = reshape(kspace,[np/2,nRcvrs,ArrayLength,ns,nv]);
1206      kspace = permute(kspace,[1 5 4 3 2]);
1207      kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1208    else
1209      kspace = reshape(kspace,[np/2,nRcvrs,nv,ArrayLength,ns]);
1210      kspace = permute(kspace,[1 3 5 4 2]);
1211      kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1212    end
1213  elseif strcmpi(seqcon,'nccnn')
1214    kspace = reshape(kspace,[np/2,etl,ns,nv/etl,ArrayLength,nRcvrs]);
1215    kspace = permute(kspace,[1 2 4 3 5 6]);
1216    kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1217  end
1218elseif AcqType==3
1219  % Reconstruct 3D data
1220  if strcmpi(seqcon,'nscsn')
1221    kspace = reshape(kspace,[np/2,etl,nv/etl,nRcvrs,ArrayLength*nv2]);
1222    kspace = permute(kspace,[1 2 3 5 4]);
1223    kspace = reshape(kspace,[np/2,nv,nv2,ArrayLength,nRcvrs]);
1224        end
1225        if strcmpi(seqcon,'ccccn')
1226                kspace = reshape(kspace,[np/2,etl,nRcvrs,ne,nv/etl,ArrayLength*nv2]);
1227                kspace = permute(kspace,[1 2 5 6 4 3]);
1228                kspace = reshape(kspace,[np/2,nv,nv2,ne*ArrayLength,nRcvrs]);
1229        end
1230else
1231  % Reconstruct nD data ... to be written
1232end
1233
1234
1235% Sort data using phasetable ------------------
1236if Dat.Sorting && ~isempty(Dat.phasetable)
1237  Dat.phasetable = Dat.phasetable.';
1238  kspace(:,Dat.phasetable(:),:,:,:)=kspace;
1239end
1240
1241% Sort interleaved 2D data using pss
1242[sorted_pss,I_pss]=sort(procpar.pss);
1243if ~ismember(procpar.pss,sorted_pss,'rows') && Dat.SortPSS
1244        kspace = kspace(:,:,I_pss,:,:);
1245        procpar.pss_orig = procpar.pss;
1246        procpar.pss = sorted_pss;
1247end
1248
1249% ==================================================================
1250% Fourier transform k-space
1251% ==================================================================
1252function [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar)
1253
1254msg_out = '';
1255data = [];
1256
1257
1258% Check dimensions
1259if procpar.nv==0
1260  % 1D-image
1261  AcqType = 1;
1262elseif procpar.nv2==0
1263  % 2D-image
1264  AcqType = 2;
1265elseif procpar.nv3==0
1266  % 3D-image
1267  AcqType = 3;
1268else
1269  AcqType = 4;
1270end
1271
1272% Check number of receivers
1273if isfield(procpar,'rcvrs') && length(procpar.rcvrs{1})>1
1274  nRcvrs = length(find(procpar.rcvrs{1}=='y'));
1275else
1276  nRcvrs = 1;
1277end
1278
1279% If zeropadding is requested, calculate the padded size
1280if Dat.ZeroPadding~=0
1281  if Dat.ZeroPadding==1
1282    % Zeropad to square
1283    if AcqType==1
1284      padSize = procpar.np/2;
1285    elseif AcqType==2
1286      padSize = ones(1,2)*procpar.np/2;
1287      padSize(3) = size(kspace,3);
1288    else
1289      padSize = ones(1,3)*procpar.np/2;
1290    end
1291  else
1292    % Zeropadding is on "auto", i.e. zeropad to FOV
1293    lpe = procpar.lpe;
1294    lpe2 = procpar.lpe2;
1295    lro = procpar.lro;
1296    if AcqType==2
1297      % 2D data
1298      padSize = [procpar.np/2 ...
1299        procpar.np/2*(lpe/lro) ...
1300        size(kspace,3)];
1301    elseif AcqType==3 && lpe2~=0
1302      % 3D data
1303      padSize = [procpar.np/2 ...
1304        procpar.np/2*(lpe/lro) ...
1305        procpar.np/2*(lpe2/lro)];
1306    end
1307  end
1308  ks_sz = [size(kspace,1) ...
1309    size(kspace,2) ...
1310    size(kspace,3)];
1311        padSize = round(padSize);
1312  if any(padSize>ks_sz)
1313    %kspace(padSize) = 0;
1314                kspace(padSize(1),padSize(2),padSize(3)) = 0;
1315  end
1316else
1317  padSize = [size(kspace,1) ...
1318    size(kspace,2) ...
1319    size(kspace,3)];
1320end
1321
1322% Allocate space for Fourier transformed data
1323if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1324  data_sz = [padSize,size(kspace,4),size(kspace,5)+1];
1325  data = zeros(data_sz,Dat.precision);
1326else
1327  data = zeros(size(kspace),Dat.precision);
1328end
1329%data = [];
1330%if strcmpi(Dat.precision,'single')
1331%  data = single(data);
1332%end
1333
1334% Fourier transform data
1335if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1336  ind = [2:size(data,5)];
1337else
1338  ind = [1:size(data,5)];
1339end
1340if AcqType==1
1341  data(:,:,:,:,ind) = abs(fftshift(ifft(kspace,[],1),1));
1342elseif AcqType==2
1343  data(:,:,:,:,ind) = abs(fftshift(fftshift(ifft(ifft(kspace,[],1),[],2),1),2));
1344elseif AcqType==3
1345  data(:,:,:,:,ind) = abs(fftshift(fftshift(fftshift(ifft(ifft(ifft(kspace,[],1),[],2),[],3),1),2),3));
1346end
1347
1348% Calculate sum-of-squares image
1349if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1350  % Calculate sum-of-squares
1351  data(:,:,:,:,1) = sqrt(sum(data(:,:,:,:,ind).^2,5));
1352  data=abs(data);
1353  if Dat.SumOfSquares==1
1354    % Remove individual receiver data
1355    data=data(:,:,:,:,1);
1356  end
1357end
1358
1359% ==================================================================
1360% Find custom functions for VNMR k-space reconstruction
1361% ==================================================================
1362function recon_func=l_GetReconFunc(recon_dir)
1363
1364recon_func = {};
1365
1366if ~isdir(recon_dir)
1367  return
1368end
1369
1370dir_struct=dir(recon_dir);
1371recon_files = {dir_struct(~[dir_struct(:).isdir]).name};
1372if isempty(recon_files)
1373  return
1374end
1375
1376% Remove files that don't have .m extension
1377ind = regexpi(recon_files,'\.m$');
1378recon_files = {recon_files{~cellfun('isempty',ind)}};
1379
1380currentDir = pwd;
1381try
1382  cd(recon_dir);
1383  for ii=1:length(recon_files)
1384    recon_func{ii}=str2func(recon_files{ii}(1:end-2));
1385  end
1386  cd(currentDir);
1387catch
1388  cd(currentDir);
1389end
1390
1391
1392
Note: See TracBrowser for help on using the repository browser.

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