source: aedes_readbruker.m @ 200

Last change on this file since 200 was 200, checked in by tjniskan, 7 years ago
  • Fixed a Windows related bug in aedes_export_gui.m
  • Fixed reading of 2dseq data with more that 4 dimensions

M aedes_readbruker.m
M aedes_export_gui.m
M aedes_revision.m

File size: 26.4 KB
Line 
1function DATA = aedes_readbruker(filename,varargin)
2% AEDES_READBRUKER - Read Bruker file formats (raw FID-files and
3%                    reconstructed 2DSEQ files)
4%   
5%
6% Synopsis:
7%        DATA=aedes_readbruker(filename,'PropertyName1',value1,'PropertyName2',value2,...)
8%        DATA=aedes_readbruker(filename,'header')
9%
10% Description:
11%        The function reads Bruker FID- and 2DSEQ files and returns a structure
12%        DATA with fields DataFormat, HDR, FTDATA, KSPACE, PROCPAR, and
13%        PHASETABLE. The fields of the DATA structure are constructed as
14%        follows:
15%       
16%        DATA
17%          |-> DataFormat         (string identifier for data format 'bruker_raw'/'bruker_reco')
18%          |-> HDR
19%                |-> FileHeader   (data file header)
20%                |-> BlockHeaders (Varian FID file data block headers)
21%                |-> fname        (file name)
22%                |-> fpath        (file path)
23%                |-> Param        (parameter values used by AEDES_READBRUKER to read the data)
24%          |-> FTDATA             (real fourier transformed image data)
25%          |-> KSPACE             (complex k-space, empty by default)
26%          |-> PROCPAR            (parameters from procpar file, NOTE: Varian FID files only)
27%          |-> PHASETABLE         (phasetable)
28%
29%        The DATA structure is returned as empty (without the above described fields)
30%        if an error has occurred during reading.
31%       
32%        The first input argument is either a string containing full path to
33%        the FID/2DSEQ-file or an empty string. Only the data file headers
34%        can be read by giving a string 'header' as the second input
35%        argument.
36%
37%        By default the k-space is not returned, i.e. the field KSPACE is
38%        empty. The returned data can be adjusted by using the 'return'
39%        property and values 1, 2, or 3 (see below for more information).
40%        The return option is ignored for 2DSEQ files.
41%
42%        The supported property-value pairs in AEDES_READBRUKER (property strings
43%        are not case sensitive):
44%
45%        Property:        Value:                Description:
46%        *********        ******                ************
47%        'Return'      :  [ {1} | 2 | 3 | 4 ]   % 1=return only ftdata (default)
48%                                               % 2=return only k-space
49%                                               % 3=return both ftdata & kspace
50%                                               % 4=return raw kspace
51%
52%        'wbar'        : [ {'on'} | 'off' ]     % Show/hide waitbar
53%
54%        'Precision'   : [{'double'}|'single']  % The precision for the
55%                                               % reconstructed and scaled
56%                                               % data. (default='double')
57%
58%        '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                        else
201                                error('Value for Force_4D parameter can be either ''on'' or ''off''.')
202                        end
203                otherwise
204                        error('Unknown parameter "%s".',param)
205        end
206end
207
208if isRawData
209        [hdr,msg] = l_ReadHeaderFid(filename);
210        if isempty(hdr)
211                error(msg);
212        end
213        [data,kspace,pe1_table,msg] = l_ReadDataFid(filename,hdr,Dat);
214        DATA.DataFormat = 'bruker_raw';
215else
216        [hdr,msg] = l_ReadHeader2dseq(filename);
217        if isempty(hdr)
218                error(msg);
219        end
220        [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat);
221        DATA.DataFormat = 'bruker_reco';
222        pe1_table = [];
223end
224if isempty(data) && isempty(kspace)
225        DATA = [];
226        error(msg);
227end
228
229DATA.HDR.FileHeader = hdr;
230DATA.HDR.BlockHeader = [];
231DATA.HDR.fname = [fn,fe];
232DATA.HDR.fpath = [fp,filesep];
233DATA.FTDATA = data;
234if Dat.return == 1
235        DATA.KSPACE = [];
236        clear kspace;
237else
238        DATA.KSPACE = kspace;
239end
240DATA.PROCPAR = [];
241DATA.PHASETABLE = pe1_table;
242
243
244% - Subfunctions -
245
246%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
247% Read header files for raw data
248%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
249function [hdr,msg] = l_ReadHeaderFid(filename)
250
251hdr = [];
252msg = '';
253
254% Parse filename
255[fp,fn,fe] = fileparts(filename);
256
257% Read Header files -------------------------
258
259% Read ACQP file
260acqp_file = [fp,filesep,'acqp'];
261if exist(acqp_file,'file')~=2
262        msg = sprintf('Cannot find ACQP file in "%s".',fp);
263        return
264end
265try
266        hdr.acqp = aedes_readjcamp(acqp_file);
267catch
268        hdr = [];
269        msg = sprintf('Error while reading "%s". The error was "%s".',...
270                acqp_file,lasterr);
271        return
272end
273
274% Read Method file
275method_file = [fp,filesep,'method'];
276if exist(method_file,'file')~=2
277        msg = sprintf('Cannot find METHOD file in "%s".',fp);
278        return
279end
280try
281        hdr.method = aedes_readjcamp(method_file);
282catch
283        hdr = [];
284        msg = sprintf('Error while reading "%s". The error was "%s".',...
285                method_file,lasterr);
286        return
287end
288
289
290
291%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
292% Read header files for reconstructed 2dseq data
293%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
294function [hdr,msg] = l_ReadHeader2dseq(filename)
295
296hdr = [];
297msg = '';
298
299% Parse filename
300[fp,fn,fe] = fileparts(filename);
301
302% Read Header files -------------------------
303
304% Read visu_pars file
305visu_pars_file = [fp,filesep,'visu_pars'];
306if exist(visu_pars_file,'file')==2
307        try
308                hdr.visu_pars = aedes_readjcamp(visu_pars_file);
309        catch
310                hdr = [];
311                msg = sprintf('Error while reading "%s". The error was "%s".',...
312                        visu_pars_file,lasterr);
313                return
314        end
315else
316        % Warn if visu_pars file was not found
317        warning('AEDES_READBRUKER:VisuParsNotFound',...
318                'VISU_PARS file was not not found in "%s". The 2dseq file might not be read correctly.',fp);
319end
320
321% Read d3proc file. This is however deprecated and the Visu parameters
322% should be used instead...
323d3proc_file = [fp,filesep,'d3proc'];
324if exist(d3proc_file,'file')==2
325        try
326                hdr.d3proc = aedes_readjcamp(d3proc_file);
327        catch
328                warning('Could not read D3PROC file "%s".',d3proc_file);
329        end
330end
331
332% Read reco file
333reco_file = [fp,filesep,'reco'];
334if exist(reco_file,'file')==2
335        try
336                hdr.reco = aedes_readjcamp(reco_file);
337        catch
338                hdr = [];
339                msg = sprintf('Error while reading "%s". The error was "%s".',...
340                        reco_file,lasterr);
341                return
342        end
343end
344
345% Read procs file
346procs_file = [fp,filesep,'procs'];
347if exist(procs_file,'file')==2
348        try
349                hdr.procs = aedes_readjcamp(procs_file);
350        catch
351                warning('Could not read PROCS file "%s".',procs_file);
352        end
353end
354
355% Read also acqp and method files if they exist
356ind = find(fp==filesep);
357if length(ind) >= 2
358        fp_fid = fp(1:ind(end-1));
359       
360        % Try to read ACQP file
361        acqp_file = [fp_fid,filesep,'acqp'];
362        if exist(acqp_file,'file')==2
363                hdr.acqp = aedes_readjcamp(acqp_file);
364        end
365       
366        % Try to read Method file
367        method_file = [fp_fid,filesep,'method'];
368        if exist(method_file,'file')==2
369                hdr.method = aedes_readjcamp(method_file);
370        end
371end
372
373% Check that either VISU_PARS, RECO or D3PROC file is found
374if ~isfield(hdr,'visu_pars') && ~isfield(hdr,'reco') && ...
375                ~isfield(hdr,'d3proc')
376        hdr = [];
377        msg = 'Could not read VISU_PARS, RECO or D3PROC files.';
378        return
379end
380
381% Check for Functional Imaging Tool files
382if exist([fp,filesep,'fun'],'dir') == 7
383       
384        % To be implemented...
385       
386end
387
388return
389
390%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
391% Read raw Bruker data (fid files)
392%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
393function [data,kspace,pe1_table,msg] = l_ReadDataFid(filename,hdr,Dat)
394
395data = [];
396kspace = [];
397msg = '';
398
399% Get word type
400if strcmpi(hdr.acqp.GO_raw_data_format,'GO_32BIT_SGN_INT')
401        precision = ['int32=>',Dat.precision];
402elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_SGN_INT')
403        precision = ['int16=>',Dat.precision];
404elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_FLOAT')
405        precision = ['single=>',Dat.precision];
406else
407        msg = sprintf('Unknown word type "%s".',hdr.acqp.GO_raw_data_format);
408        return
409end
410
411% Get byteorder
412if strcmpi(hdr.acqp.BYTORDA,'little')
413        byteorder = 'ieee-le';
414else
415        byteorder = 'ieee-be';
416end
417
418% Get block size
419if strcmpi(hdr.acqp.GO_block_size,'continuous')
420        fixed_blocksize = false;
421elseif strcmpi(hdr.acqp.GO_block_size,'Standard_KBlock_format')
422        fixed_blocksize = true;
423else
424        msg = sprintf('Unknown block size identifier "%s".',hdr.acqp.GO_block_size);
425        return
426end
427
428% Open fid file for reading
429fid = fopen(filename,'r',byteorder);
430if fid < 0
431        hdr.acqp = [];
432        hdr.method = [];
433        msg = sprintf('Could not open file "%s" for reading.',filename);
434        return
435end
436
437% Read fid file
438tmp = fread(fid,inf,precision);
439fclose(fid);
440kspace = complex(tmp(1:2:end),tmp(2:2:end));
441
442% Handle block size
443if strcmp(hdr.acqp.GO_block_size,'Standard_KBlock_Format')
444        blk_size = round(length(kspace)/1024)*1024;
445else
446        blk_size = hdr.acqp.ACQ_size(1)/2;
447end
448kspace = reshape(kspace,blk_size,[]);
449
450% Reconstruct k-space
451if Dat.return ~= 4
452        % Issue a warning that reconstruction of raw data is experimental
453        warning('AEDES_READBRUKER:RawDataReconstructionExperimental',...
454                'The reconstruction of raw bruker data is EXPERIMENTAL and may not work correctly.');
455       
456        [kspace,pe1_table,pe2_table,msg] = l_ReconstructKspace(kspace,hdr,Dat);
457else
458        [pe1_table,pe2_table] = l_GetPeTable(hdr);
459end
460if isempty(kspace)
461        kspace = [];
462        return
463end
464
465% Fourier transform k-space
466if any(Dat.return==[1 3])
467        [data,msg] = l_doFFT(kspace,hdr,Dat);
468        if isempty(data)
469                kspace = [];
470                data = [];
471                return
472        end
473else
474        data = [];
475end
476
477
478%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
479% Read reconstructed Bruker data (2dseq files)
480%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
481function [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat)
482
483data = [];
484kspace = [];
485msg = '';
486
487% Check which files to use
488isVisuPars = false;
489isReco = false;
490isD3proc = false;
491if isfield(hdr,'visu_pars')
492        isVisuPars = true;
493end
494if isfield(hdr,'reco')
495        isReco = true;
496end
497if isfield(hdr,'d3proc')
498        isD3proc = true;
499end
500
501% Get word type
502if isVisuPars
503        % Use visu_pars information first
504        if strcmpi(hdr.visu_pars.VisuCoreWordType,'_8BIT_UNSGN_INT')
505                precision = ['uint8=>',Dat.precision];
506        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_16BIT_SGN_INT')
507                precision = ['int16=>',Dat.precision];
508        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_SGN_INT')
509                precision = ['int32=>',Dat.precision];
510        elseif strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_FLOAT')
511                precision = ['single=>',Dat.precision];
512        else
513                msg = sprintf('Unknown word type "%s".',hdr.visu_pars.VisuCoreWordType);
514                return
515        end
516elseif isReco
517        % If visu_pars was not found, use reco information
518        if strcmpi(hdr.reco.RECO_wordtype,'_8BIT_UNSGN_INT')
519                precision = ['uint8=>',Dat.precision];
520        elseif strcmpi(hdr.reco.RECO_wordtype,'_16BIT_SGN_INT')
521                precision = ['int16=>',Dat.precision];
522        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_SGN_INT')
523                precision = ['int32=>',Dat.precision];
524        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
525                precision = ['single=>',Dat.precision];
526        else
527                msg = sprintf('Unknown word type "%s".',hdr.reco.RECO_wordtype);
528                return
529        end
530else
531        % No VISU_PARS or RECO information available. Fall back to the depricated
532        % D3PROC information and try to guess the word type from file size...
533        sz = [hdr.d3proc.IM_SIX,...
534                hdr.d3proc.IM_SIY,...
535                hdr.d3proc.IM_SIZ,...
536                hdr.d3proc.IM_SIT];
537        d=dir(filename);
538        nBytes = d.bytes;
539        bytesPerPixel = nBytes/prod(sz);
540        if bytesPerPixel == 1
541                precision = ['uint8=>',Dat.precision];
542        elseif bytesPerPixel == 2
543                precision = ['int16=>',Dat.precision];
544        else
545                precision = ['int32=>',Dat.precision];
546        end
547end
548
549% Get byteorder
550if isVisuPars
551        if strcmpi(hdr.visu_pars.VisuCoreByteOrder,'littleEndian')
552                byteorder = 'ieee-le';
553        else
554                byteorder = 'ieee-be';
555        end
556elseif isReco
557        if strcmpi(hdr.reco.RECO_byte_order,'littleEndian')
558                byteorder = 'ieee-le';
559        else
560                byteorder = 'ieee-be';
561        end
562else
563        % Use little endian if no information is available
564        byteorder = 'ieee-le';
565end
566
567% Get image type
568if isVisuPars
569        if isfield(hdr.visu_pars,'VisuCoreFrameType')
570                if strcmpi(hdr.visu_pars.VisuCoreFrameType,'COMPLEX_IMAGE')
571                        isDataComplex = true;
572                elseif strcmpi(hdr.visu_pars.VisuCoreFrameType,'MAGNITUDE_IMAGE')
573                        isDataComplex = false;
574                else
575                        msg = sprintf('Unknown image type "%s".',hdr.visu_pars.VisuCoreFrameType);
576                        return
577                end
578        else
579                % If VisuCoreFrameType parameter was not found in visu_pars, assume
580                % MAGNITUDE...
581                isDataComplex = false;
582        end
583elseif isReco
584        if strcmpi(hdr.reco.RECO_image_type,'COMPLEX_IMAGE')
585                isDataComplex = true;
586        elseif strcmpi(hdr.reco.RECO_image_type,'MAGNITUDE_IMAGE')
587                isDataComplex = false;
588        else
589                msg = sprintf('Unknown image type "%s".',hdr.reco.RECO_image_type);
590                return
591        end
592else
593        % Assume magnitude image if information is not available
594        isDataComplex = false;
595end
596
597% Open 2dseq file for reading
598fid = fopen(filename,'r',byteorder);
599if fid < 0
600        msg = sprintf('Could not open 2dseq file "%s" for reading.');
601        return
602end
603
604% Show waitbar
605if Dat.showWbar
606        wbh = aedes_calc_wait('Reading data from Bruker 2dseq file.');
607end
608
609% Check if word type is float
610isWordFloat = false;
611if isVisuPars && strcmpi(hdr.visu_pars.VisuCoreWordType,'_32BIT_FLOAT')
612        isWordFloat = true;
613elseif isReco && strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
614        isWordFloat = true;
615end
616
617% - Read data ---------------------------------
618[data,count] = fread(fid,inf,precision);
619fclose(fid);
620
621% Try to get image dimensions
622im_size = ones(1,4);
623if isVisuPars
624        for ii=1:length(hdr.visu_pars.VisuCoreSize)
625                im_size(ii) = hdr.visu_pars.VisuCoreSize(ii);
626        end
627        im_size(length(hdr.visu_pars.VisuCoreSize)+1) = hdr.visu_pars.VisuCoreFrameCount;
628else
629        % Fall back to the depricated d3proc information, if visu_pars
630        % information is not available...
631        im_size = [hdr.d3proc.IM_SIX,...
632                hdr.d3proc.IM_SIY,...
633                hdr.d3proc.IM_SIZ,...
634                hdr.d3proc.IM_SIT];
635end
636
637% Reshape to correct size
638data = reshape(data,im_size);
639
640% Map integer data to single or double values
641slope = [];
642offset = [];
643if isfield(hdr,'visu_pars') && ~isWordFloat
644        % The slope values in visu_pars are multiplicative inverses of the slope
645        % values in RECO file...
646        slope = hdr.visu_pars.VisuCoreDataSlope;
647        offset = hdr.visu_pars.VisuCoreDataOffs;
648        nDim = hdr.visu_pars.VisuCoreDim;
649elseif isfield(hdr,'reco') && ~isWordFloat
650        % Use RECO information if VISU_PARS is not available
651        slope = 1/hdr.reco.RECO_map_slope;
652        offset = hdr.reco.RECO_map_offset;
653        nDim = length(hdr.reco.RECO_size);
654elseif ~isWordFloat
655        warning('AEDES_READBRUKER:ScalingInformationNotFound',...
656                'Scaling information for 2dseq data was not found. Integer data is not scaled.');
657end
658
659% Calibrate data
660if ~isempty(slope) && ~isempty(offset)
661        for ii=1:length(slope)
662                if any(nDim == [1 2])
663                        data(:,:,ii) = data(:,:,ii)*slope(ii)+offset(ii);
664                else
665                        data(:,:,:,ii) = data(:,:,:,ii)*slope(ii)+offset(ii);
666                end
667        end
668end
669
670% Reshape frames using VisuFGOrderDesc information
671if isVisuPars && isfield(hdr.visu_pars,'VisuFGOrderDesc')
672        if hdr.visu_pars.VisuCoreFrameCount > 1 && ...
673                        size(hdr.visu_pars.VisuFGOrderDesc,1)>1
674                nTotalDims = hdr.visu_pars.VisuCoreDim+hdr.visu_pars.VisuFGOrderDescDim;
675                if size(hdr.visu_pars.VisuFGOrderDesc,1)==1
676                        frame_size = [hdr.visu_pars.VisuCoreSize,...
677                                hdr.visu_pars.VisuFGOrderDesc{1,1}];
678                else
679                        FG_size = [hdr.visu_pars.VisuFGOrderDesc{:,1}];
680                        frame_size = [hdr.visu_pars.VisuCoreSize,FG_size];
681                       
682                        % Force data to be 4D...
683                        if Dat.Force4D
684                                frame_size = [frame_size(1:3) prod(frame_size(4:end))];
685                        end
686                end
687                data = reshape(data,frame_size);
688        end
689end
690
691% Permute to correct orientation
692%data = permute(data,[2 1 3 4]);
693nDim = ndims(data);
694data = permute(data,[2 1 3:nDim]);
695
696if isDataComplex       
697        kspace = complex(data(:,:,:,1),data(:,:,:,2));
698        data = [];
699end
700
701% Close waitbar
702if Dat.showWbar
703        close(wbh);
704end
705
706%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
707% Reconstruct k-space
708%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
709function [kspace,pe1_table,pe2_table,msg] = l_ReconstructKspace(kspace,hdr,Dat)
710%
711% NOTE: The reconstruction of raw bruker data is EXPERIMENTAL and may work
712% correctly.
713%
714
715
716msg = '';
717
718% Don't try to reconstruct 1D data
719nDims = hdr.acqp.ACQ_dim;
720if nDims == 1
721        pe1_table = l_GetPeTable(hdr);
722        pe2_table = [];
723        return
724end
725
726% Handle EPI data
727if strncmpi(hdr.acqp.PULPROG,'EPI',3)
728        kspace = [];
729        msg = 'EPI reconstruction is not supported.';
730        return
731        %epi_matrix_size = hdr.method.PVM_Matrix;
732        %kspace = reshape(kspace,epi_matrix_size(1),epi_matrix_size(2),NSLICES,NR);
733end
734
735NI = hdr.acqp.NI; % Number of objects
736NSLICES = hdr.acqp.NSLICES; % Number of slices
737NR = hdr.acqp.NR; % Number of repetitions
738phase_factor = hdr.acqp.ACQ_phase_factor; % scans belonging to a single image
739im_size = hdr.acqp.ACQ_size;im_size(1)=im_size(1)/2;
740order = hdr.acqp.ACQ_obj_order;
741
742% Get number of receivers
743if isfield(hdr,'method') && isfield(hdr.method,'PVM_EncNReceivers')
744        nRcvrs = hdr.method.PVM_EncNReceivers;
745else
746        nRcvrs = 1;
747end
748
749% Get phase table
750usePeTable = true;
751[pe1_table,pe2_table] = l_GetPeTable(hdr);
752if isempty(pe1_table) && isempty(pe2_table)
753        usePeTable = false;
754end
755
756% Reshape data so that all echoes are in correct planes
757if nDims == 2
758        % Handle 2D data
759        kspace = reshape(kspace,im_size(1),...
760                nRcvrs,phase_factor,NI,im_size(2)/phase_factor,NR);
761        kspace = permute(kspace,[1 3 5 4 6 2]);
762        kspace = reshape(kspace,im_size(1),im_size(2),NI,NR,nRcvrs);
763elseif nDims == 3
764        % Handle 3D data
765        kspace = reshape(kspace,im_size(1),...
766                nRcvrs,phase_factor,im_size(2)/phase_factor,im_size(3),NR);
767        kspace = permute(kspace,[1 3 4 5 6 2]);
768        kspace = reshape(kspace,im_size(1),im_size(2),im_size(3),NR,nRcvrs);
769else
770        kspace = [];
771        msg = sprintf('The number of dimensions "%d" is not supported.',nDims);
772        return
773end
774
775% Sort echoes
776if usePeTable
777        if ~isempty(pe1_table)
778                kspace = kspace(:,pe1_table,:,:,:);
779        end
780        if nDims == 3 && ~isempty(pe2_table)
781                kspace = kspace(:,:,pe2_table,:,:);
782        end
783end
784
785% Sort object order
786if nDims == 2
787        kspace(:,:,order+1,:,:) = kspace;
788end
789
790% Permute to correct orientation
791kspace = flipdim(flipdim(kspace,1),2);
792
793%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
794% Check that the required reconstruction parameters are found in acqp
795%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
796function [ok,msg] = l_CheckAcqpParams(hdr)
797
798ok = false;
799msg = '';
800
801% Minimal list of ACQP parameters for data reconstruction
802minParamList = {...
803        'ACQ_dim',...
804        'ACQ_size',...
805        'NI',...
806        'NR',...
807        'BYTORDA',...
808        'ACQ_word_size',...
809        'ACQ_scan_size',...
810        'ACQ_scan_shift',...
811        'ACQ_phase_factor',...
812        'ACQ_rare_factor',...
813        'ACQ_phase_encoding_mode',...
814        'ACQ_phase_enc_start',...
815        'ACQ_spatial_size_1',...
816        'ACQ_spatial_size_2',...
817        'ACQ_spatial_phase_1',...
818        'ACQ_spatial_phase_2',...
819        'ACQ_obj_order',...
820        'ACQ_mod',...
821        'DSPFVS',...
822        'DSPFIRM',...
823        'DECIM'};
824
825param = fieldnames(hdr.acqp);
826
827
828
829%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
830% Fourier transform data
831%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
832function [data,msg] = l_doFFT(kspace,hdr,Dat)
833
834data = [];
835msg = '';
836
837% Check dimensions
838if hdr.acqp.ACQ_dim == 1
839  % 1D-image
840  AcqType = 1;
841elseif hdr.acqp.ACQ_dim == 2
842  % 2D-image
843  AcqType = 2;
844elseif hdr.acqp.ACQ_dim == 3
845  % 3D-image
846  AcqType = 3;
847else
848  AcqType = 4;
849end
850
851% Get number of receivers
852nRcvrs = size(kspace,5);
853
854% Do Zero filling if requested
855if any(AcqType == [2 3])
856        switch Dat.zerofill
857                case 'auto'
858                        % Zeropadding is on "auto", i.e. zeropad to FOV
859                        fov1 = hdr.acqp.ACQ_fov(1);
860                        fov2 = hdr.acqp.ACQ_fov(2);
861                        if AcqType == 3
862                                fov3 = hdr.acqp.ACQ_fov(3);
863                        end
864                        if AcqType==2
865                                % 2D data
866                                padSize = [hdr.acqp.ACQ_size(1)/2 ...
867                                        hdr.acqp.ACQ_size(1)/2*(fov2/fov1) ...
868                                        size(kspace,3)];
869                        elseif AcqType==3
870                                % 3D data
871                                padSize = [hdr.acqp.ACQ_size(1)/2 ...
872                                        hdr.acqp.ACQ_size(1)/2*(fov2/fov1) ...
873                                        hdr.acqp.ACQ_size(1)/2*(fov3/fov1)];
874                        end
875                case 'on'
876                        % Zeropad to square
877                        if AcqType==1
878                                padSize = hdr.acqp.ACQ_size(1)/2;
879                        elseif AcqType==2
880                                padSize = ones(1,2)*hdr.acqp.ACQ_size(1)/2;
881                                padSize(3) = size(kspace,3);
882                        else
883                                padSize = ones(1,3)*hdr.acqp.ACQ_size(1)/2;
884                        end
885                case 'off'
886                        padSize = [size(kspace,1) ...
887                                size(kspace,2) ...
888                                size(kspace,3)];
889                otherwise
890                        warning('Unknown zerofill option "%s". Skipping zerofilling...',...
891                                Dat.zerofill);
892                        padSize = [size(kspace,1) ...
893                                size(kspace,2) ...
894                                size(kspace,3)];
895        end
896       
897        % Pad with zeros
898        ks_sz = [size(kspace,1) ...
899                size(kspace,2) ...
900                size(kspace,3)];
901        if any(padSize>ks_sz)
902                kspace(padSize(1),padSize(2),padSize(3))=0;
903        end
904end
905
906% Allocate space for Fourier transformed data
907if nRcvrs>1 && Dat.SumOfSquares==1
908  data_sz = [padSize,size(kspace,4),size(kspace,5)];
909  data = zeros(data_sz,Dat.precision);
910else
911  data = zeros(size(kspace),Dat.precision);
912end
913
914
915if AcqType==1
916  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fft(kspace,[],1),1));
917elseif AcqType==2
918  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2));
919elseif AcqType==3
920  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fftshift(fftshift(fft(fft(fft(kspace,[],1),[],2),[],3),1),2),3));
921end
922
923% Calculate sum-of-squares image
924if nRcvrs>1 && Dat.SumOfSquares==1
925  % Calculate sum-of-squares
926  data = sqrt(sum(data(:,:,:,:,1:size(data,5)).^2,5));
927  data=abs(data);
928elseif nRcvrs>1 && Dat.SumOfSquares==2
929        % Move individual receiver data to 4th dimension if possible...
930        if size(data,4) == 1
931                data = reshape(data,size(data,1),size(data,2),size(data,3),size(data,5));
932        end
933end
934
935
936%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
937% Calculate phase table indexes
938%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
939function [pt1,pt2] = l_GetPeTable(hdr)
940
941if ~isfield(hdr,'acqp') || ~isfield(hdr,'method')
942        error('Cannot determine phase table, because ACQP and/or METHODS fields are missing.')
943end
944
945pt1 = [];
946pt2 = [];
947
948if isfield(hdr.acqp,'ACQ_spatial_phase_1')
949        [s,pt1] = sort(hdr.acqp.ACQ_spatial_phase_1);
950elseif isfield(hdr.method,'PVM_EncSteps1')
951        pt1 = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
952end
953
954if isfield(hdr.method,'PVM_EncSteps2')
955        pt2 = hdr.method.PVM_EncSteps2+(-min(hdr.method.PVM_EncSteps2(:))+1);
956end
957
958
959
960
961
962
963
964
965
966
Note: See TracBrowser for help on using the repository browser.

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