source: aedes_readvnmr.m @ 133

Last change on this file since 133 was 128, checked in by tjniskan, 9 years ago
  • Added first alpha version of the new VNMR-data reader

aedes_readvnmr.m (the old aedes_readfid is becoming more and more
cluttered with endless patches here and there)

  • Added the directory for custom codes for reconstruction VNMR data.
  • Added first alpha version of custom code for reconstruction EPI data

from INOVA system

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

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

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