source: aedes_readfid.m @ 119

Last change on this file since 119 was 119, checked in by tjniskan, 9 years ago
  • Added support for triple-reference EPI data.
  • Added two new properties for aedes_readfid:
    • FileChunkSize? for reducing memory requirements when reading large fid-files (default FileChunkSize?=500 MB)
    • EPIBlockSize for reducing memory requirements when reading large multireceiver EPI (default=100 volumes)
  • Tried to speed up smoothing in fMRI analysis by performing the

operation in frequency space. Still needs some work...

M misclib/fmri_smooth.m
M misclib/fmri_analysis.m
M aedes_readfid.m
M aedes_revision.m

File size: 59.5 KB
Line 
1function [DATA,msg_out] = aedes_readfid(filename,varargin)
2% AEDES_READFID - Read VNMR (Varian) FID-files
3%   
4%
5% Synopsis:
6%        [DATA,msg]=aedes_readfid(filename,'PropertyName1',value1,'PropertyName2',value2,...)
7%        [DATA,msg]=aedes_readfid(filename,'header')
8%        [DATA,msg]=aedes_readfid(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_READFID 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_READFID (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%        'FileChunkSize' : [ megabytes ]        % Chunk size in megabytes for
62%                                               % reading fid-files
63%                                               % (default=500 MB).
64%                                               % Lowering the chunk size
65%                                               % helps to conserve memory
66%                                               % when reading large files,
67%                                               % but is slower. This
68%                                               % property has no effect if
69%                                               % FastRead='off'.
70%                                               
71%
72%        'OutputFile'  :  filename              % Writes the slices into
73%                                               % NIfTI-format files
74%                                               % using the given filename.
75%
76%        'DCcorrection':  [ {'off'} | 'on' ]    % 'on'=use DC correction
77%                                               % 'off'=don't use DC correction
78%                                               % (default)
79%
80%        'Procpar'     :  (procpar-structure)   % Input procpar
81%                                               % structure. If this
82%                                               % property is not set the
83%                                               % procpar structure
84%                                               % will be read using
85%                                               % AEDES_READPROCPAR
86%
87%        'ZeroPadding' :  ['off'|'on'|{'auto'}] % 'off' = force
88%                                               % zeropadding off
89%                                               % 'on' = force
90%                                               % zeropadding on (force
91%                                               % images to be square)
92%                                               % 'auto' = autoselect
93%                                               % using relative FOV
94%                                               % dimensions (default)
95%
96%        'PhaseTable'  : (custom phase table)   % User-specified
97%                                               % line order for
98%                                               % k-space. The phase table
99%                                               % is obtained from the file
100%                                               % specified in
101%                                               % procpar.petable if
102%                                               % necessary.
103%
104%        'sorting'      : [ 'off' | {'on'} ]    % 'off' =disable k-space
105%                                               % sorting, i.e. ignore the
106%                                               % phase table and/or arrays.
107%                                               % 'on' =sort k-space using
108%                                               % phase table and/or array
109%                                               % information if necessary
110%                                               % (default)
111%
112%        'wbar'        : [ {'on'} | 'off' ]     % 'on'  = show waitbar
113%                                               % 'off' = don't show waitbar
114%
115%        'FlipKspace'  : [ {'off'} | 'LR' |
116%                             'UD' | 'LRUD' ]   % 'off' = don't flip (default)
117%                                               % 'LR' = left/right
118%                                               % 'UD' = up/down
119%                                               % 'LRUD' = left/right and
120%                                               % up/down
121%
122%        'FlipInd'     : [ {'all'} | 'alt' |
123%                          (custom vector)  ]   % 'all' = flip all slices
124%                                               % (default)
125%                                               % 'alt' = flip every
126%                                               % second slice
127%                                               % (custom vector) = A
128%                                               % custom vector containing
129%                                               % indices to the flipped slices
130%
131%        'Precision'   : ['single'|{'double'}]  % Precision of the
132%                                               % outputted data.
133%                                               % 'single' = 32-bit float
134%                                               % 'double' = 64-bit float
135%
136%        'OrientImages': [ {'on'} | 'off' ]     % Orient FTDATA after
137%                                               % reading the data using
138%                                               % PROCPAR.orient property.
139%
140%        'RemoveEPIphaseIm' : [{'on'}|'off']    % Remove the phase image
141%                                               % from EPI data. This
142%                                               % option only affects EPI
143%                                               % images.
144%
145%        'EPIBlockSize' : [integer value]       % Block size (number of
146%                                               % volumes) used for
147%                                               % processing multireceiver
148%                                               % EPI data. Default=100
149%
150% Examples:
151%        [DATA,msg]=aedes_readfid(filename)             % Read image data from 'filename'
152%        [DATA,msg]=aedes_readfid(filename,'header')    % Read only data file header
153%
154%        [DATA,msg]=aedes_readfid(filename,'return',1)  % Return only image data (default)
155%        [DATA,msg]=aedes_readfid(filename,'return',2)  % Return only k-space
156%        [DATA,msg]=aedes_readfid(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_readfid(filename,'header')
160%        [DATA,msg]=aedes_readfid(DATA.HDR,'return',3)
161%
162%        % Read VNMR-data with default options and save slices in NIfTI
163%        % format
164%        DATA=aedes_readfid('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_readfid('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) 2006 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
192
193
194if nargout==2
195  Dat.ShowErrors = false;
196  msg_out='';
197else
198  Dat.ShowErrors = true;
199end
200
201%% ---------------------
202% Defaults
203Dat.ReturnKSpace  = false;
204Dat.ReturnFTData  = true;
205Dat.DCcorrection  = false;
206Dat.ZeroPadding = 2;
207Dat.Sorting = true;
208Dat.UsePhaseTable = true;
209Dat.FastDataRead = true;
210Dat.precision = 'single';
211Dat.FileChunkSize = 500; % Chunk size (in MB) for FastRead
212
213%% Other defaults
214Dat.ShowWaitbar   = true;
215procpar=[];
216Dat.phasetable=[];
217Dat.FlipKspace = 0;
218Dat.FlipInd = 'all';
219Dat.OutputFile = false;
220Dat.ReturnRawKspace = false;
221Dat.ReorientEPI = false;
222Dat.RemoveEPIphaseIm = true;
223Dat.EPIBlockSize = 100;
224Dat.OrientImages = true;
225
226%% -------------------------------------------------------------------
227
228
229%% Set data format label
230DATA.DataFormat = 'vnmr';
231
232%% Parse input arguments
233if nargin==0 || isempty(filename)
234 
235  %% Get default directory
236  try
237    defaultdir = getpref('Aedes','GetDataFileDir');
238  catch
239    defaultdir = [pwd,filesep];
240  end
241 
242  %% If no input arguments are given, ask for a file
243  [f_name, f_path, f_index] = uigetfile( ...
244       {'fid;FID','Varian FID-files (fid)'; ...
245        '*.*', 'All Files (*.*)'}, ...
246        'Select VNMR (Varian) FID file',defaultdir);
247  if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed
248    DATA=[];
249    msg_out='Canceled';
250    return
251  end
252  ReadHdr = true;
253  ReadData = true;
254  filename=[f_path,f_name];
255 
256  %% Set default directory to preferences
257  setpref('Aedes','GetDataFileDir',f_path)
258 
259end
260
261if nargin==1
262  if isstruct(filename)
263    hdr = filename;
264    filename = [hdr.fpath,hdr.fname];
265    ReadHdr = false;
266    ReadData = true;
267% $$$     ReturnKSpace = false;
268% $$$     ReturnFTData = true;
269   
270  elseif ischar(filename)
271    ReadHdr = true;
272    ReadData = true;
273% $$$     ReturnKSpace = false;
274% $$$     ReturnFTData = true;
275  end
276elseif nargin==2
277  if strcmpi(varargin{1},'header')
278    ReadHdr = true;
279    ReadData = false;
280% $$$     ReturnKSpace = 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 'phasetable'
372      Dat.phasetable = varargin{ii+1};
373     
374      case 'sorting'
375          if strcmpi(varargin{ii+1},'on')
376        Dat.UsePhaseTable = true;
377        Dat.Sorting = true;
378      else
379        Dat.UsePhaseTable = false;
380        Dat.Sorting = false;
381          end
382     
383         case 'fastread'
384     if strcmpi(varargin{ii+1},'on')
385       Dat.FastDataRead = true;
386     else
387       Dat.FastDataRead = false;
388     end
389     
390      case 'filechunksize'
391        Dat.FileChunkSize = round(varargin{ii+1});
392       
393      case 'wbar'
394        if strcmpi(varargin{ii+1},'on')
395          Dat.ShowWaitbar = 1;
396        else
397          Dat.ShowWaitbar = 0;
398        end
399         
400         case 'precision'
401           if strcmpi(varargin{ii+1},'single')
402                 Dat.precision = 'single';
403           end
404               
405          case 'flipkspace'
406       
407        flipkspace = varargin{ii+1};
408        if strcmpi(flipkspace,'off')
409          Dat.FlipKspace=0;
410        elseif strcmpi(flipkspace,'LR')
411          Dat.FlipKspace=1;
412        elseif strcmpi(flipkspace,'UD')
413          Dat.FlipKspace=2;
414        elseif strcmpi(flipkspace,'LRUD')
415          Dat.FlipKspace=3;
416        else
417          DATA=[];
418          if ~Dat.ShowErrors
419            msg_out = 'Bad value for property FlipKspace!';
420          else
421            error('Bad value for property FlipKspace!')
422          end
423          return
424        end
425       
426      case 'flipind'
427        Dat.FlipInd = varargin{ii+1};
428       
429      case 'outputfile'
430        Dat.OutputFile = varargin{ii+1};
431       
432      case 'reorientepi'
433        if strcmpi(varargin{ii+1},'on')
434          Dat.ReorientEPI = true;
435        end
436       
437      case 'removeepiphaseim'
438        if strcmpi(varargin{ii+1},'on')
439          Dat.RemoveEPIphaseIm = true;
440        end
441      case 'epiblocksize'
442        Dat.EPIBlockSize = round(varargin{ii+1});
443      case 'orientimages'
444        if strcmpi(varargin{ii+1},'off')
445          Dat.OrientImages = false;
446        end
447     otherwise
448      DATA=[];
449      if ~Dat.ShowErrors
450        msg_out = ['Unknown property "' varargin{ii} '"'];
451      else
452        error(['Unknown property "' varargin{ii} '"'])
453      end
454      return
455    end
456  end
457end
458
459% Parse filename
460[fpath,fname,fext]=fileparts(filename);
461if ~strcmpi(fname,'fid')
462  if isempty(fname)
463    fpath = [fpath,filesep];
464  else
465        if isempty(fpath)
466          fpath = [pwd,filesep,fname,fext,filesep];
467        else
468          fpath = [fpath,filesep,fname,fext,filesep];
469        end
470  end
471  fname = 'fid';
472else
473  fpath = [fpath,filesep];
474end
475
476%% Read procpar if not given as input argument
477if isempty(procpar) || nargin~=2
478  [procpar,proc_msg]=aedes_readprocpar([fpath,'procpar']);
479  if isempty(procpar)
480    DATA=[];
481    if ~Dat.ShowErrors
482      msg_out=proc_msg;
483    else
484      error(proc_msg)
485    end
486    return
487  end
488end
489
490%% Read phasetable if necessary
491if isfield(procpar,'petable') && isempty(Dat.phasetable) && ...
492    ~isempty(procpar.petable{1}) && ~strcmpi(procpar.petable{1},'n') && ...
493    Dat.Sorting
494  % Look in preferences for tablib-directory
495  try
496    tabpath = getpref('Aedes','TabLibDir');
497  catch
498    % If the table path was not found in preferences, try to use aedes
499    % directory
500    tmp_path = which('aedes');
501    if ~isempty(tmp_path)
502      aedes_path=fileparts(tmp_path);
503      tabpath = [aedes_path,filesep,'tablib',filesep];
504    else
505      % If all else fails, look in the current directory
506      tabpath = [pwd,filesep,'tablib',filesep];
507    end
508  end
509  [phasetable,msg] = aedes_readphasetable([tabpath,procpar.petable{1}]);
510 
511  % If petable cannot be found, check if it is in procpar...
512  if isempty(phasetable) && isfield(procpar,'pe_table')
513    phasetable = {'t1',procpar.pe_table(:)};
514  elseif isempty(phasetable) && isfield(procpar,'pelist') && ...
515      ~isempty(procpar.pelist) && isnumeric(procpar.pelist)
516    phasetable = {'t1',reshape(procpar.pelist,procpar.etl,[]).'};
517  end
518 
519  % If the sequence is a fast spin echo, try to construct the phasetable
520  if isempty(phasetable) && isfield(procpar,'petable') && ...
521      strncmpi(procpar.petable{1},'fse',3)
522    err_msg = 'Could not find phasetable!';
523    if ~Dat.ShowErrors
524      msg_out=err_msg;
525    else
526      error(err_msg)
527    end
528    return
529    %phasetable = {'t1',l_CreateFsemsPhaseTable(procpar)};
530  end
531 
532  %% Abort and throw error if phasetable cannot be read and the FID-file
533  % has not been sorted
534  if isempty(phasetable) && ~isfield(procpar,'flash_converted')
535    DATA=[];
536    if ~Dat.ShowErrors
537      msg_out={['Could not find the required phase table "' ...
538                procpar.petable{1} '" in the following folders'],...
539               ['"' tabpath '".']};
540    else
541      error('Could not find the required phase table "%s" in %s',...
542            procpar.petable{1},tabpath)
543    end
544    return
545  elseif ( isempty(phasetable) && isfield(procpar,'flash_converted') ) || ...
546      ~Dat.UsePhaseTable
547    %% If the FID-file has been sorted, the reading can continue but
548    % throw a warning.
549    fprintf(1,'Warning: aedes_readfid: Could not find phasetable "%s" in "%s"!\n',procpar.petable{1},tabpath)
550    Dat.phasetable=[];
551  else
552    Dat.phasetable = phasetable{1,2};
553  end
554end
555
556% Convert phasetable indices to MATLAB indices
557if ~isempty(Dat.phasetable)
558  Dat.phasetable=Dat.phasetable+(-min(min(Dat.phasetable))+1);
559else
560  Dat.UsePhaseTable = false;
561end
562
563%% Open FID-file
564[file_fid,msg]=fopen([fpath,fname],'r','ieee-be');
565if file_fid < 0
566  DATA=[];
567  if ~Dat.ShowErrors
568    msg_out=['Could not open file "' filename '" for reading.'];
569  else
570    error('Could not open file "%s" for reading.',filename)
571  end
572  return
573end
574
575% Read only header
576if ReadHdr && ~ReadData
577  [hdr,msg_out]=l_ReadDataFileHeader(file_fid);
578  if ~isempty(msg_out)
579    DATA=[];
580    fclose(file_fid);
581    if Dat.ShowErrors
582      error(msg_out)
583    end
584    return
585  else
586    DATA.HDR.FileHeader = hdr.FileHeader;
587    DATA.FTDATA=[];
588    DATA.KSPACE=[];
589    DATA.PROCPAR=[];
590    DATA.PHASETABLE=[];
591  end
592elseif ~ReadHdr && ReadData % Header structure given as input argument
593                            % Read only data.
594  [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar);
595  if ~isempty(msg_out)
596    DATA=[];
597    fclose(file_fid);
598    if Dat.ShowErrors
599      error(msg_out)
600    end
601    return
602  else
603    DATA.HDR.FileHeader=hdr.FileHeader;
604    DATA.HDR.BlockHeaders = hdr.BlockHeaders;
605    DATA.FTDATA=data;
606    DATA.KSPACE=kspace;
607    DATA.PROCPAR=procpar;
608    DATA.PHASETABLE=Dat.phasetable;
609  end
610elseif ReadHdr && ReadData  % Read headers and data
611  [hdr,msg_out]=l_ReadDataFileHeader(file_fid);
612  [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar);
613  if ~isempty(msg_out)
614    DATA=[];
615    fclose(file_fid);
616    if Dat.ShowErrors
617      error(msg_out)
618    end
619    return
620  else
621    DATA.HDR.FileHeader=hdr.FileHeader;
622    DATA.HDR.BlockHeaders = hdr.BlockHeaders;
623    DATA.FTDATA=data;
624    DATA.KSPACE=kspace;
625    DATA.PROCPAR=procpar;
626    DATA.PHASETABLE=Dat.phasetable;
627  end
628end
629
630% Set file name and path to the HDR structure
631DATA.HDR.fname=fname;
632DATA.HDR.fpath=fpath;
633
634% Set parameter values
635DATA.HDR.Param.ReturnKSpace = Dat.ReturnKSpace;
636DATA.HDR.Param.ReturnFTData = Dat.ReturnFTData;
637if Dat.ZeroPadding==0
638  DATA.HDR.Param.ZeroPadding = 'off';
639elseif Dat.ZeroPadding==1
640  DATA.HDR.Param.ZeroPadding = 'on';
641else
642  DATA.HDR.Param.ZeroPadding = 'auto';
643end
644if ~Dat.DCcorrection
645  DATA.HDR.Param.DCcorrection = 'off';
646else
647  DATA.HDR.Param.DCcorrection = 'on';
648end
649if Dat.Sorting==0
650  DATA.HDR.Param.Sorting = 'off';
651else
652  DATA.HDR.Param.Sorting = 'on';
653end
654if Dat.FlipKspace==0
655  DATA.HDR.Param.FlipKspace = 'off';
656elseif Dat.FlipKspace==1
657  DATA.HDR.Param.FlipKspace = 'LR';
658elseif Dat.FlipKspace==2
659  DATA.HDR.Param.FlipKspace = 'UD';
660elseif Dat.FlipKspace==3
661  DATA.HDR.Param.FlipKspace = 'LRUD';
662end
663DATA.HDR.Param.FlipInd = Dat.FlipInd;
664
665
666% Close file
667fclose(file_fid);
668
669%% Write slices to NIfTI files
670if ischar(Dat.OutputFile)
671  if isempty(Dat.OutputFile)
672    [fp,fn,fe]=fileparts(DATA.HDR.fpath(1:end-1));
673    fprefix=fn;
674    niipath = [pwd,filesep];
675  else
676    [fp,fn,fe]=fileparts(Dat.OutputFile);
677    if strcmpi(fe,'.nii')
678      fprefix=fn;
679    else
680      fprefix=[fn,fe];
681    end
682    if isempty(fp)
683      niipath = [pwd,filesep];
684    else
685      niipath = [fp,filesep];
686    end
687  end
688 
689  % Create file names
690  filenames={};
691  for ii=1:size(DATA.FTDATA,3)
692    filenames{ii}=sprintf('%s%03d%s',[fprefix,'_'],ii,'.nii');
693  end
694  nFiles = length(filenames);
695   
696  h=aedes_wbar(0,sprintf('Saving slice 1/%d in NIfTI format...',nFiles));
697  for ii=1:nFiles
698    aedes_wbar(ii/nFiles,h,sprintf('Saving slice %d/%d in NIfTI format...',ii, ...
699                             nFiles))
700    [done,msg]=aedes_write_nifti(DATA.FTDATA(:,:,ii),...
701                           [niipath,filenames{ii}],'DataType','single',...
702                           'FileType',2);
703  end
704  delete(h)
705 
706  if ~done
707    warning('Error occurred while writing NIfTI-files. Could not write file(s)!')
708  end
709 
710end
711return
712
713%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
714% Read Data File (Main) Header
715%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
716function [hdr,msg_out]=l_ReadDataFileHeader(file_fid)
717
718msg_out='';
719%% Read Data File Header
720try
721  FH.nblocks = fread(file_fid,1,'long');   % Number of blocks in file
722  FH.ntraces = fread(file_fid,1,'long');   % Number of traces per block
723  FH.np = fread(file_fid,1,'long');        % Number of elements per trace
724  FH.ebytes = fread(file_fid,1,'long');    % Number of bytes per element
725  FH.tbytes = fread(file_fid,1,'long');    % Number of bytes per trace
726  FH.bbytes = fread(file_fid,1,'long');    % Number of bytes per block
727  FH.vers_id = fread(file_fid,1,'short');  % Software version, file_id status bits
728  FH.status = fread(file_fid,1,'short');   % Status of whole file
729  FH.nbheaders = fread(file_fid,1,'long'); % Number of block headers per block
730 
731  hdr.FileHeader = FH;
732 
733  %% Parse status bits
734  hdr.FileHeader.status=[];
735  tmp_str = {'S_DATA',...          % 0 = no data, 1 = data
736             'S_SPEC',...          % 0 = FID, 1 = spectrum
737             'S_32',...            % 0 = 16 bit, 1 = 32 bit integer
738             'S_FLOAT',...         % 0 = integer, 1 = floating point
739             'S_COMPLEX',...       % 0 = real, 1 = complex
740             'S_HYPERCOMPLEX',...  % 1 = hypercomplex
741             '',...                % Unused bit
742             'S_ACQPAR',...        % 0 = not Acqpar, 1 = Acqpar
743             'S_SECND',...         % 0 = first FT, 1 = second FT
744             'S_TRANSF',...        % 0 = regular, 1 = transposed
745             '',...                % Unused bit
746             'S_NP',...            % 1 = np dimension is active
747             'S_NF',...            % 1 = nf dimension is active
748             'S_NI',...            % 1 = ni dimension is active
749             'S_NI2',...           % 1 = ni2 dimension is active
750             ''...                 % Unused bit
751            };
752  status_bits = fliplr(double(dec2bin(FH.status,16))==49);
753  for ii=1:length(tmp_str)
754    if ~isempty(tmp_str{ii})
755      hdr.FileHeader.status.(tmp_str{ii}) = status_bits(ii);
756    end
757  end
758catch
759  hdr=[];
760  msg_out=['Error occurred while reading file header from "' filename '".'];
761  return
762end
763
764
765%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
766% Read Block Headers and Data
767%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
768function [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar)
769
770hdr.BlockHeaders = [];
771%hdr.HypercmplxBHeaders = [];
772data=[];
773kspace=[];
774count = [];
775msg_out='';
776
777if isfield(procpar,'seqcon')
778  procpar.seqcon=procpar.seqcon{1};
779else
780  procpar.seqcon='nnnnn';
781end
782
783%% Determine acquisition type from seqcon parameter
784if all(procpar.seqcon=='n')          % 1D spectroscopy
785  Dat.AcqType=1;
786elseif all(procpar.seqcon(4:5)=='n') % 2D imaging
787  Dat.AcqType=2;
788elseif procpar.seqcon(5)=='n'        % 3D imaging
789  Dat.AcqType=3;
790else                         % 4D imaging
791  Dat.AcqType=4;
792end
793
794%% Double-check that AcqType is not 1D spectral
795if isfield(procpar,'nv') && procpar.nv==0
796  Dat.AcqType=1;
797end
798
799%% Use some heuristical approach to "triple-check" that the data is not
800%  1D spectral
801if ~isempty(strfind(procpar.seqfil{1},'STEAM')) || ...
802      ~isempty(strfind(procpar.seqfil{1},'steam')) || ...
803      ~isempty(strfind(procpar.seqfil{1},'LASER')) || ...
804      ~isempty(strfind(procpar.seqfil{1},'laser')) || ...
805      ~isempty(strfind(procpar.seqfil{1},'PRESS')) || ...
806      ~isempty(strfind(procpar.seqfil{1},'press')) || ...
807      ~isempty(strfind(procpar.seqfil{1},'1PULSE')) || ...
808      ~isempty(strfind(procpar.seqfil{1},'1pulse')) || ...
809      ~isempty(strfind(procpar.seqfil{1},'2PULSE')) || ...
810      ~isempty(strfind(procpar.seqfil{1},'2pulse'))
811  Dat.AcqType=1;
812end
813
814UsePhaseTable=Dat.UsePhaseTable;
815
816
817%% If the data has been converted with flashc, the seqcon parameter needs
818% to be changed
819if isfield(procpar,'flash_converted')
820  if isfield(procpar,'ni') && procpar.ni>1
821    procpar.seqcon(3)='s';
822  elseif ((procpar.seqcon(2)=='c') && (procpar.seqcon(3)=='c'))
823    procpar.seqcon(2)='s';
824  end
825  UsePhaseTable=false; % Do not try to order data if flashc has been run
826end
827
828%% Set number of transients
829if ~isfield(procpar,'nt')
830  nt=1;
831else
832  nt=procpar.nt;
833end
834
835%% Determine if the arrayed acquisition was used
836if isfield(procpar,'array') && ~isempty(procpar.array{1})
837 
838  if length(procpar.array)==1 && ~iscell(procpar.array{1}) && ...
839      strcmp(procpar.array{1},'pad') && all(procpar.pad==0)
840    % Skip the array parsing if the array is a dummy "pad"-array...
841    Dat.isAcqArrayed = false;
842    Dat.ArrayLength = 1;
843  else
844    Dat.isAcqArrayed = true;
845    Dat.ArrayLength = [];
846   
847    % Determine array length
848    for ii=1:length(procpar.array)
849      if iscell(procpar.array{ii})
850        Dat.ArrayLength(ii) = length(procpar.(procpar.array{ii}{1}));
851      else
852        Dat.ArrayLength(ii) = length(procpar.(procpar.array{ii}));
853      end
854    end
855    Dat.ArrayLength = prod(Dat.ArrayLength);
856  end
857else
858  Dat.isAcqArrayed = false;
859  Dat.ArrayLength = 1;
860end
861
862%% Determine if the data is EPI data
863if ( isfield(procpar,'readres') && isfield(procpar,'phaseres') ) || ...
864    ( isfield(procpar,'apptype') && strcmp(procpar.apptype{1},'imEPI') )
865  Dat.isEPIdata = true;
866else
867  Dat.isEPIdata = false;
868end
869
870BlockHeadStatusLabels = {'S_DATA',...       % 0 = no data, 1 = data
871                    'S_SPEC',...          % 0 = FID, 1 = spectrum
872                    'S_32',...            % 0 = 16 bit, 1 = 32 bit integer
873                    'S_FLOAT',...         % 0 = integer, 1 = floating point
874                    'S_COMPLEX',...       % 0 = real, 1 = complex
875                    'S_HYPERCOMPLEX',...  % 1 = hypercomplex
876                    '',...                % Unused bit
877                    'MORE_BLOCKS',...     % 0 = absent, 1 = present
878                    'NP_CMPLX',...        % 0 = real, 1 = complex
879                    'NF_CMPLX',...        % 0 = real, 1 = complex
880                    'NI_CMPLX',...        % 0 = real, 1 = complex
881                    'NI2_CMPLX',...       % 0 = real, 1 = complex
882                    '',...                % Unused bit
883                    '',...                % Unused bit
884                    '',...                % Unuesd bit
885                    ''...                 % Unused bit
886                   };
887
888BlockHeadModeLabels = {'NP_PHMODE',...   % 1 = ph mode
889                    'NP_AVMODE',...    % 1 = av mode
890                    'NP_PWRMODE',...   % 1 = pwr mode
891                    '',...             % Unused bit
892                    'NF_PHMODE',...    % 1 = ph mode
893                    'NF_AVMODE',...    % 1 = av mode
894                    'NF_PWRMODE',...   % 1 = pwr mode
895                    '',...             % Unused bit
896                    'NI_PHMODE',...    % 1 = ph mode
897                    'NI_AVMODE',...    % 1 = av mode
898                    'NI_PWRMODE',...   % 1 = pwr mode
899                    '',...             % Unused bit
900                    'NI2_PHMODE',...   % 1 = ph mode
901                    'NI2_AVMODE',...   % 1 = av mode
902                    'NI2_PWRMODE',...  % 1 = pwr mode
903                    ''...              % Unused bit
904                   };
905
906
907% The nbheaders -field is not correct in some cases. Thus, this field
908% cannot be trusted and the real nbheaders has to be calculated from the
909% byte counts.
910tmp_bytes=hdr.FileHeader.bbytes-hdr.FileHeader.tbytes*hdr.FileHeader.ntraces;
911header_bytes=28;
912if rem(tmp_bytes,header_bytes)~=0
913  msg_out = 'Block header byte count does not match with file header';
914  return
915else
916  nbheaders = tmp_bytes/28;
917end
918%nbheaders = hdr.FileHeader.nbheaders;
919
920%% Allocate space for k-space
921% kspace=zeros(hdr.FileHeader.np/2,...
922%              hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);
923
924if ~Dat.FastDataRead
925  if any(Dat.AcqType==[1 2])
926    switch procpar.seqcon(2:3)
927      case {'cc','sc'}
928        kspace = complex(zeros(hdr.FileHeader.np/2,...
929          hdr.FileHeader.ntraces,...
930          hdr.FileHeader.nblocks,Dat.precision));
931      otherwise
932        kspace = complex(zeros(hdr.FileHeader.np/2,...
933          hdr.FileHeader.nblocks,...
934          hdr.FileHeader.ntraces,Dat.precision));
935    end
936  else
937    kspace = complex(zeros(hdr.FileHeader.np/2,...
938      hdr.FileHeader.ntraces,...
939      hdr.FileHeader.nblocks,Dat.precision));
940  end
941else
942  %kspace = [];
943   kspace = complex(zeros(hdr.FileHeader.np/2*hdr.FileHeader.ntraces,...
944     hdr.FileHeader.nblocks,Dat.precision));
945end
946
947%% - The older robust (but also slower) way for reading the data.
948%% When the blocksize is relatively small, this is also quite fast.
949if ~Dat.FastDataRead
950
951  % Initialize waitbar
952  if Dat.ShowWaitbar
953        wb_h = aedes_wbar(0/hdr.FileHeader.nblocks,...
954          {['Reading ',num2str(Dat.AcqType),'D VNMR data (seqcon: "' procpar.seqcon '")'],...
955          ['(Processing data block ' ...
956          num2str(0) '/' num2str(hdr.FileHeader.nblocks) ')']});
957  end
958
959  %% Read data and block headers
960  for ii=1:hdr.FileHeader.nblocks
961        %% Update waitbar
962        if Dat.ShowWaitbar
963          aedes_wbar(ii/hdr.FileHeader.nblocks,...
964                wb_h,...
965                {['Reading ',num2str(Dat.AcqType),'D VNMR data (seqcon: "' procpar.seqcon '")'],...
966                ['(Processing data block ' ...
967                num2str(ii) '/' num2str(hdr.FileHeader.nblocks) ')']})
968        end
969
970        %% Read block header and hypercomplex header
971        for kk=1:nbheaders
972          %% Read block header
973          if kk==1
974                hdr.BlockHeaders.scale = fread(file_fid,1,'short');  % Scaling factor
975                tmp_status = fread(file_fid,1,'short'); % Status of data in block
976                hdr.BlockHeaders.status = [];
977                hdr.BlockHeaders.index = fread(file_fid,1,'short');  % Block index
978                tmp_mode = fread(file_fid,1,'short');   % Mode of data in block
979                hdr.BlockHeaders.mode = [];
980                hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
981                hdr.BlockHeaders.lpval = fread(file_fid,1,'float');  % f2 (2D-f1) left phase in phasefile
982                hdr.BlockHeaders.rpval = fread(file_fid,1,'float');  % f2 (2D-f1) right phase in phasefile
983                hdr.BlockHeaders.lvl = fread(file_fid,1,'float');    % level drift correction
984                hdr.BlockHeaders.tlt = fread(file_fid,1,'float');    % tilt drift correction
985
986                %% Parse status and mode bits
987                status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
988                mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
989                for tt=1:length(BlockHeadStatusLabels)
990                  if ~isempty(BlockHeadStatusLabels{tt})
991                        hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
992                  end
993                  if ~isempty(BlockHeadModeLabels{tt})
994                        hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
995                  end
996                end
997
998
999          else %% Read hypercomplex header
1000                fread(file_fid,1,'short'); % Spare
1001                hdr.BlockHeaders.HCHstatus = fread(file_fid,1,'short');
1002                fread(file_fid,1,'short'); % Spare
1003                fread(file_fid,1,'short'); % Spare
1004                fread(file_fid,1,'long'); % Spare
1005                hdr.BlockHeaders.HCHlpval1 = fread(file_fid,1,'float');
1006                hdr.BlockHeaders.HCHrpval1 = fread(file_fid,1,'float');
1007                fread(file_fid,1,'float'); % Spare
1008                fread(file_fid,1,'float'); % Spare
1009          end
1010        end
1011       
1012        %% Check block index to be sure about the data type
1013        if hdr.BlockHeaders.index~=ii
1014          fprintf(1,'Warning: Index mismatch in "%s" block %d\n',fopen(file_fid),ii);
1015         
1016          % Use information from the file header instead of the BlockHeader if
1017          % there is a mismatch in blockheader index...
1018          useFileHeader = true;
1019        else
1020          useFileHeader = false;
1021        end
1022       
1023        %% Determine data precision
1024        if useFileHeader
1025          if hdr.FileHeader.status.S_FLOAT==1
1026                prec_str = ['single=>',Dat.precision];
1027          elseif hdr.FileHeader.status.S_32==1 ...
1028                  && hdr.FileHeader.status.S_FLOAT==0
1029                prec_str = ['int32=>',Dat.precision];
1030          elseif hdr.FileHeader.status.S_32==0 ...
1031                  && hdr.FileHeader.status.S_FLOAT==0
1032                prec_str = ['int16=>',Dat.precision];
1033          end
1034         
1035        else
1036          if hdr.BlockHeaders.status.S_FLOAT==1
1037                prec_str = ['single=>',Dat.precision];
1038          elseif hdr.BlockHeaders.status.S_32==1 ...
1039                  && hdr.BlockHeaders.status.S_FLOAT==0
1040                prec_str = ['int32=>',Dat.precision];
1041          elseif hdr.BlockHeaders.status.S_32==0 ...
1042                  && hdr.BlockHeaders.status.S_FLOAT==0
1043                prec_str = ['int16=>',Dat.precision];
1044          end
1045        end
1046
1047        % Read k-space
1048        tmp=fread(file_fid,...
1049          [hdr.FileHeader.np,hdr.FileHeader.ntraces],...
1050          prec_str);
1051
1052
1053        % Store complex block and perform DC correction
1054        if ~Dat.DCcorrection || ( nt(1)>1 )
1055          complex_block = complex(tmp(1:2:end,:),tmp(2:2:end,:));
1056        else
1057          complex_block = complex(tmp(1:2:end,:)-hdr.BlockHeaders.lvl,...
1058                tmp(2:2:end,:)-hdr.BlockHeaders.tlt);
1059        end
1060
1061
1062        %% Store and order k-space values
1063        if any(Dat.AcqType==[1 2])
1064          switch procpar.seqcon(2:3)
1065                case {'cc','sc'}
1066                  kspace(:,:,ii) = complex_block;
1067                otherwise
1068                  kspace(:,ii,:) = complex_block;
1069          end
1070        else
1071          kspace(:,:,ii) = complex_block;
1072        end
1073
1074        % Do not save blockheaders by default. When reading data files with a lot of
1075        % blocks (e.g. over 1000) the processing of the DATA structure can be
1076        % slowed down considerably. If you for some reason want to save also the
1077        % block headers in the DATA structure comment out the code line below.
1078        hdr.BlockHeaders = [];
1079  end % for ii=1:hdr.
1080
1081  if Dat.ShowWaitbar
1082        delete(wb_h)
1083  end
1084 
1085else
1086  %% -------------------------------------------------------------------
1087  %% Fast Method for reading data. This may fail with some datas and can
1088  %% also require a relatively large amount of memory. This
1089  %% method should be used for EPI datas that contain a large number
1090  %% of block headers...
1091 
1092  % Check the size of the FID-file
1093  d=dir(fopen(file_fid));
1094  file_sz = d.bytes/1024/1024; % File size in MB
1095  if file_sz<Dat.FileChunkSize
1096    nBlocks = 1;
1097  else
1098    nBlocks = ceil(file_sz/Dat.FileChunkSize); % Read data in 500 MB blocks
1099  end
1100 
1101  % Initialize waitbar
1102  if Dat.ShowWaitbar
1103    if nBlocks==1
1104      wb_h = aedes_calc_wait(['Reading ',num2str(Dat.AcqType),...
1105        'D VNMR data (seqcon: "' procpar.seqcon '")']);
1106    else
1107      wb_h = aedes_wbar(1/nBlocks,{['Reading ',num2str(Dat.AcqType),...
1108        'D VNMR data (seqcon: "' procpar.seqcon '")'],...
1109        sprintf('Reading block 0/%d',nBlocks)});
1110    end
1111  end
1112 
1113  % The first block header is read and it is assumed that the values in
1114  % the other block headers don't change. 
1115  hdr.BlockHeaders.scale = fread(file_fid,1,'short');  % Scaling factor
1116  tmp_status = fread(file_fid,1,'short'); % Status of data in block
1117  hdr.BlockHeaders.status = [];
1118  hdr.BlockHeaders.index = fread(file_fid,1,'short');  % Block index
1119  tmp_mode = fread(file_fid,1,'short');   % Mode of data in block
1120  hdr.BlockHeaders.mode = [];
1121  hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID
1122  hdr.BlockHeaders.lpval = fread(file_fid,1,'float');  % f2 (2D-f1) left phase in phasefile
1123  hdr.BlockHeaders.rpval = fread(file_fid,1,'float');  % f2 (2D-f1) right phase in phasefile
1124  hdr.BlockHeaders.lvl = fread(file_fid,1,'float');    % level drift correction
1125  hdr.BlockHeaders.tlt = fread(file_fid,1,'float');    % tilt drift correction
1126 
1127  %% Parse status and mode bits
1128  status_bits = fliplr(double(dec2bin(tmp_status,16))==49);
1129  mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49);
1130  for tt=1:length(BlockHeadStatusLabels)
1131    if ~isempty(BlockHeadStatusLabels{tt})
1132      hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt);
1133    end
1134    if ~isempty(BlockHeadModeLabels{tt})
1135      hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt);
1136    end
1137  end
1138 
1139  %% Determine data precision
1140  if hdr.BlockHeaders.status.S_FLOAT==1
1141    prec_str = ['single=>',Dat.precision];
1142    prec = 4; % Precision in bytes
1143  elseif hdr.BlockHeaders.status.S_32==1 ...
1144      && hdr.BlockHeaders.status.S_FLOAT==0
1145    prec_str = ['int32=>',Dat.precision];
1146    prec = 4;
1147  elseif hdr.BlockHeaders.status.S_32==0 ...
1148      && hdr.BlockHeaders.status.S_FLOAT==0
1149    prec_str = ['int16=>',Dat.precision];
1150    prec = 2;
1151  end
1152 
1153  % Seek the file back to the beginning of the first block header
1154  fseek(file_fid,-28,0);
1155 
1156  % Determine the number of values that will result from block header(s)
1157  nVals = (nbheaders*28)/prec;
1158 
1159  nbh = floor(hdr.FileHeader.nblocks/nBlocks);
1160  szh = nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces;
1161  for ii=1:nBlocks
1162    if nBlocks~=1
1163      aedes_wbar(ii/nBlocks,wb_h,{['Reading ',num2str(Dat.AcqType),...
1164        'D VNMR data (seqcon: "' procpar.seqcon '")'],...
1165        sprintf('Reading block %d/%d',ii,nBlocks)});
1166    end
1167   
1168    % Read the whole data including block headers etc...
1169    if ii==nBlocks
1170      tmp = fread(file_fid,inf,prec_str);
1171    else
1172      tmp = fread(file_fid,nbh*szh,prec_str);
1173    end
1174    tmp=reshape(tmp,nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces,[]);
1175    tmp(1:nVals,:)=[];
1176   
1177    if ii==nBlocks
1178      inds = ((ii-1)*nbh+1):size(kspace,2);
1179    else
1180      inds = ((ii-1)*nbh+1):ii*nbh;
1181    end
1182   
1183    % Do DC-correction if necessary
1184    if ~Dat.DCcorrection || ( nt(1)>1 )
1185      kspace(:,inds)=complex(tmp(1:2:end,:,:),tmp(2:2:end,:,:));
1186    else
1187      kspace(:,inds)=complex(tmp(1:2:end,:,:)-hdr.BlockHeaders.lvl,...
1188        tmp(2:2:end,:,:)-hdr.BlockHeaders.tlt);
1189    end
1190  end
1191  clear tmp
1192 
1193  % Transform to 3D matrix
1194  kspace = reshape(kspace,hdr.FileHeader.np/2,...
1195    hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);
1196 
1197  %% Store and order k-space values
1198  if any(Dat.AcqType==[1 2])
1199    switch procpar.seqcon(2:3)
1200      case {'cc','sc'}
1201        %kspace(:,:,ii) = complex_block;
1202      otherwise
1203        %kspace(:,ii,:) = complex_block;
1204        kspace = permute(kspace,[1 3 2]);
1205    end
1206  else
1207    %kspace(:,:,ii) = complex_block;
1208  end
1209 
1210  if Dat.ShowWaitbar
1211    delete(wb_h)
1212    pause(0.1)
1213  end
1214 
1215end
1216
1217% Remove singleton dimensions from kspace
1218kspace = squeeze(kspace);
1219
1220% Check if raw kspace should be returned
1221if Dat.ReturnRawKspace
1222  %% Delete waitbar
1223  if Dat.ShowWaitbar && ishandle(wb_h)
1224        delete(wb_h)
1225  end
1226  return
1227end
1228
1229%% Support for RASER sequence ---------------------------
1230if isfield(procpar,'teType')
1231 
1232  if strcmpi(procpar.sing_sh,'y') && strcmpi(procpar.teType,'c')
1233       
1234    % Separate reference scans from the data
1235    if isfield(procpar,'refscan') && strcmp(procpar.refscan,'y')
1236      kspace=permute(reshape(kspace,procpar.np/2,procpar.ne/2,2,[]),[1 2 4 3]);
1237      noise = kspace(:,:,:,1);
1238      kspace = kspace(:,:,:,2);
1239      % DC-correction
1240      dc_level = squeeze(mean(mean(noise)));
1241      for ii=1:size(kspace,3)
1242        kspace(:,:,ii)=kspace(:,:,ii)-dc_level(ii);
1243      end
1244    end
1245   
1246    if strcmpi(procpar.readType,'s')
1247      sw = procpar.sw;
1248      dwell = cumsum(ones(1,size(kspace,1))*(1/sw));
1249    else
1250      error('This feature has not been implemented yet!!!')
1251    end
1252    dwell = dwell*1000;
1253    %kspace = permute(kspace,[2 1 3]);
1254   
1255   
1256   
1257    % Phase correction
1258    dims = size(kspace);
1259    nI = size(kspace,3);
1260    nv = size(kspace,4);
1261    ne = size(kspace,2);
1262    np = size(kspace,1);
1263    g = 42.58;
1264    G = procpar.gvox2;
1265    Ds = procpar.vox2/100;
1266    pos2 = procpar.pos2/10;
1267    fudge = 0;
1268    fudge2 = 1;
1269   
1270    tt_dwell=repmat(dwell(:),1,ne);
1271    tt_ne = repmat((1:ne),length(dwell),1);
1272    tt_np = repmat((1:np)',1,ne);
1273   
1274    i=sqrt(-1);
1275    phi = (-2.*pi.*g.*G.*Ds.*tt_dwell.*(tt_ne./ne - 1/2 - pos2/Ds))+fudge.*tt_np;
1276    phi = repmat(phi,[1,1,size(kspace,3)]);
1277    kspace = abs(kspace).*exp(i*(angle(kspace)+phi));
1278   
1279    % Flip every other column
1280    sz(1)=size(kspace,1);
1281    sz(2)=size(kspace,2);
1282    sz(3)=size(kspace,3);
1283    sz(4)=size(kspace,4);
1284    kspace=reshape(kspace,sz(1),prod(sz(2:4)));
1285    kspace(:,1:2:end) = flipud(kspace(:,1:2:end));
1286    kspace = reshape(kspace,sz);
1287   
1288    %kspace = reshape(permute(kspace,[2 1 3]),sz(2),[],1);
1289    %kspace(:,1:2:end) = flipud(kspace(:,1:2:end));
1290    %kspace=permute(reshape(kspace,sz(2),[],sz(3)),[2 1 3]);
1291   
1292  else
1293    data = [];
1294    error('RASER sequence of this type has not been implemeted yet!')
1295  end
1296 
1297  % Reshape into 4D matrix
1298  kspace = reshape(kspace,size(kspace,1),size(kspace,2),[],size(kspace,3));
1299 
1300  % Reorient if requested
1301  if Dat.OrientImages
1302    kspace = flipdim(kspace,2);
1303  end
1304 
1305  data = abs(fftshift(fft(fftshift(fft(kspace,[],3),3),[],1),1));
1306 
1307  % Delete kspace if not returned
1308  if ~Dat.ReturnKSpace
1309    kspace=[];
1310  end
1311 
1312  return
1313 
1314end
1315
1316%% Support for fat/water measurements ----------------------
1317if isfield(procpar,'seqfil') && strcmpi(procpar.seqfil,'ge3d_csi2')
1318 
1319  % Split kspace into fat and water (in 4th dimesion)
1320  kspace=reshape(permute(reshape(kspace,256,2,[]),[1 3 4 2]),256,128,[],2);
1321 
1322  % Fourier transform data
1323  if any(Dat.ZeroPadding==[1 2])
1324    data_sz = [procpar.np/2,procpar.nf,procpar.nv2,2];
1325    data = zeros(data_sz,class(kspace));
1326    data(:,:,:,1) = abs(fftshift(fftn(kspace(:,:,:,1),data_sz(1:3))));
1327    data(:,:,:,2) = abs(fftshift(fftn(kspace(:,:,:,2),data_sz(1:3))));
1328  else
1329    data = zeros(size(kspace),class(kspace));
1330    data(:,:,:,1) = abs(fftshift(fftn(kspace(:,:,:,1))));
1331    data(:,:,:,2) = abs(fftshift(fftn(kspace(:,:,:,2))));
1332  end
1333 
1334  % Delete kspace if not returned
1335  if ~Dat.ReturnKSpace
1336    kspace=[];
1337  end
1338 
1339  return
1340end
1341
1342% Check the number of receivers (PI stuff)
1343if isfield(procpar,'rcvrs') && ~isempty(procpar.rcvrs) && ...
1344    length(procpar.rcvrs{1})>1
1345  % Multiple receivers used
1346  if isfield(procpar,'apptype') && strcmp(procpar.apptype{1},'imEPI')
1347    % Store multiple receiver data in EPI measurements in 5th dimension
1348    % and calculate sum-of-squares image
1349    nRcv = length(find(procpar.rcvrs{1}=='y'));
1350    if ~isfield(procpar,'epiref_type') || strcmpi(procpar.epiref_type,'single')
1351      nRef = 1;
1352    elseif strcmpi(procpar.epiref_type,'triple')
1353      nRef = 3;
1354    end
1355    nVols = size(kspace,3)/nRcv-1;
1356    data = zeros(procpar.nv,procpar.np/2,procpar.ns,nVols+1,'single');
1357    kssz=size(kspace);
1358    blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
1359    nBlocks = ceil((size(kspace,3)/nRcv-nRef)/blksz);
1360    lnum = length(num2str(nBlocks));
1361    lnumstr = num2str(lnum);
1362    bsl = lnum*2+1;
1363    fprintf(1,'Processing data in blocks of %d volumes\n',blksz)
1364    fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks);
1365    for ii=1:nBlocks
1366      fprintf(1,repmat('\b',1,bsl));
1367      fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks);
1368      tmp_data = [];
1369      for kk=1:nRcv
1370        inds = [kk:nRcv:nRcv*nRef ((ii-1)*blksz*nRcv+nRcv*nRef+kk):nRcv:min((nRcv*ii*blksz+kk),kssz(3))];
1371        tmp_kspace = l_ReconstructKspace(kspace(:,:,inds),procpar,Dat);
1372        tmp_data(:,:,:,:,kk) = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
1373      end
1374      if ii==1
1375        data(:,:,:,1:size(tmp_data,4)) = sqrt(sum(tmp_data.*conj(tmp_data),5));
1376      elseif ii==nBlocks
1377        tmp_data(:,:,:,1:nRef,:)=[];
1378        data(:,:,:,((ii-1)*blksz+1+nRef):(nVols+1)) = sqrt(sum(tmp_data.*conj(tmp_data),5));
1379      else
1380        tmp_data(:,:,:,1:nRef,:)=[];
1381        data(:,:,:,((ii-1)*blksz+1+nRef):(ii*blksz+nRef)) = sqrt(sum(tmp_data.*conj(tmp_data),5));
1382      end
1383    end
1384    fprintf(1,'\n')
1385   
1386%     for kk=1:(size(kspace,3)/nRcv-1)
1387%       for ii=1:nRcv
1388%         tmp_kspace = l_ReconstructKspace(kspace(:,:,[ii kk*nRcv+ii]),procpar,Dat);
1389%         tmp_data(:,:,:,:,ii) = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
1390%         
1391%         %tmp_kspace = l_ReconstructKspace(kspace(:,:,ii:nRcv:end),procpar,Dat);
1392%         %kspace2(:,:,:,:,ii)=tmp_kspace;
1393%       end
1394%       if kk==1
1395%         data = sqrt(mean(tmp_data.*conj(tmp_data),5));
1396%       else
1397%         data(:,:,:,kk+1) = sqrt(mean(tmp_data(:,:,:,2,:).*conj(tmp_data(:,:,:,2,:)),5));
1398%       end
1399%       fprintf(1,repmat('\b',1,bsl));
1400%       fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],kk,nVols);
1401%     end
1402%     fprintf(1,'\n')
1403%     %kspace = kspace2;
1404%     %data = fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2);
1405%     %data = sqrt(mean(data.*conj(data),5));
1406%     %kspace2=[];
1407   
1408    % Remove reference image if requested
1409    if Dat.isEPIdata && Dat.RemoveEPIphaseIm
1410      %data = data(:,:,:,2:end);
1411      data(:,:,:,1:nRef)=[];
1412    end
1413   
1414    if Dat.OrientImages && ~isempty(procpar) && ...
1415        isfield(procpar,'orient') && any(Dat.AcqType==[2 3 4])
1416      orient = procpar.orient{1};
1417      if any(strcmpi(orient,{'xyz','trans90','cor90','sag90'}))
1418        data = flipdim(aedes_rot3d(data,1,3),2);
1419      elseif strcmpi(orient,'oblique')
1420        data = flipdim(flipdim(data,1),2);
1421      else
1422        data = flipdim(flipdim(data,1),2);
1423      end
1424    end
1425   
1426  else
1427    nRcv = length(find(procpar.rcvrs{1}=='y'));
1428    data = [];
1429    kspace2 = [];
1430    for ii=1:nRcv
1431      tmp_kspace = l_ReconstructKspace(kspace(:,:,ii),procpar,Dat);
1432      kspace2(:,:,:,:,ii)=tmp_kspace;
1433      if Dat.ReturnFTData
1434        tmp_data = l_CalculateFFT(tmp_kspace,procpar,Dat);
1435        data(:,:,:,:,ii) = tmp_data;
1436      end
1437    end
1438    kspace = kspace2;
1439    kspace2=[];
1440  end
1441else
1442  % Only one receiver
1443  kspace = l_ReconstructKspace(kspace,procpar,Dat);
1444  data = l_CalculateFFT(kspace,procpar,Dat);
1445end
1446
1447% Delete kspace if not returned
1448if ~Dat.ReturnKSpace
1449  kspace=[];
1450end
1451
1452%===============================================
1453% Reconstruct kspace for different sequences
1454%===============================================
1455function kspace=l_ReconstructKspace(kspace,procpar,Dat)
1456
1457
1458%% Flip images for certain measurements
1459if Dat.FlipKspace~=0
1460  if ischar(Dat.FlipInd)
1461    if strcmpi(Dat.FlipInd,'all')
1462      flipind = 1:size(kspace,3);
1463    elseif strcmpi(Dat.FlipInd,'alt')
1464      flipind = 2:2:size(kspace,3);
1465    end
1466  else
1467    flipind = Dat.FlipInd;
1468  end
1469 
1470  for ii=flipind
1471    if Dat.FlipKspace==1
1472      kspace(:,:,ii)=fliplr(kspace(:,:,ii));
1473    elseif Dat.FlipKspace==2
1474      kspace(:,:,ii)=flipud(kspace(:,:,ii));
1475    else
1476      kspace(:,:,ii)=flipud(fliplr(kspace(:,:,ii)));
1477    end
1478  end
1479end
1480
1481% Handle arrayed and RARE sequences
1482if Dat.isAcqArrayed && Dat.AcqType~=1 && Dat.Sorting && ~Dat.isEPIdata
1483  %if ~strcmpi(procpar.seqcon(2:3),'cc') && ( isempty(Dat.phasetable) | ...
1484  %        isfield(procpar,'flash_converted') )
1485  if (( isempty(Dat.phasetable) && ~strcmpi(procpar.seqcon(2:3),'sc') ) && ...
1486      ~(strcmpi(procpar.seqcon(2:3),'cc') && Dat.ArrayLength==size(kspace,3))) || ...
1487      isfield(procpar,'flash_converted')
1488                                 
1489    % Sort uncompressed arrayed data
1490    ks_order = 1:procpar.nv*Dat.ArrayLength;
1491    tmp = 0:procpar.nv:(procpar.ne-1)*procpar.nv;
1492    tmp=repmat(tmp,procpar.nv*Dat.ArrayLength,1);
1493    ks_order = reshape(ks_order,Dat.ArrayLength,procpar.nv).';
1494    ks_order = ks_order(:);
1495    ks_order = repmat(ks_order,1,procpar.ne);
1496    ks_order = ks_order+tmp;
1497    ks_order = ks_order.';
1498   
1499    if procpar.ns==1
1500      kspace=kspace(:,ks_order);
1501    else
1502      kspace = kspace(:,ks_order,:);
1503    end
1504   
1505    % Reshape into 3D matrix
1506    kspace=reshape(kspace,[size(kspace,1) procpar.nv Dat.ArrayLength* ...
1507                        procpar.ns]);
1508  else
1509    %% Sort RARE type data
1510    if Dat.UsePhaseTable && ~isempty(Dat.phasetable)
1511      %Dat.phasetable = Dat.phasetable';
1512      if Dat.isAcqArrayed && Dat.ArrayLength>1 && strcmpi(procpar.seqcon(2:3),'cs')
1513        kspace = permute(reshape(reshape(kspace,size(kspace,1),[]),size(kspace,1),Dat.ArrayLength,[]),[1 3 2]);
1514      else
1515        Dat.phasetable = Dat.phasetable.';
1516      end
1517      kspace(:,Dat.phasetable(:),:) = kspace;
1518    end
1519  end
1520elseif Dat.isEPIdata
1521  %% Support for EPI-measurements
1522 
1523 
1524  % EPI data is measured with old VNMR
1525  if isfield(procpar,'readres') && isfield(procpar,'phaseres')
1526   
1527    % Number of slices
1528    tmp_ns=length(procpar.pss);
1529   
1530    if isfield(procpar,'navecho') && strcmpi(procpar.navecho{1},'y')
1531      tmp_nv = procpar.nv-procpar.nseg;
1532    else
1533      tmp_nv = procpar.nv;
1534    end
1535    kspace = reshape(kspace,[size(kspace,1) ...
1536      size(kspace,2)/tmp_nv tmp_nv]);
1537    kspace = permute(kspace,[1 3 2]);
1538   
1539    % Reshape to 4D matrix
1540    kspace = reshape(kspace,[size(kspace,1) size(kspace,2) ...
1541      tmp_ns size(kspace,3)/tmp_ns]);
1542   
1543   
1544  elseif isfield(procpar,'apptype') && strcmp(procpar.apptype{1},'imEPI')
1545    % If EPI data is measured with new VNMRj system
1546   
1547    % Number of slices
1548    ns = length(procpar.pss);
1549   
1550    % Reshape kspace to 4D matrix
1551    if isfield(procpar,'fract_ky') && procpar.fract_ky~=procpar.nv/2
1552      kspace = reshape(kspace,procpar.np/2,procpar.nv/2+procpar.fract_ky,...
1553        ns,[]);
1554      kspace(procpar.np/2,procpar.nv,1)=0;
1555    else
1556      kspace = reshape(kspace,procpar.np/2,procpar.nv,...
1557        ns,[]);
1558    end
1559   
1560    % Flip even kspace lines
1561    if ~Dat.FastDataRead
1562      for tt=2:2:size(kspace,2)
1563        kspace(:,tt,:,:) = flipdim(kspace(:,tt,:,:),1);
1564      end
1565    else
1566      kspace(:,2:2:end,:,:) = flipdim(kspace(:,2:2:end,:,:),1);
1567    end
1568   
1569   
1570    % Sort centric measurements
1571    if isfield(procpar,'ky_order') && strcmpi(procpar.ky_order,'c')
1572      kspace(:,Dat.phasetable(:),:,:)=kspace;
1573    end
1574    %kspace = flipdim(kspace,1);
1575   
1576    % EPI phase correction -------------------------
1577    if ~isfield(procpar,'epiref_type') || strcmpi(procpar.epiref_type,'single')
1578      % Single reference pointwise phase correction
1579     
1580      % Get the reference image
1581      ref_im = kspace(:,:,:,1);
1582     
1583      % Do phase correction.
1584      phase_e = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
1585      for kk=2:size(kspace,4)
1586        kspace(:,:,:,kk) = ifft(fft(kspace(:,:,:,kk),[],1).*phase_e,[],1);
1587      end
1588    elseif strcmpi(procpar.epiref_type,'triple')
1589      % Triple reference pointwise phase correction
1590     
1591      % Get the reference images
1592      ref1 = kspace(:,:,:,1);
1593      phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
1594      ref2 = flipdim(kspace(:,:,:,3),1);
1595      phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
1596      im1 = flipdim(kspace(:,:,:,2),1);
1597     
1598      % Correct phase for reversed read gradient image
1599      rev_phase = fft(im1,[],1).*phase2;
1600      for kk=4:size(kspace,4)
1601        kspace(:,:,:,kk)=ifft(rev_phase+fft(kspace(:,:,:,kk),[],1).*phase1);
1602      end
1603    end
1604  end
1605 
1606else
1607  % This section contains various "chewing gum and ironwire" type
1608  % patches that hopefully work...
1609 
1610  if strcmp(procpar.seqcon,'nncnn')
1611    kspace = permute(kspace,[1 3 2]);
1612  elseif strcmp(procpar.seqcon,'nccnn') && length(size(kspace))==2 && ...
1613      procpar.ns>=1 && Dat.AcqType~=1
1614    if ~isempty(Dat.phasetable)
1615      kssz = size(kspace);
1616      phsz = size(Dat.phasetable);
1617      kspace=permute(reshape(kspace,procpar.np/2,...
1618        phsz(2),...
1619        kssz(2)/phsz(2)),[1 3 2]);
1620      kspace = permute(reshape(kspace,...
1621        procpar.np/2,procpar.ns,kssz(2)/procpar.ns),[1 3 2]);
1622      if Dat.Sorting
1623        kspace(:,Dat.phasetable(:),:)=kspace;
1624      end
1625    else
1626      kspace = reshape(kspace,...
1627        [procpar.np/2 procpar.ns ...
1628        size(kspace,2)/procpar.ns]);
1629      kspace=permute(kspace,[1 3 2]);
1630    end
1631   
1632  elseif strcmp(procpar.seqcon,'nscsn') && length(size(kspace))==3 && ...
1633          Dat.AcqType~=1
1634        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1635        % Support for 3D fast spin-echo
1636        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1637        if ~isempty(Dat.phasetable)
1638          for ii=1:size(kspace,3)
1639                tmp = kspace(:,:,ii);
1640                kssz = size(tmp);
1641                phsz = size(Dat.phasetable);
1642                tmp=permute(reshape(tmp,procpar.np/2,...
1643                  phsz(2),...
1644                  kssz(2)/phsz(2)),[1 3 2]);
1645                tmp = permute(reshape(tmp,...
1646                  procpar.np/2,procpar.ns,kssz(2)/procpar.ns),[1 3 2]);
1647                if Dat.Sorting
1648                  tmp(:,Dat.phasetable(:),:)=tmp;
1649                end
1650                kspace(:,:,ii)=tmp;
1651          end
1652        end
1653       
1654       
1655  elseif Dat.AcqType~=1 && isfield(procpar,'nv')
1656    if length(size(kspace))==2 && ...
1657              ( size(kspace,1)~=(procpar.np/2) || size(kspace,2)~=procpar.nv )
1658      %% Reshape the kspace to the appropriate size
1659      kspace=reshape(kspace,[procpar.np/2 procpar.nv size(kspace,2)/ ...
1660                          procpar.nv]);
1661    elseif length(size(kspace))==3
1662      kspace = reshape(kspace,procpar.np/2,procpar.nv,[]);
1663      if Dat.Sorting && ~isempty(Dat.phasetable)
1664        kspace(:,Dat.phasetable(:),:) = kspace;
1665      end
1666    end
1667  end
1668 
1669 
1670  %%% Support for Teemu's ASE3D-data
1671  %if strcmpi(procpar.seqcon,'ncccn')
1672  %  %% Reshape the kspace to the appropriate size
1673  %  kspace=reshape(kspace,[procpar.np/2 procpar.nv size(kspace,2)/procpar.nv]);
1674  %end
1675end
1676
1677
1678
1679%=========================================
1680% Fourier Transform Data
1681%=========================================
1682function data=l_CalculateFFT(kspace,procpar,Dat)
1683data=[];
1684
1685% Return image/spectral data --------------------------------
1686if Dat.ReturnFTData
1687  %% Fourier transform spectral data
1688  if Dat.AcqType==1
1689        wb_h = aedes_wbar(0,'Fourier transforming spectral data');
1690    %data=zeros(size(kspace),'single');
1691        data=kspace;
1692    sz=size(kspace,2)*size(kspace,3);
1693    count=1;
1694    for kk=1:size(kspace,3)
1695      for ii=1:size(kspace,2)
1696        %% Update waitbar
1697        if Dat.ShowWaitbar
1698          aedes_wbar(count/sz,...
1699               wb_h,...
1700               {'Fourier transforming spectral data',...
1701                ['(Processing spectrum ' num2str(count) '/' num2str(sz) ')']})
1702                end
1703        data(:,ii,kk) = abs(fftshift(fft(data(:,ii,kk))));
1704        count=count+1;
1705      end
1706    end
1707 
1708  %% Fourier transform image data
1709  else
1710   
1711    %% Zeropadding auto
1712    if Dat.ZeroPadding==2
1713     
1714      if Dat.isEPIdata && isfield(procpar,'readres') && isfield(procpar,'phaseres')
1715       
1716        %% Determine the image size for EPI data for procpar "readres"
1717        %% and "phaseres" fields
1718        data_sz = [procpar.readres procpar.phaseres];
1719        if isequal(data_sz,[size(kspace,1),size(kspace,2)])
1720          DoZeroPadding=false;
1721        else
1722          DoZeroPadding=true;
1723        end
1724        data_sz = [procpar.readres ...
1725          procpar.phaseres size(kspace,3) size(kspace,4)];
1726      else
1727       
1728        %% Check if zeropadding is necessary.
1729        ks_sz_sorted = sort([size(kspace,1),size(kspace,2)]);
1730        fov_sz_sorted = sort([procpar.lro,procpar.lpe]);
1731        sliceRelativeDim = ks_sz_sorted(2)/ks_sz_sorted(1);
1732        FOVrelativeDim = fov_sz_sorted(2)/fov_sz_sorted(1);
1733        if FOVrelativeDim~=sliceRelativeDim
1734          ind=find([size(kspace,1),size(kspace,2)]==ks_sz_sorted(1));
1735          data_sz = size(kspace);
1736          data_sz(ind) = round((fov_sz_sorted(1)/fov_sz_sorted(2))*ks_sz_sorted(2));
1737          DoZeroPadding=true;
1738        else
1739          data_sz = size(kspace);
1740          DoZeroPadding=false;
1741        end
1742      end
1743     
1744      %% Force zeropadding on. Force images to be square.
1745    elseif Dat.ZeroPadding==1
1746      ks_sz_max = max(size(kspace,1),size(kspace,2));
1747      data_sz = size(kspace);
1748      data_sz(1:2)=ks_sz_max;
1749      DoZeroPadding=true;
1750     
1751      %% Force zeropadding off
1752    else
1753      data_sz = size(kspace);
1754      DoZeroPadding=false;
1755    end
1756   
1757    %% Fourier transform 2D and EPI image data
1758    if Dat.AcqType==2 || Dat.isEPIdata
1759      sz=size(kspace,3);
1760     
1761      if DoZeroPadding
1762        data=zeros(data_sz,class(kspace));
1763        data(1:size(kspace,1),1:size(kspace,2),:,:)=kspace;
1764      else
1765        data=kspace;
1766      end
1767      if ~Dat.ReturnKSpace
1768        kspace = [];
1769      end
1770      if Dat.ShowWaitbar
1771        [wb_h,txh]=aedes_calc_wait('Fourier transforming image data...');
1772      end
1773     
1774      %data=permute(fftshift(fftshift(abs(fft(permute(fft(data),[2 1 3 4]))),1),2),[2 1 3 4]);
1775      if ~Dat.FastDataRead
1776        for tt=1:size(data,4)
1777          data(:,:,:,tt)=fftshift(fftshift(abs(fft(fft(data(:,:,:,tt),[],1),[],2)),1),2);
1778        end
1779      else
1780        data=fftshift(fftshift(abs(fft(fft(data,[],1),[],2)),1),2);
1781      end
1782    else
1783      %% Fourier transform 3D image data
1784      if Dat.ShowWaitbar
1785        [wb_h,txh]=aedes_calc_wait('Fourier transforming 3D image data...');
1786      end
1787     
1788      if DoZeroPadding
1789        % Check zeropadding in the 3rd dimension if zeropadding is set to
1790        % "auto"
1791        if Dat.ZeroPadding==2
1792          if isfield(procpar,'lpe2') && isfield(procpar,'lpe')
1793            data_sz(3) = max(1,round((procpar.lpe2/procpar.lpe)*data_sz(2)*Dat.ArrayLength));
1794          end
1795        end
1796        data=zeros(data_sz,class(kspace));
1797        data = abs(fftshift(fftn(kspace,data_sz)));
1798      else
1799        data=zeros(data_sz,class(kspace));
1800        data = abs(fftshift(fftn(kspace)));
1801      end
1802    end
1803   
1804    % Remove reference image if requested
1805    if Dat.isEPIdata && Dat.RemoveEPIphaseIm
1806      data = data(:,:,:,2:end);
1807    end
1808   
1809    % Reorient images and remove the reference image if requested
1810    if Dat.OrientImages && ~isempty(procpar) && ...
1811        isfield(procpar,'orient') && any(Dat.AcqType==[2 3 4])
1812      orient = procpar.orient{1};
1813      if any(strcmpi(orient,{'xyz','trans90','cor90','sag90'}))
1814        data = flipdim(aedes_rot3d(data,1,3),2);
1815      elseif strcmpi(orient,'oblique')
1816        data = flipdim(flipdim(data,1),2);
1817      else
1818        data = flipdim(flipdim(data,1),2);
1819      end
1820    end
1821   
1822   
1823   
1824% $$$   %% Fourier transform 4D image data
1825% $$$   elseif Dat.AcqType==4
1826% $$$     data = abs(fftshift(fftshift(fft2(kspace),2),1));
1827  end
1828 
1829  %% Delete waitbar
1830  if Dat.ShowWaitbar && ishandle(wb_h)
1831    delete(wb_h)
1832  end
1833end
1834
1835
1836
1837
1838
1839%========================================================
1840% A function for creating phasetables for fse and fsems
1841%========================================================
1842function phasetable=l_CreateFsemsPhaseTable(procpar)
1843
1844phasetable = [];
1845if ~isfield(procpar,'etl') || ~isfield(procpar,'kzero')
1846  return
1847end
1848etl = procpar.etl;
1849kzero = procpar.kzero;
1850nv = procpar.nv;
1851
1852t = (-nv/2+1):(nv/2);
1853t = t(:);
1854tt = flipud(t);
1855tt = tt(:);
1856
1857phasetable = [reshape(t(1:nv/2),[],etl);...
1858  flipud(reshape(tt(1:nv/2),[],etl))];
1859phasetable = circshift(fliplr(phasetable),[0 kzero-1]);
1860
1861
1862%% - EOF - %%
1863return
Note: See TracBrowser for help on using the repository browser.

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