source: aedes_readfid.m

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

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

  • Fixed overwriting ROIs in aedes_roi_copy_gui.m and added additional

information in the ROI copy button tooltips.

  • Some further iterations to the new readvnmr-function.

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

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

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