source: aedes_readbruker.m @ 196

Last change on this file since 196 was 196, checked in by tjniskan, 7 years ago
  • Fixed a small bug in aedes_readbruker

M aedes_readbruker.m
M aedes_revision.m

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

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