source: aedes_readbruker.m @ 199

Last change on this file since 199 was 199, checked in by tjniskan, 7 years ago
  • Fixed abug with 1D data in Brujer data support.

M aedes_readbruker.m
M aedes_revision.m

File size: 25.2 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 visu_pars file
280visu_pars_file = [fp,filesep,'visu_pars'];
281if exist(visu_pars_file,'file')==2
282        try
283                hdr.visu_pars = aedes_readjcamp(visu_pars_file);
284        catch
285                hdr = [];
286                msg = sprintf('Error while reading "%s". The error was "%s".',...
287                        visu_pars_file,lasterr);
288                return
289        end
290else
291        % Warn if visu_pars file was not found
292        warning('AEDES_READBRUKER:VisuParsNotFound',...
293                'VISU_PARS file was not not found in "%s". The 2dseq file might not be read correctly.',fp);
294end
295
296% Read d3proc file. This is however deprecated and the Visu parameters
297% should be used instead...
298d3proc_file = [fp,filesep,'d3proc'];
299if exist(d3proc_file,'file')==2
300        try
301                hdr.d3proc = aedes_readjcamp(d3proc_file);
302        catch
303                warning('Could not read D3PROC file "%s".',d3proc_file);
304        end
305end
306
307% Read reco file
308reco_file = [fp,filesep,'reco'];
309if exist(reco_file,'file')==2
310        try
311                hdr.reco = aedes_readjcamp(reco_file);
312        catch
313                hdr = [];
314                msg = sprintf('Error while reading "%s". The error was "%s".',...
315                        reco_file,lasterr);
316                return
317        end
318end
319
320% Read procs file
321procs_file = [fp,filesep,'procs'];
322if exist(procs_file,'file')==2
323        try
324                hdr.procs = aedes_readjcamp(procs_file);
325        catch
326                warning('Could not read PROCS file "%s".',procs_file);
327        end
328end
329
330% Read also acqp and method files if they exist
331ind = find(fp==filesep);
332if length(ind) >= 2
333        fp_fid = fp(1:ind(end-1));
334       
335        % Try to read ACQP file
336        acqp_file = [fp_fid,filesep,'acqp'];
337        if exist(acqp_file,'file')==2
338                hdr.acqp = aedes_readjcamp(acqp_file);
339        end
340       
341        % Try to read Method file
342        method_file = [fp_fid,filesep,'method'];
343        if exist(method_file,'file')==2
344                hdr.method = aedes_readjcamp(method_file);
345        end
346end
347
348% Check that either VISU_PARS, RECO or D3PROC file is found
349if ~isfield(hdr,'visu_pars') && ~isfield(hdr,'reco') && ...
350                ~isfield(hdr,'d3proc')
351        hdr = [];
352        msg = 'Could not read VISU_PARS, RECO or D3PROC files.';
353        return
354end
355
356% Check for Functional Imaging Tool files
357if exist([fp,filesep,'fun'],'dir') == 7
358       
359        % To be implemented...
360       
361end
362
363return
364
365%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
366% Read raw Bruker data (fid files)
367%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
368function [data,kspace,pe1_table,msg] = l_ReadDataFid(filename,hdr,Dat)
369
370data = [];
371kspace = [];
372msg = '';
373
374% Get word type
375if strcmpi(hdr.acqp.GO_raw_data_format,'GO_32BIT_SGN_INT')
376        precision = ['int32=>',Dat.precision];
377elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_SGN_INT')
378        precision = ['int16=>',Dat.precision];
379elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_FLOAT')
380        precision = ['single=>',Dat.precision];
381else
382        msg = sprintf('Unknown word type "%s".',hdr.acqp.GO_raw_data_format);
383        return
384end
385
386% Get byteorder
387if strcmpi(hdr.acqp.BYTORDA,'little')
388        byteorder = 'ieee-le';
389else
390        byteorder = 'ieee-be';
391end
392
393% Get block size
394if strcmpi(hdr.acqp.GO_block_size,'continuous')
395        fixed_blocksize = false;
396elseif strcmpi(hdr.acqp.GO_block_size,'Standard_KBlock_format')
397        fixed_blocksize = true;
398else
399        msg = sprintf('Unknown block size identifier "%s".',hdr.acqp.GO_block_size);
400        return
401end
402
403% Open fid file for reading
404fid = fopen(filename,'r',byteorder);
405if fid < 0
406        hdr.acqp = [];
407        hdr.method = [];
408        msg = sprintf('Could not open file "%s" for reading.',filename);
409        return
410end
411
412% Read fid file
413tmp = fread(fid,inf,precision);
414fclose(fid);
415kspace = complex(tmp(1:2:end),tmp(2:2:end));
416
417% Handle block size
418if strcmp(hdr.acqp.GO_block_size,'Standard_KBlock_Format')
419        blk_size = round(length(kspace)/1024)*1024;
420else
421        blk_size = hdr.acqp.ACQ_size(1)/2;
422end
423kspace = reshape(kspace,blk_size,[]);
424
425% Reconstruct k-space
426if Dat.return ~= 4
427        % Issue a warning that reconstruction of raw data is experimental
428        warning('AEDES_READBRUKER:RawDataReconstructionExperimental',...
429                'The reconstruction of raw bruker data is EXPERIMENTAL and may not work correctly.');
430       
431        [kspace,pe1_table,pe2_table,msg] = l_ReconstructKspace(kspace,hdr,Dat);
432else
433        [pe1_table,pe2_table] = l_GetPeTable(hdr);
434end
435if isempty(kspace)
436        kspace = [];
437        return
438end
439
440% Fourier transform k-space
441if any(Dat.return==[1 3])
442        [data,msg] = l_doFFT(kspace,hdr,Dat);
443        if isempty(data)
444                kspace = [];
445                data = [];
446                return
447        end
448else
449        data = [];
450end
451
452
453%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
454% Read reconstructed Bruker data (2dseq files)
455%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
456function [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat)
457
458data = [];
459kspace = [];
460msg = '';
461
462% Check which files to use
463isVisuPars = false;
464isReco = false;
465isD3proc = false;
466if isfield(hdr,'visu_pars')
467        isVisuPars = true;
468end
469if isfield(hdr,'reco')
470        isReco = true;
471end
472if isfield(hdr,'d3proc')
473        isD3proc = true;
474end
475
476% Get word type
477if isVisuPars
478        % Use visu_pars information first
479        if strcmpi(hdr.visu_pars.VisuCoreWordType,'_8BIT_UNSGN_INT')
480                precision = ['uint8=>',Dat.precision];
481        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_16BIT_SGN_INT')
482                precision = ['int16=>',Dat.precision];
483        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_SGN_INT')
484                precision = ['int32=>',Dat.precision];
485        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_FLOAT')
486                precision = ['single=>',Dat.precision];
487        else
488                msg = sprintf('Unknown word type "%s".',hdr.visu_pars.VisuCoreWordType);
489                return
490        end
491elseif isReco
492        % If visu_pars was not found, use reco information
493        if strcmpi(hdr.reco.RECO_wordtype,'_8BIT_UNSGN_INT')
494                precision = ['uint8=>',Dat.precision];
495        elseif strcmpi(hdr.reco.RECO_wordtype,'_16BIT_SGN_INT')
496                precision = ['int16=>',Dat.precision];
497        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_SGN_INT')
498                precision = ['int32=>',Dat.precision];
499        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
500                precision = ['single=>',Dat.precision];
501        else
502                msg = sprintf('Unknown word type "%s".',hdr.reco.RECO_wordtype);
503                return
504        end
505else
506        % No VISU_PARS or RECO information available. Fall back to the depricated
507        % D3PROC information and try to guess the word type from file size...
508        sz = [hdr.d3proc.IM_SIX,...
509                hdr.d3proc.IM_SIY,...
510                hdr.d3proc.IM_SIZ,...
511                hdr.d3proc.IM_SIT];
512        d=dir(filename);
513        nBytes = d.bytes;
514        bytesPerPixel = nBytes/prod(sz);
515        if bytesPerPixel == 1
516                precision = ['uint8=>',Dat.precision];
517        elseif bytesPerPixel == 2
518                precision = ['int16=>',Dat.precision];
519        else
520                precision = ['int32=>',Dat.precision];
521        end
522end
523
524% Get byteorder
525if isVisuPars
526        if strcmpi(hdr.visu_pars.VisuCoreByteOrder,'littleEndian')
527                byteorder = 'ieee-le';
528        else
529                byteorder = 'ieee-be';
530        end
531elseif isReco
532        if strcmpi(hdr.reco.RECO_byte_order,'littleEndian')
533                byteorder = 'ieee-le';
534        else
535                byteorder = 'ieee-be';
536        end
537else
538        % Use little endian if no information is available
539        byteorder = 'ieee-le';
540end
541
542% Get image type
543if isVisuPars
544        if isfield(hdr.visu_pars,'VisuCoreFrameType')
545                if strcmpi(hdr.visu_pars.VisuCoreFrameType,'COMPLEX_IMAGE')
546                        isDataComplex = true;
547                elseif strcmpi(hdr.visu_pars.VisuCoreFrameType,'MAGNITUDE_IMAGE')
548                        isDataComplex = false;
549                else
550                        msg = sprintf('Unknown image type "%s".',hdr.visu_pars.VisuCoreFrameType);
551                        return
552                end
553        else
554                % If VisuCoreFrameType parameter was not found in visu_pars, assume
555                % MAGNITUDE...
556                isDataComplex = false;
557        end
558elseif isReco
559        if strcmpi(hdr.reco.RECO_image_type,'COMPLEX_IMAGE')
560                isDataComplex = true;
561        elseif strcmpi(hdr.reco.RECO_image_type,'MAGNITUDE_IMAGE')
562                isDataComplex = false;
563        else
564                msg = sprintf('Unknown image type "%s".',hdr.reco.RECO_image_type);
565                return
566        end
567else
568        % Assume magnitude image if information is not available
569        isDataComplex = false;
570end
571
572% Open 2dseq file for reading
573fid = fopen(filename,'r',byteorder);
574if fid < 0
575        msg = sprintf('Could not open 2dseq file "%s" for reading.');
576        return
577end
578
579% Show waitbar
580if Dat.showWbar
581        wbh = aedes_calc_wait('Reading data from Bruker 2dseq file.');
582end
583
584% Check if word type is float
585isWordFloat = false;
586if isVisuPars && strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_FLOAT')
587        isWordFloat = true;
588elseif isReco && strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
589        isWordFloat = true;
590end
591
592% - Read data ---------------------------------
593[data,count] = fread(fid,inf,precision);
594fclose(fid);
595
596% Try to get image dimensions
597im_size = ones(1,4);
598if isVisuPars
599        for ii=1:length(hdr.visu_pars.VisuCoreSize)
600                im_size(ii) = hdr.visu_pars.VisuCoreSize(ii);
601        end
602        im_size(length(hdr.visu_pars.VisuCoreSize)+1) = hdr.visu_pars.VisuCoreFrameCount;
603else
604        % Fall back to the depricated d3proc information, if visu_pars
605        % information is not available...
606        im_size = [hdr.d3proc.IM_SIX,...
607                hdr.d3proc.IM_SIY,...
608                hdr.d3proc.IM_SIZ,...
609                hdr.d3proc.IM_SIT];
610end
611
612% Reshape to correct size
613data = reshape(data,im_size);
614
615% Map integer data to single or double values
616slope = [];
617offset = [];
618if isfield(hdr,'visu_pars') && ~isWordFloat
619        % The slope values in visu_pars are multiplicative inverses of the slope
620        % values in RECO file...
621        slope = hdr.visu_pars.VisuCoreDataSlope;
622        offset = hdr.visu_pars.VisuCoreDataOffs;
623        nDim = hdr.visu_pars.VisuCoreDim;
624elseif isfield(hdr,'reco') && ~isWordFloat
625        % Use RECO information if VISU_PARS is not available
626        slope = 1/hdr.reco.RECO_map_slope;
627        offset = hdr.reco.RECO_map_offset;
628        nDim = length(hdr.reco.RECO_size);
629elseif ~isWordFloat
630        warning('AEDES_READBRUKER:ScalingInformationNotFound',...
631                'Scaling information for 2dseq data was not found. Integer data is not scaled.');
632end
633
634% Calibrate data
635if ~isempty(slope) && ~isempty(offset)
636        for ii=1:length(slope)
637                if any(nDim == [1 2])
638                        data(:,:,ii) = data(:,:,ii)*slope(ii)+offset(ii);
639                else
640                        data(:,:,:,ii) = data(:,:,:,ii)*slope(ii)+offset(ii);
641                end
642        end
643end
644
645% Reshape frames using VisuFGOrderDesc information
646if isVisuPars && isfield(hdr.visu_pars,'VisuFGOrderDesc')
647        if hdr.visu_pars.VisuCoreFrameCount > 1 && ...
648                        size(hdr.visu_pars.VisuFGOrderDesc,1)>1
649                nTotalDims = hdr.visu_pars.VisuCoreDim+hdr.visu_pars.VisuFGOrderDescDim;
650                frame_size = [hdr.visu_pars.VisuCoreSize,...
651                        hdr.visu_pars.VisuFGOrderDesc{1,1},hdr.visu_pars.VisuFGOrderDesc{2,1}];
652                data = reshape(data,frame_size);
653        end
654end
655
656% Permute to correct orientation
657data = permute(data,[2 1 3 4]);
658
659if isDataComplex       
660        kspace = complex(data(:,:,:,1),data(:,:,:,2));
661        data = [];
662end
663
664% Close waitbar
665if Dat.showWbar
666        close(wbh);
667end
668
669%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
670% Reconstruct k-space
671%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
672function [kspace,pe1_table,pe2_table,msg] = l_ReconstructKspace(kspace,hdr,Dat)
673%
674% NOTE: The reconstruction of raw bruker data is EXPERIMENTAL and may work
675% correctly.
676%
677
678
679msg = '';
680
681% Don't try to reconstruct 1D data
682nDims = hdr.acqp.ACQ_dim;
683if nDims == 1
684        pe1_table = l_GetPeTable(hdr);
685        pe2_table = [];
686        return
687end
688
689% Handle EPI data
690if strncmpi(hdr.acqp.PULPROG,'EPI',3)
691        kspace = [];
692        msg = 'EPI reconstruction is not supported.';
693        return
694        %epi_matrix_size = hdr.method.PVM_Matrix;
695        %kspace = reshape(kspace,epi_matrix_size(1),epi_matrix_size(2),NSLICES,NR);
696end
697
698NI = hdr.acqp.NI; % Number of objects
699NSLICES = hdr.acqp.NSLICES; % Number of slices
700NR = hdr.acqp.NR; % Number of repetitions
701phase_factor = hdr.acqp.ACQ_phase_factor; % scans belonging to a single image
702im_size = hdr.acqp.ACQ_size;im_size(1)=im_size(1)/2;
703order = hdr.acqp.ACQ_obj_order;
704
705% Get number of receivers
706if isfield(hdr,'method') && isfield(hdr.method,'PVM_EncNReceivers')
707        nRcvrs = hdr.method.PVM_EncNReceivers;
708else
709        nRcvrs = 1;
710end
711
712% Get phase table
713usePeTable = true;
714[pe1_table,pe2_table] = l_GetPeTable(hdr);
715if isempty(pe1_table) && isempty(pe2_table)
716        usePeTable = false;
717end
718
719% Reshape data so that all echoes are in correct planes
720if nDims == 2
721        % Handle 2D data
722        kspace = reshape(kspace,im_size(1),...
723                nRcvrs,phase_factor,NI,im_size(2)/phase_factor,NR);
724        kspace = permute(kspace,[1 3 5 4 6 2]);
725        kspace = reshape(kspace,im_size(1),im_size(2),NI,NR,nRcvrs);
726elseif nDims == 3
727        % Handle 3D data
728        kspace = reshape(kspace,im_size(1),...
729                nRcvrs,phase_factor,im_size(2)/phase_factor,im_size(3),NR);
730        kspace = permute(kspace,[1 3 4 5 6 2]);
731        kspace = reshape(kspace,im_size(1),im_size(2),im_size(3),NR,nRcvrs);
732else
733        kspace = [];
734        msg = sprintf('The number of dimensions "%d" is not supported.',nDims);
735        return
736end
737
738% Sort echoes
739if usePeTable
740        if ~isempty(pe1_table)
741                kspace = kspace(:,pe1_table,:,:,:);
742        end
743        if nDims == 3 && ~isempty(pe2_table)
744                kspace = kspace(:,:,pe2_table,:,:);
745        end
746end
747
748% Sort object order
749if nDims == 2
750        kspace(:,:,order+1,:,:) = kspace;
751end
752
753% Permute to correct orientation
754kspace = flipdim(flipdim(kspace,1),2);
755
756%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
757% Check that the required reconstruction parameters are found in acqp
758%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
759function [ok,msg] = l_CheckAcqpParams(hdr)
760
761ok = false;
762msg = '';
763
764% Minimal list of ACQP parameters for data reconstruction
765minParamList = {...
766        'ACQ_dim',...
767        'ACQ_size',...
768        'NI',...
769        'NR',...
770        'BYTORDA',...
771        'ACQ_word_size',...
772        'ACQ_scan_size',...
773        'ACQ_scan_shift',...
774        'ACQ_phase_factor',...
775        'ACQ_rare_factor',...
776        'ACQ_phase_encoding_mode',...
777        'ACQ_phase_enc_start',...
778        'ACQ_spatial_size_1',...
779        'ACQ_spatial_size_2',...
780        'ACQ_spatial_phase_1',...
781        'ACQ_spatial_phase_2',...
782        'ACQ_obj_order',...
783        'ACQ_mod',...
784        'DSPFVS',...
785        'DSPFIRM',...
786        'DECIM'};
787
788param = fieldnames(hdr.acqp);
789
790
791
792%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
793% Fourier transform data
794%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
795function [data,msg] = l_doFFT(kspace,hdr,Dat)
796
797data = [];
798msg = '';
799
800% Check dimensions
801if hdr.acqp.ACQ_dim == 1
802  % 1D-image
803  AcqType = 1;
804elseif hdr.acqp.ACQ_dim == 2
805  % 2D-image
806  AcqType = 2;
807elseif hdr.acqp.ACQ_dim == 3
808  % 3D-image
809  AcqType = 3;
810else
811  AcqType = 4;
812end
813
814% Get number of receivers
815nRcvrs = size(kspace,5);
816
817% Do Zero filling if requested
818if any(AcqType == [2 3])
819        switch Dat.zerofill
820                case 'auto'
821                        % Zeropadding is on "auto", i.e. zeropad to FOV
822                        fov1 = hdr.acqp.ACQ_fov(1);
823                        fov2 = hdr.acqp.ACQ_fov(2);
824                        if AcqType == 3
825                                fov3 = hdr.acqp.ACQ_fov(3);
826                        end
827                        if AcqType==2
828                                % 2D data
829                                padSize = [hdr.acqp.ACQ_size(1)/2 ...
830                                        hdr.acqp.ACQ_size(1)/2*(fov2/fov1) ...
831                                        size(kspace,3)];
832                        elseif AcqType==3
833                                % 3D data
834                                padSize = [hdr.acqp.ACQ_size(1)/2 ...
835                                        hdr.acqp.ACQ_size(1)/2*(fov2/fov1) ...
836                                        hdr.acqp.ACQ_size(1)/2*(fov3/fov1)];
837                        end
838                case 'on'
839                        % Zeropad to square
840                        if AcqType==1
841                                padSize = hdr.acqp.ACQ_size(1)/2;
842                        elseif AcqType==2
843                                padSize = ones(1,2)*hdr.acqp.ACQ_size(1)/2;
844                                padSize(3) = size(kspace,3);
845                        else
846                                padSize = ones(1,3)*hdr.acqp.ACQ_size(1)/2;
847                        end
848                case 'off'
849                        padSize = [size(kspace,1) ...
850                                size(kspace,2) ...
851                                size(kspace,3)];
852                otherwise
853                        warning('Unknown zerofill option "%s". Skipping zerofilling...',...
854                                Dat.zerofill);
855                        padSize = [size(kspace,1) ...
856                                size(kspace,2) ...
857                                size(kspace,3)];
858        end
859       
860        % Pad with zeros
861        ks_sz = [size(kspace,1) ...
862                size(kspace,2) ...
863                size(kspace,3)];
864        if any(padSize>ks_sz)
865                kspace(padSize(1),padSize(2),padSize(3))=0;
866        end
867end
868
869% Allocate space for Fourier transformed data
870if nRcvrs>1 && Dat.SumOfSquares==1
871  data_sz = [padSize,size(kspace,4),size(kspace,5)];
872  data = zeros(data_sz,Dat.precision);
873else
874  data = zeros(size(kspace),Dat.precision);
875end
876
877
878if AcqType==1
879  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fft(kspace,[],1),1));
880elseif AcqType==2
881  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2));
882elseif AcqType==3
883  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fftshift(fftshift(fft(fft(fft(kspace,[],1),[],2),[],3),1),2),3));
884end
885
886% Calculate sum-of-squares image
887if nRcvrs>1 && Dat.SumOfSquares==1
888  % Calculate sum-of-squares
889  data = sqrt(sum(data(:,:,:,:,1:size(data,5)).^2,5));
890  data=abs(data);
891elseif nRcvrs>1 && Dat.SumOfSquares==2
892        % Move individual receiver data to 4th dimension if possible...
893        if size(data,4) == 1
894                data = reshape(data,size(data,1),size(data,2),size(data,3),size(data,5));
895        end
896end
897
898
899%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
900% Calculate phase table indexes
901%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
902function [pt1,pt2] = l_GetPeTable(hdr)
903
904if ~isfield(hdr,'acqp') || ~isfield(hdr,'method')
905        error('Cannot determine phase table, because ACQP and/or METHODS fields are missing.')
906end
907
908pt1 = [];
909pt2 = [];
910
911if isfield(hdr.acqp,'ACQ_spatial_phase_1')
912        [s,pt1] = sort(hdr.acqp.ACQ_spatial_phase_1);
913elseif isfield(hdr.method,'PVM_EncSteps1')
914        pt1 = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
915end
916
917if isfield(hdr.method,'PVM_EncSteps2')
918        pt2 = hdr.method.PVM_EncSteps2+(-min(hdr.method.PVM_EncSteps2(:))+1);
919end
920
921
922
923
924
925
926
927
928
929
Note: See TracBrowser for help on using the repository browser.

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