source: aedes_readbruker.m @ 194

Last change on this file since 194 was 194, checked in by tjniskan, 7 years ago
  • aedes_readjcamp now correctly reads matrices from JCAMP files into cell arrays.
  • aedes_readbruker now tries to read information from acqp and method files when opening a 2dseq file.
  • aedes_readbruker now reads also visu_pars file and uses that to read the 2dseq file if reco file is not found.

M aedes_readbruker.m
M aedes_readjcamp.m
M aedes_revision.m

File size: 19.4 KB
Line 
1function DATA = aedes_readbruker(filename,varargin)
2% AEDES_READBRUKER - Read Bruker file formats (raw FID-files and
3%                    reconstructed 2DSEQ files)
4%   
5%
6% Synopsis:
7%        DATA=aedes_readbruker(filename,'PropertyName1',value1,'PropertyName2',value2,...)
8%        DATA=aedes_readbruker(filename,'header')
9%
10% Description:
11%        The function reads Bruker FID- and 2DSEQ files and returns 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 'bruker_raw'/'bruker_reco')
18%          |-> HDR
19%                |-> FileHeader   (data file header)
20%                |-> BlockHeaders (Varian FID file data block headers)
21%                |-> fname        (file name)
22%                |-> fpath        (file path)
23%                |-> Param        (parameter values used by AEDES_READBRUKER 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, NOTE: Varian FID files only)
27%          |-> PHASETABLE         (phasetable)
28%
29%        The DATA structure is returned as empty (without the above described fields)
30%        if an error has occurred during reading.
31%       
32%        The first input argument is either a string containing full path to
33%        the FID/2DSEQ-file or an empty string. Only the data file headers
34%        can be read by giving a string 'header' as the second input
35%        argument.
36%
37%        By default the k-space is not returned, i.e. the field KSPACE is
38%        empty. The returned data can be adjusted by using the 'return'
39%        property and values 1, 2, or 3 (see below for more information).
40%        The return option is ignored for 2DSEQ files.
41%
42%        The supported property-value pairs in AEDES_READBRUKER (property strings
43%        are not case sensitive):
44%
45%        Property:        Value:                Description:
46%        *********        ******                ************
47%        'Return'      :  [ {1} | 2 | 3 | 4 ]   % 1=return only ftdata (default)
48%                                               % 2=return only k-space
49%                                               % 3=return both ftdata & kspace
50%                                               % 4=return raw kspace
51%
52%        'wbar'        : [ {'on'} | 'off' ]     % Show/hide waitbar
53%
54%        'Precision'   : [{'double'}|'single']  % The precision for the
55%                                               % reconstructed and scaled
56%                                               % data. (default='double')
57%
58%        NOTE: The support for reconstructing raw bruker data is
59%        EXPERIMENTAL and should not normally be used. If you don't need
60%        kspace information USE THE 2DSEQ FILES INSTEAD OF THE FID FILES!
61%
62% Examples:
63%        DATA=aedes_readbruker(filename)   % Read image data from 'filename'
64%
65% See also:
66%        AEDES, AEDES_READVNMR
67
68% This function is a part of Aedes - A graphical tool for analyzing
69% medical images
70%
71% Copyright (C) 2011 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uef.fi>
72%
73% Department of Applied Physics, Department of Neurobiology
74% University of Eastern Finland, Kuopio, FINLAND
75%
76% This program may be used under the terms of the GNU General Public
77% License version 2.0 as published by the Free Software Foundation
78% and appearing in the file LICENSE.TXT included in the packaging of
79% this program.
80%
81% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
82% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
83
84% Defaults
85DATA = [];
86isRawData = true;
87Dat.showWbar = true;
88Dat.precision = 'double';
89Dat.return = 1;
90Dat.zerofill = 'auto';
91
92if nargin == 0 || isempty(filename)
93        [fname,fpath,findex] = uigetfile({'fid;fid.*;2dseq','Bruker file formats (fid, 2dseq)';...
94                '*.*','All Files (*.*)'},'Select Bruker file formats');
95        if isequal(fname,0)
96                % Canceled
97                return
98        end
99        filename = [fpath,fname];
100end
101
102% Parse filename
103[fp,fn,fe] = fileparts(filename);
104if strcmpi(fn,'fid')
105        isRawData = true;
106        data_format = aedes_getdataformat(filename);
107        if strcmpi(data_format,'vnmr')
108                error('aedes_readbruker cannot read Varian FID-files. Use aedes_readvnmr instead.');
109        end
110elseif strcmpi(fn,'2dseq')
111        isRawData = false;
112else
113        error('Only Bruker FID and 2DSEQ files are supported.');
114end
115
116% Check varargin length
117if rem(length(varargin),2)==1
118        error('Input arguments may only include filename and param/value pairs.')
119end
120
121% Parse varargin
122for ii=1:2:length(varargin)
123        param = varargin{ii};
124        if ~ischar(param) || isempty(param)
125                error('Parameters must be strings.')
126        end
127        value = varargin{ii+1};
128        switch lower(param)
129                case 'wbar'
130                        Dat.wbar = value;
131                case 'precision'
132                        Dat.precision = value;
133                case 'return'
134                        if ~isnumeric(value) || ~isscalar(value) || isempty(value) || ...
135                                        isnan(value) || ~isreal(value) || (value-floor(value)) ~= 0 || ...
136                                        value < 1 || value > 4
137                                error('Invalid return value. Return value can be an integer in the range 1..4')
138                        end
139                        Dat.return = value;
140                otherwise
141                        error('Unknown parameter "%s".',param)
142        end
143end
144
145if isRawData
146        [hdr,msg] = l_ReadHeaderFid(filename);
147        if isempty(hdr)
148                error(msg);
149        end
150        [data,kspace,pe1_table,msg] = l_ReadDataFid(filename,hdr,Dat);
151        DATA.DataFormat = 'bruker_raw';
152else
153        [hdr,msg] = l_ReadHeader2dseq(filename);
154        if isempty(hdr)
155                error(msg);
156        end
157        [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat);
158        DATA.DataFormat = 'bruker_reco';
159        pe1_table = [];
160end
161if isempty(data) && isempty(kspace)
162        DATA = [];
163        error(msg);
164end
165
166DATA.HDR.FileHeader = hdr;
167DATA.HDR.BlockHeader = [];
168DATA.HDR.fname = [fn,fe];
169DATA.HDR.fpath = [fp,filesep];
170DATA.FTDATA = data;
171if Dat.return == 1
172        DATA.KSPACE = [];
173        clear kspace;
174else
175        DATA.KSPACE = kspace;
176end
177DATA.PROCPAR = [];
178DATA.PHASETABLE = pe1_table;
179
180
181% - Subfunctions -
182
183%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184% Read header files for raw data
185%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
186function [hdr,msg] = l_ReadHeaderFid(filename)
187
188hdr = [];
189msg = '';
190
191% Parse filename
192[fp,fn,fe] = fileparts(filename);
193
194% Read Header files -------------------------
195
196% Read ACQP file
197acqp_file = [fp,filesep,'acqp'];
198if exist(acqp_file,'file')~=2
199        msg = sprintf('Cannot find ACQP file in "%s".',fp);
200        return
201end
202try
203        hdr.acqp = aedes_readjcamp(acqp_file);
204catch
205        hdr = [];
206        msg = sprintf('Error while reading "%s". The error was "%s".',...
207                acqp_file,lasterr);
208        return
209end
210
211% Read Method file
212method_file = [fp,filesep,'method'];
213if exist(method_file,'file')~=2
214        msg = sprintf('Cannot find METHOD file in "%s".',fp);
215        return
216end
217try
218        hdr.method = aedes_readjcamp(method_file);
219catch
220        hdr = [];
221        msg = sprintf('Error while reading "%s". The error was "%s".',...
222                method_file,lasterr);
223        return
224end
225
226
227
228%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
229% Read header files for reconstructed 2dseq data
230%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
231function [hdr,msg] = l_ReadHeader2dseq(filename)
232
233hdr = [];
234msg = '';
235
236% Parse filename
237[fp,fn,fe] = fileparts(filename);
238
239% Read Header files -------------------------
240
241% Read d3proc file
242d3proc_file = [fp,filesep,'d3proc'];
243if exist(d3proc_file,'file')~=2
244        hdr = [];
245        msg = sprintf('Cannot find D3PROC file in "%s".',fp);
246        return
247end
248try
249        hdr.d3proc = aedes_readjcamp(d3proc_file);
250catch
251        hdr = [];
252        msg = sprintf('Error while reading "%s". The error was "%s".',...
253                d3proc_file,lasterr);
254        return
255end
256
257% Read procs file
258procs_file = [fp,filesep,'procs'];
259if exist(procs_file,'file')~=2
260        hdr = [];
261        msg = sprintf('Cannot find PROCS file in "%s".',fp);
262        return
263end
264try
265        hdr.procs = aedes_readjcamp(procs_file);
266catch
267        hdr = [];
268        msg = sprintf('Error while reading "%s". The error was "%s".',...
269                procs_file,lasterr);
270        return
271end
272
273% Read reco file
274reco_file = [fp,filesep,'reco'];
275if exist(reco_file,'file')==2
276        try
277                hdr.reco = aedes_readjcamp(reco_file);
278        catch
279                hdr = [];
280                msg = sprintf('Error while reading "%s". The error was "%s".',...
281                        reco_file,lasterr);
282                return
283        end
284end
285
286% Read visu_pars file
287visu_pars_file = [fp,filesep,'visu_pars'];
288if exist(visu_pars_file,'file')==2
289        try
290                hdr.visu_pars = aedes_readjcamp(visu_pars_file);
291        catch
292                hdr = [];
293                msg = sprintf('Error while reading "%s". The error was "%s".',...
294                        reco_file,lasterr);
295                return
296        end
297end
298
299% Read also acqp and method files if they exist
300ind = find(fp==filesep);
301if length(ind) >= 2
302        fp_fid = fp(1:ind(end-1));
303       
304        % Try to read ACQP file
305        acqp_file = [fp_fid,filesep,'acqp'];
306        if exist(acqp_file,'file')==2
307                hdr.acqp = aedes_readjcamp(acqp_file);
308        end
309       
310        % Try to read Method file
311        method_file = [fp_fid,filesep,'method'];
312        if exist(method_file,'file')==2
313                hdr.method = aedes_readjcamp(method_file);
314        end
315end
316
317% Check that either RECO or VISU_PARS file is found
318if exist(reco_file,'file')~=2 && exist(visu_pars_file,'file')~=2
319        warning('AEDES_READBRUKER:RecoVisuParsNotFound','Cannot find RECO and VISU_PARS files in "%s".',fp);
320end
321
322% Check for Functional Imaging Tool files
323if exist([fp,filesep,'fun'],'dir') == 7
324       
325        % To be implemented...
326       
327end
328
329return
330
331%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
332% Read raw Bruker data (fid files)
333%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
334function [data,kspace,pe1_table,msg] = l_ReadDataFid(filename,hdr,Dat)
335
336data = [];
337kspace = [];
338msg = '';
339
340% Get word type
341if strcmpi(hdr.acqp.GO_raw_data_format,'GO_32BIT_SGN_INT')
342        precision = ['int32=>',Dat.precision];
343elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_SGN_INT')
344        precision = ['int16=>',Dat.precision];
345elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_FLOAT')
346        precision = ['single=>',Dat.precision];
347else
348        msg = sprintf('Unknown word type "%s".',hdr.acqp.GO_raw_data_format);
349        return
350end
351
352% Get byteorder
353if strcmpi(hdr.acqp.BYTORDA,'little')
354        byteorder = 'ieee-le';
355else
356        byteorder = 'ieee-be';
357end
358
359% Get block size
360if strcmpi(hdr.acqp.GO_block_size,'continuous')
361        fixed_blocksize = false;
362elseif strcmpi(hdr.acqp.GO_block_size,'Standard_KBlock_format')
363        fixed_blocksize = true;
364else
365        msg = sprintf('Unknown block size identifier "%s".',hdr.acqp.GO_block_size);
366        return
367end
368
369% Open fid file for reading
370fid = fopen(filename,'r',byteorder);
371if fid < 0
372        hdr.acqp = [];
373        hdr.method = [];
374        msg = sprintf('Could not open file "%s" for reading.',filename);
375        return
376end
377
378% Read fid file
379tmp = fread(fid,inf,precision);
380fclose(fid);
381kspace = complex(tmp(1:2:end),tmp(2:2:end));
382
383% Handle block size
384if strcmp(hdr.acqp.GO_block_size,'Standard_KBlock_Format')
385        blk_size = round(length(kspace)/1024)*1024;
386else
387        blk_size = hdr.acqp.ACQ_size(1)/2;
388end
389kspace = reshape(kspace,blk_size,[]);
390
391% Reconstruct k-space
392if Dat.return ~= 4
393        % Issue a warning that reconstruction of raw data is experimental
394        warning('AEDES_READBRUKER:RawDataReconstructionExperimental',...
395                'The reconstruction of raw bruker data is EXPERIMENTAL and may work correctly.');
396       
397        [kspace,pe1_table,msg] = l_ReconstructKspace(kspace,hdr,Dat);
398else
399        pe1_table = l_GetPeTable(hdr);
400end
401if isempty(kspace)
402        kspace = [];
403        return
404end
405
406% Fourier transform k-space
407if any(Dat.return==[1 3])
408        [data,msg] = l_doFFT(kspace,hdr,Dat);
409        if isempty(data)
410                kspace = [];
411                data = [];
412                return
413        end
414else
415        data = [];
416end
417
418
419%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
420% Read reconstructed Bruker data (2dseq files)
421%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
422function [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat)
423
424data = [];
425kspace = [];
426msg = '';
427
428% Get word type
429if isfield(hdr,'reco')
430        if strcmpi(hdr.reco.RECO_wordtype,'_8BIT_UNSGN_INT')
431                precision = ['uint8=>',Dat.precision];
432        elseif strcmpi(hdr.reco.RECO_wordtype,'_16BIT_SGN_INT')
433                precision = ['int16=>',Dat.precision];
434        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_SGN_INT')
435                precision = ['int32=>',Dat.precision];
436        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
437                precision = ['single=>',Dat.precision];
438        else
439                msg = sprintf('Unknown word type "%s".',hdr.reco.RECO_wordtype);
440                return
441        end
442elseif isfield(hdr,'visu_pars')
443        if strcmpi(hdr.visu_pars.VisuCoreWordType,'_8BIT_UNSGN_INT')
444                precision = ['uint8=>',Dat.precision];
445        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_16BIT_SGN_INT')
446                precision = ['int16=>',Dat.precision];
447        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_SGN_INT')
448                precision = ['int32=>',Dat.precision];
449        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_FLOAT')
450                precision = ['single=>',Dat.precision];
451        else
452                msg = sprintf('Unknown word type "%s".',hdr.visu_pars.VisuCoreWordType);
453                return
454        end
455else
456        % Try to guess the word type from file size if reco file was not found...
457        sz = [hdr.d3proc.IM_SIX,...
458                hdr.d3proc.IM_SIY,...
459                hdr.d3proc.IM_SIZ,...
460                hdr.d3proc.IM_SIT];
461        d=dir(filename);
462        nBytes = d.bytes;
463        bytesPerPixel = nBytes/prod(sz);
464        if bytesPerPixel == 1
465                precision = ['uint8=>',Dat.precision];
466        elseif bytesPerPixel == 2
467                precision = ['int16=>',Dat.precision];
468        else
469                precision = ['int32=>',Dat.precision];
470        end
471end
472
473% Get byteorder
474if isfield(hdr,'reco')
475        if strcmpi(hdr.reco.RECO_byte_order,'littleEndian')
476                byteorder = 'ieee-le';
477        else
478                byteorder = 'ieee-be';
479        end
480elseif isfield(hdr,'visu_pars')
481        if strcmpi(hdr.visu_pars.VisuCoreByteOrder,'littleEndian')
482                byteorder = 'ieee-le';
483        else
484                byteorder = 'ieee-be';
485        end
486else
487        % Use little endian if no information is available
488        byteorder = 'ieee-le';
489end
490
491% Get image type
492if isfield(hdr,'reco')
493        if strcmpi(hdr.reco.RECO_image_type,'COMPLEX_IMAGE')
494                isDataComplex = true;
495        elseif strcmpi(hdr.reco.RECO_image_type,'MAGNITUDE_IMAGE')
496                isDataComplex = false;
497        else
498                msg = sprintf('Unknown image type "%s".',hdr.reco.RECO_image_type);
499                return
500        end
501else
502        % Assume magnitude image if information is not available
503        isDataComplex = false;
504end
505
506% Open 2dseq file for reading
507fid = fopen(filename,'r',byteorder);
508if fid < 0
509        msg = sprintf('Could not open 2dseq file "%s" for reading.');
510        return
511end
512
513% Show waitbar
514if Dat.showWbar
515        wbh = aedes_calc_wait('Reading data from Bruker 2dseq file.');
516end
517
518% Check if word type is float
519isWordFloat = false;
520if isfield(hdr,'reco') && strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
521        isWordFloat = true;
522elseif isfield(hdr,'visu_pars') && strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_FLOAT')
523        isWordFloat = true;
524end
525
526% Read data
527
528% - Read magnitude data ---------------------------------
529[data,count] = fread(fid,inf,precision);
530fclose(fid);
531
532% Reshape and permute to correct size and orientation
533data = reshape(data,[hdr.d3proc.IM_SIX,...
534        hdr.d3proc.IM_SIY,...
535        hdr.d3proc.IM_SIZ,...
536        hdr.d3proc.IM_SIT]);
537data = permute(data,[2 1 3 4]);
538
539% Map integer data to single or double values
540slope = [];
541offset = [];
542if isfield(hdr,'reco') && ~isWordFloat
543        slope = hdr.reco.RECO_map_slope;
544        offset = hdr.reco.RECO_map_offset;
545        nDim = length(hdr.reco.RECO_size);
546elseif isfield(hdr,'visu_pars') && ~isWordFloat
547        % Try to get the scaling information from visu_pars if reco file could
548        % not be found...
549        slope = hdr.visu_pars.VisuCoreDataSlope;
550        offset = hdr.visu_pars.VisuCoreDataOffs;
551        nDim = length(hdr.visu_pars.VisuCoreSize);
552elseif ~isWordFloat
553        warning('AEDES_READBRUKER:ScalingInformationNotFound',...
554                'Scaling information for 2dseq data was not found. Integer data is not scaled.');
555end
556
557if ~isempty(slope) && ~isempty(offset)
558        for ii=1:length(slope)
559                if any(nDim == [1 2])
560                        data(:,:,ii) = data(:,:,ii)/slope(ii)+offset(ii);
561                else
562                        data(:,:,:,ii) = data(:,:,:,ii)/slope(ii)+offset(ii);
563                end
564        end
565end
566
567if ~isDataComplex
568       
569        % The length of RECO_transposition should be same as NI.
570        if isfield(hdr,'reco')
571                NI = length(hdr.reco.RECO_transposition);
572                if NI~=hdr.d3proc.IM_SIZ
573                        data = reshape(data,hdr.d3proc.IM_SIY,...
574                                hdr.d3proc.IM_SIX,NI,[]);
575                end
576        end
577       
578else
579        data = reshape(data,hdr.reco.RECO_size(1),...
580                hdr.reco.RECO_size(2),...
581                hdr.reco.RECO_size(3),[]);
582       
583        kspace = complex(data(:,:,:,1),data(:,:,:,2));
584        data = [];
585end
586
587% Close waitbar
588if Dat.showWbar
589        close(wbh);
590end
591
592%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
593% Reconstruct k-space
594%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
595function [kspace,pe1_table,msg] = l_ReconstructKspace(kspace,hdr,Dat)
596
597msg = '';
598
599% Don't try to reconstruct 1D data
600nDims = hdr.acqp.ACQ_dim;
601if nDims == 1
602        pe1_table = l_GetPeTable(hdr);
603        return
604end
605NI = hdr.acqp.NI; % Number of objects
606NSLICES = hdr.acqp.NSLICES; % Number of slices
607NR = hdr.acqp.NR; % Number of repetitions
608phase_factor = hdr.acqp.ACQ_phase_factor; % scans belonging to a single image
609im_size = hdr.acqp.ACQ_size;im_size(1)=im_size(1)/2;
610order = hdr.acqp.ACQ_obj_order;
611
612% Get phase table
613usePeTable = true;
614pe1_table = l_GetPeTable(hdr);
615if isempty(pe1_table)
616        usePeTable = false;
617end
618
619% Reshape data so that all echoes are in correct planes
620kspace = reshape(kspace,im_size(1),...
621        phase_factor,NI,im_size(2)/phase_factor,NR);
622kspace = permute(kspace,[1 2 4 3 5]);
623kspace = reshape(kspace,im_size(1),im_size(2),NI,NR);
624
625% Handle EPI data
626if strncmpi(hdr.acqp.PULPROG,'EPI',3)
627        epi_matrix_size = hdr.method.PVM_Matrix;
628        kspace = reshape(kspace,epi_matrix_size(1),epi_matrix_size(2),NSLICES,NR);
629end
630
631% Sort echoes
632if usePeTable
633        kspace = kspace(:,pe1_table,:,:);
634end
635
636% Sort object order
637kspace(:,:,order+1,:) = kspace;
638
639% Permute to correct orientation
640kspace = flipdim(flipdim(kspace,1),2);
641
642%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
643% Check that the required reconstruction parameters are found in acqp
644%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
645function [ok,msg] = l_CheckAcqpParams(hdr)
646
647ok = false;
648msg = '';
649
650% Minimal list of ACQP parameters for data reconstruction
651minParamList = {...
652        'ACQ_dim',...
653        'ACQ_size',...
654        'NI',...
655        'NR',...
656        'BYTORDA',...
657        'ACQ_word_size',...
658        'ACQ_scan_size',...
659        'ACQ_scan_shift',...
660        'ACQ_phase_factor',...
661        'ACQ_rare_factor',...
662        'ACQ_phase_encoding_mode',...
663        'ACQ_phase_enc_start',...
664        'ACQ_spatial_size_1',...
665        'ACQ_spatial_size_2',...
666        'ACQ_spatial_phase_1',...
667        'ACQ_spatial_phase_2',...
668        'ACQ_obj_order',...
669        'ACQ_mod',...
670        'DSPFVS',...
671        'DSPFIRM',...
672        'DECIM'};
673
674param = fieldnames(hdr.acqp);
675
676
677
678%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
679% Fourier transform data
680%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
681function [data,msg] = l_doFFT(kspace,hdr,Dat)
682
683data = [];
684msg = '';
685
686% Check dimensions
687if hdr.acqp.ACQ_dim == 1
688  % 1D-image
689  AcqType = 1;
690elseif hdr.acqp.ACQ_dim == 2
691  % 2D-image
692  AcqType = 2;
693elseif hdr.acqp.ACQ_dim == 3
694  % 3D-image
695  AcqType = 3;
696else
697  AcqType = 4;
698end
699
700% Do Zero filling if requested
701if any(AcqType == [2 3])
702        switch Dat.zerofill
703                case 'auto'
704                       
705                case 'on'
706                       
707                case 'off'
708                       
709                otherwise
710                        warning('Unknown zerofill option "%s". Skipping zerofilling...',...
711                                Dat.zerofill);
712        end
713end
714
715if AcqType==1
716        %data = flipdim(abs(fftshift(fft(kspace,[],1),1)),1);
717        data = abs(fftshift(fft(flipdim(kspace,1),[],1),1));
718  %data = circshift(flipdim(abs(fftshift(fft(kspace,[],1),1)),1),[1 0 0]);
719elseif AcqType==2
720  data = abs(fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2));
721elseif AcqType==3
722  data = abs(fftshift(fftshift(fftshift(fft(fft(fft(kspace,[],1),[],2),[],3),1),2),3));
723end
724
725%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
726% Calculate phase table indexes
727%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
728function pt = l_GetPeTable(hdr)
729
730if ~isfield(hdr,'acqp') || ~isfield(hdr,'method')
731        error('Cannot determine phase table, because ACQP and/or METHODS fields are missing.')
732end
733
734if isfield(hdr.acqp,'ACQ_spatial_phase_1')
735        [s,pt] = sort(hdr.acqp.ACQ_spatial_phase_1);
736elseif isfield(hdr.method,'PVM_EncSteps1')
737        pt = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
738else
739        pt = [];
740end
741
742
743
744
745
746
747
748
749
750
751
Note: See TracBrowser for help on using the repository browser.

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