source: aedes_readvnmr.m @ 162

Last change on this file since 162 was 161, checked in by tjniskan, 8 years ago
  • Added support for a certain multiecho 3D sequences in aedes_readvnmr.m
  • Fixed a small bug in histograms plugin

M aedes_readvnmr.m
M plugins/histograms.m
M aedes_revision.m

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

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