source: aedes_readvnmr.m @ 138

Last change on this file since 138 was 138, checked in by tjniskan, 9 years ago
  • Added option to VNMR preferences to use the new version of the VNMR

read function

  • Added read support for certain CINE sequences
  • Added a command line preference option for inputting data aspect

ratio manually in created mosaic images

M aedes_readfidprefs.m
M aedes_data_read.m
M aedes_readvnmr.m
M aedes_export_gui.m
M aedes_revision.m
A vnmr_recon/cine_recon.m

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

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