source: aedes_readvnmr.m

Last change on this file was 209, checked in by tjniskan, 4 years ago
  • Fixed a bug small in aedes_readvnmr
  • Added support in readvnmr for mge3d data

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

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

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