source: aedes_readvnmr.m @ 148

Last change on this file since 148 was 148, checked in by tjniskan, 9 years ago
  • Added aedes_roifill.m for doing binary flood fill operation. This is

much slower than using the imfill-function, but does not dependend on
Image Processing toolbox. The goal is to make Aedes run with just the
base Matlab install. DICOM support will probably be the last obstacle
in achieving this, because it's going to be a great pain in the a to
code...

  • Fixed overwriting ROIs in aedes_roi_copy_gui.m and added additional

information in the ROI copy button tooltips.

  • Some further iterations to the new readvnmr-function.

M aedes_roi_copy_gui.m
M aedes_getmatlabversion.m
M aedes_readvnmr.m
M aedes.m
M aedes_helpabout.m
M aedes_readfid.m
M aedes_revision.m
A aedes_roifill.m

File size: 42.0 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 
1152  % Reconstruct 2D data
1153  if strcmpi(seqcon,'ncsnn')
1154    kspace = reshape(kspace,[np/2,etl,ns,ArrayLength,nv/etl,nRcvrs]);
1155    kspace = permute(kspace,[1 2 5 3 4 6]);
1156    kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1157  elseif strcmpi(seqcon,'nscnn')
1158    if isfield(procpar,'flash_converted')
1159      kspace = reshape(kspace,[np/2,nRcvrs,ArrayLength,nv,ns]);
1160      kspace = permute(kspace,[1 4 5 3 2]);
1161      kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1162    else
1163      kspace = reshape(kspace,[np/2,nRcvrs,nv,ArrayLength,ns]);
1164      kspace = permute(kspace,[1 3 5 4 2]);
1165      kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1166    end
1167  elseif strcmpi(seqcon,'nccnn')
1168    kspace = reshape(kspace,[np/2,etl,nRcvrs,ns,nv/etl,ArrayLength]);
1169    kspace = permute(kspace,[1 2 5 4 6 3]);
1170    kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]);
1171  end
1172elseif AcqType==3
1173  % Reconstruct 3D data
1174  if strcmpi(seqcon,'nscsn')
1175    kspace = reshape(kspace,[np/2,etl,nv/etl,nRcvrs,ArrayLength*nv2]);
1176    kspace = permute(kspace,[1 2 3 5 4]);
1177    kspace = reshape(kspace,[np/2,nv,nv2,ArrayLength,nRcvrs]);
1178  end
1179else
1180  % Reconstruct nD data ...
1181end
1182
1183
1184% Sort data using phasetable ------------------
1185if Dat.Sorting && ~isempty(Dat.phasetable)
1186  Dat.phasetable = Dat.phasetable.';
1187  kspace(:,Dat.phasetable(:),:,:,:)=kspace;
1188end
1189
1190
1191
1192% ==================================================================
1193% Fourier transform k-space
1194% ==================================================================
1195function [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar)
1196
1197msg_out = '';
1198data = [];
1199
1200
1201% Check dimensions
1202if procpar.nv==0
1203  % 1D-image
1204  AcqType = 1;
1205elseif procpar.nv2==0
1206  % 2D-image
1207  AcqType = 2;
1208elseif procpar.nv3==0
1209  % 3D-image
1210  AcqType = 3;
1211else
1212  AcqType = 4;
1213end
1214
1215% Check number of receivers
1216if isfield(procpar,'rcvrs') && length(procpar.rcvrs{1})>1
1217  nRcvrs = length(find(procpar.rcvrs{1}=='y'));
1218else
1219  nRcvrs = 1;
1220end
1221
1222% If zeropadding is requested, calculate the padded size
1223if Dat.ZeroPadding~=0
1224  if Dat.ZeroPadding==1
1225    % Zeropad to square
1226    if AcqType==1
1227      padSize = procpar.np/2;
1228    elseif AcqType==2
1229      padSize = ones(1,2)*procpar.np/2;
1230      padSize(3) = size(kspace,3);
1231    else
1232      padSize = ones(1,3)*procpar.np/2;
1233    end
1234  else
1235    % Zeropadding is on "auto", i.e. zeropad to FOV
1236    lpe = procpar.lpe;
1237    lpe2 = procpar.lpe2;
1238    lro = procpar.lro;
1239    if AcqType==2
1240      % 2D data
1241      padSize = [procpar.np/2 ...
1242        procpar.np/2*(lpe/lro) ...
1243        size(kspace,3)];
1244    elseif AcqType==3 && lpe2~=0
1245      % 3D data
1246      padSize = [procpar.np/2 ...
1247        procpar.np/2*(lpe/lro) ...
1248        procpar.np/2*(lpe2/lro)];
1249    end
1250  end
1251  ks_sz = [size(kspace,1) ...
1252    size(kspace,2) ...
1253    size(kspace,3)];
1254  if any(padSize>ks_sz)
1255    kspace(padSize) = 0;
1256  end
1257else
1258  padSize = [size(kspace,1) ...
1259    size(kspace,2) ...
1260    size(kspace,3)];
1261end
1262
1263% Allocate space for Fourier transformed data
1264if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1265  data_sz = [padSize,size(kspace,4),size(kspace,5)+1];
1266  data = zeros(data_sz,Dat.precision);
1267else
1268  data = zeros(size(kspace),Dat.precision);
1269end
1270%data = [];
1271%if strcmpi(Dat.precision,'single')
1272%  data = single(data);
1273%end
1274
1275% Fourier transform data
1276if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1277  ind = [2:size(data,5)];
1278else
1279  ind = [1:size(data,5)];
1280end
1281if AcqType==1
1282  data(:,:,:,:,ind) = abs(fftshift(ifft(kspace,[],1),1));
1283elseif AcqType==2
1284  data(:,:,:,:,ind) = abs(fftshift(fftshift(ifft(ifft(kspace,[],1),[],2),1),2));
1285elseif AcqType==3
1286  data(:,:,:,:,ind) = abs(fftshift(fftshift(fftshift(ifft(ifft(ifft(kspace,[],1),[],2),[],3),1),2),3));
1287end
1288
1289% Calculate sum-of-squares image
1290if nRcvrs>1 && any(Dat.SumOfSquares==[1 2])
1291  % Calculate sum-of-squares
1292  data(:,:,:,:,1) = sqrt(sum(data(:,:,:,:,ind).^2,5));
1293  data=abs(data);
1294  if Dat.SumOfSquares==1
1295    % Remove individual receiver data
1296    data=data(:,:,:,:,1);
1297  end
1298end
1299
1300% ==================================================================
1301% Find custom functions for VNMR k-space reconstruction
1302% ==================================================================
1303function recon_func=l_GetReconFunc(recon_dir)
1304
1305recon_func = {};
1306
1307if ~isdir(recon_dir)
1308  return
1309end
1310
1311dir_struct=dir(recon_dir);
1312recon_files = {dir_struct(~[dir_struct(:).isdir]).name};
1313if isempty(recon_files)
1314  return
1315end
1316
1317% Remove files that don't have .m extension
1318ind = regexpi(recon_files,'\.m$');
1319recon_files = {recon_files{~cellfun('isempty',ind)}};
1320
1321currentDir = pwd;
1322try
1323  cd(recon_dir);
1324  for ii=1:length(recon_files)
1325    recon_func{ii}=str2func(recon_files{ii}(1:end-2));
1326  end
1327  cd(currentDir);
1328catch
1329  cd(currentDir);
1330end
1331
1332
1333
Note: See TracBrowser for help on using the repository browser.

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