source: aedes_readbruker.m

Last change on this file was 201, checked in by tjniskan, 5 years ago
  • A little bug fix

M aedes_readbruker.m
M aedes_revision.m

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

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