source: aedes_readbruker.m @ 177

Last change on this file since 177 was 177, checked in by tjniskan, 8 years ago
  • Added support for reading complex data from bruker 2dseq files.
  • Fixed scaling of 3D data from 2dseq files.
  • Aedes now prompts for handling of complex data (absolute, real

or imag parts) if only complex data is read from a file.

  • Added an option for EPI phase correction in aedes_readvnmr and epi_recon

M aedes_readbruker.m
M aedes_readvnmr.m
M aedes.m
M aedes_revision.m
M vnmr_recon/epi_recon.m

File size: 16.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%
59% Examples:
60%        DATA=aedes_readbruker(filename)   % Read image data from 'filename'
61%
62% See also:
63%        AEDES, AEDES_READVNMR
64
65% This function is a part of Aedes - A graphical tool for analyzing
66% medical images
67%
68% Copyright (C) 2011 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uef.fi>
69%
70% Department of Applied Physics, Department of Neurobiology
71% University of Eastern Finland, Kuopio, FINLAND
72%
73% This program may be used under the terms of the GNU General Public
74% License version 2.0 as published by the Free Software Foundation
75% and appearing in the file LICENSE.TXT included in the packaging of
76% this program.
77%
78% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
79% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
80
81% Defaults
82DATA = [];
83isRawData = true;
84Dat.showWbar = true;
85Dat.precision = 'double';
86Dat.return = 1;
87
88if nargin == 0 || isempty(filename)
89        [fname,fpath,findex] = uigetfile({'fid;2dseq','Bruker file formats (fid, 2dseq)';...
90                '*.*','All Files (*.*)'},'Select Bruker file formats');
91        if isequal(fname,0)
92                % Canceled
93                return
94        end
95        filename = [fpath,fname];
96end
97
98% Parse filename
99[fp,fn,fe] = fileparts(filename);
100if strcmpi(fn,'fid')
101        isRawData = true;
102        data_format = aedes_getdataformat(filename);
103        if strcmpi(data_format,'vnmr')
104                error('aedes_readbruker cannot read Varian FID-files. Use aedes_readvnmr instead.');
105        end
106elseif strcmpi(fn,'2dseq')
107        isRawData = false;
108else
109        error('Only Bruker FID and 2DSEQ files are supported.');
110end
111
112% Check varargin length
113if rem(length(varargin),2)==1
114        error('Input arguments may only include filename and param/value pairs.')
115end
116
117% Parse varargin
118for ii=1:2:length(varargin)
119        param = varargin{ii};
120        if ~ischar(param) || isempty(param)
121                error('Parameters must be strings.')
122        end
123        value = varargin{ii+1};
124        switch lower(param)
125                case 'wbar'
126                        Dat.wbar = value;
127                case 'precision'
128                        Dat.precision = value;
129                case 'return'
130                        if ~isnumeric(value) || ~isscalar(value) || isempty(value) || ...
131                                        isnan(value) || ~isreal(value) || (value-floor(value)) ~= 0 || ...
132                                        value < 1 || value > 4
133                                error('Invalid return value. Return value can be an integer in the range 1..4')
134                        end
135                        Dat.return = value;
136                otherwise
137                        error('Unknown parameter "%s".',param)
138        end
139end
140
141if isRawData
142        [hdr,msg] = l_ReadHeaderFid(filename);
143        if isempty(hdr)
144                error(msg);
145        end
146        [data,kspace,pe1_table,msg] = l_ReadDataFid(filename,hdr,Dat);
147        DATA.DataFormat = 'bruker_raw';
148else
149        [hdr,msg] = l_ReadHeader2dseq(filename);
150        if isempty(hdr)
151                error(msg);
152        end
153        [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat);
154        DATA.DataFormat = 'bruker_reco';
155        pe1_table = [];
156end
157if isempty(data) && isempty(kspace)
158        DATA = [];
159        error(msg);
160end
161
162DATA.HDR.FileHeader = hdr;
163DATA.HDR.BlockHeader = [];
164DATA.HDR.fname = fn;
165DATA.HDR.fpath = [fp,filesep];
166DATA.FTDATA = data;
167DATA.KSPACE = kspace;
168DATA.PROCPAR = [];
169DATA.PHASETABLE = pe1_table;
170
171
172% - Subfunctions -
173
174%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
175% Read header files for raw data
176%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
177function [hdr,msg] = l_ReadHeaderFid(filename)
178
179hdr = [];
180msg = '';
181
182% Parse filename
183[fp,fn,fe] = fileparts(filename);
184
185% Read Header files -------------------------
186
187% Read ACQP file
188acqp_file = [fp,filesep,'acqp'];
189if exist(acqp_file,'file')~=2
190        msg = sprintf('Cannot find ACQP file in "%s".',fp);
191        return
192end
193try
194        hdr.acqp = aedes_readjcamp(acqp_file);
195catch
196        hdr = [];
197        msg = sprintf('Error while reading "%s". The error was "%s".',...
198                acqp_file,lasterr);
199        return
200end
201
202% Read Method file
203method_file = [fp,filesep,'method'];
204if exist(method_file,'file')~=2
205        msg = sprintf('Cannot find METHOD file in "%s".',fp);
206        return
207end
208try
209        hdr.method = aedes_readjcamp(method_file);
210catch
211        hdr = [];
212        msg = sprintf('Error while reading "%s". The error was "%s".',...
213                method_file,lasterr);
214        return
215end
216
217
218
219%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
220% Read header files for reconstructed 2dseq data
221%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
222function [hdr,msg] = l_ReadHeader2dseq(filename)
223
224hdr = [];
225msg = '';
226
227% Parse filename
228[fp,fn,fe] = fileparts(filename);
229
230% Read Header files -------------------------
231
232% Read d3proc file
233d3proc_file = [fp,filesep,'d3proc'];
234if exist(d3proc_file,'file')~=2
235        hdr = [];
236        msg = sprintf('Cannot find D3PROC file in "%s".',fp);
237        return
238end
239try
240        hdr.d3proc = aedes_readjcamp(d3proc_file);
241catch
242        hdr = [];
243        msg = sprintf('Error while reading "%s". The error was "%s".',...
244                d3proc_file,lasterr);
245        return
246end
247
248% Read procs file
249procs_file = [fp,filesep,'procs'];
250if exist(procs_file,'file')~=2
251        hdr = [];
252        msg = sprintf('Cannot find PROCS file in "%s".',fp);
253        return
254end
255try
256        hdr.procs = aedes_readjcamp(procs_file);
257catch
258        hdr = [];
259        msg = sprintf('Error while reading "%s". The error was "%s".',...
260                procs_file,lasterr);
261        return
262end
263
264% Read reco file
265reco_file = [fp,filesep,'reco'];
266if exist(reco_file,'file')~=2
267        warning('Cannot find RECO file in "%s"',fp);
268else
269        try
270                hdr.reco = aedes_readjcamp(reco_file);
271        catch
272                hdr = [];
273                msg = sprintf('Error while reading "%s". The error was "%s".',...
274                        reco_file,lasterr);
275                return
276        end
277end
278
279% Check for Functional Imaging Tool files
280if exist([fp,filesep,'fun'],'dir') == 7
281       
282        % To be implemented...
283       
284end
285
286return
287
288%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
289% Read raw Bruker data (fid files)
290%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
291function [data,kspace,pe1_table,msg] = l_ReadDataFid(filename,hdr,Dat)
292
293data = [];
294kspace = [];
295msg = '';
296
297% Get word type
298if strcmpi(hdr.acqp.GO_raw_data_format,'GO_32BIT_SGN_INT')
299        precision = ['int32=>',Dat.precision];
300elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_SGN_INT')
301        precision = ['int16=>',Dat.precision];
302elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_FLOAT')
303        precision = ['single=>',Dat.precision];
304else
305        msg = sprintf('Unknown word type "%s".',hdr.acqp.GO_raw_data_format);
306        return
307end
308
309% Get byteorder
310if strcmpi(hdr.acqp.BYTORDA,'little')
311        byteorder = 'ieee-le';
312else
313        byteorder = 'ieee-be';
314end
315
316% Get block size
317if strcmpi(hdr.acqp.GO_block_size,'continuous')
318        fixed_blocksize = false;
319elseif strcmpi(hdr.acqp.GO_block_size,'Standard_KBlock_format')
320        fixed_blocksize = true;
321else
322        msg = sprintf('Unknown block size identifier "%s".',hdr.acqp.GO_block_size);
323        return
324end
325
326% Open fid file for reading
327fid = fopen(filename,'r',byteorder);
328if fid < 0
329        hdr.acqp = [];
330        hdr.method = [];
331        msg = sprintf('Could not open file "%s" for reading.',filename);
332        return
333end
334
335% Read fid file
336tmp = fread(fid,inf,precision);
337fclose(fid);
338kspace = complex(tmp(1:2:end),tmp(2:2:end));
339kspace = reshape(kspace,hdr.acqp.ACQ_size(1)/2,[]);
340
341% Reconstruct k-space
342[kspace,pe1_table,msg] = l_ReconstructKspace(kspace,hdr,Dat);
343if isempty(kspace)
344        kspace = [];
345        return
346end
347
348% Fourier transform k-space
349[data,msg] = l_doFFT(kspace,hdr,Dat);
350if isempty(data)
351        kspace = [];
352        data = [];
353        return
354end
355
356
357%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
358% Read reconstructed Bruker data (2dseq files)
359%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
360function [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat)
361
362data = [];
363kspace = [];
364msg = '';
365
366% Get word type
367if isfield(hdr,'reco')
368        if strcmpi(hdr.reco.RECO_wordtype,'_8BIT_UNSGN_INT')
369                precision = ['uint8=>',Dat.precision];
370        elseif strcmpi(hdr.reco.RECO_wordtype,'_16BIT_SGN_INT')
371                precision = ['int16=>',Dat.precision];
372        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_SGN_INT')
373                precision = ['int32=>',Dat.precision];
374        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
375                precision = ['single=>',Dat.precision];
376        else
377                msg = sprintf('Unknown word type "%s".',hdr.reco.RECO_wordtype);
378                return
379        end
380else
381        % Try to guess the word type from file size if reco file was not found...
382        sz = [hdr.d3proc.IM_SIX,...
383                hdr.d3proc.IM_SIY,...
384                hdr.d3proc.IM_SIZ,...
385                hdr.d3proc.IM_SIT];
386        d=dir(filename);
387        nBytes = d.bytes;
388        bytesPerPixel = nBytes/prod(sz);
389        if bytesPerPixel == 1
390                precision = ['uint8=>',Dat.precision];
391        elseif bytesPerPixel == 2
392                precision = ['int16=>',Dat.precision];
393        else
394                precision = ['int32=>',Dat.precision];
395        end
396end
397
398% Get byteorder
399if isfield(hdr,'reco')
400        if strcmpi(hdr.reco.RECO_byte_order,'littleEndian')
401                byteorder = 'ieee-le';
402        else
403                byteorder = 'ieee-be';
404        end
405else
406        % Use little endian if no information is available
407        byteorder = 'ieee-le';
408end
409
410% Get image type
411if isfield(hdr,'reco')
412        if strcmpi(hdr.reco.RECO_image_type,'COMPLEX_IMAGE')
413                isDataComplex = true;
414        elseif strcmpi(hdr.reco.RECO_image_type,'MAGNITUDE_IMAGE')
415                isDataComplex = false;
416        else
417                msg = sprintf('Unknown image type "%s".',hdr.reco.RECO_image_type);
418                return
419        end
420else
421        % Assume magnitude image if information is not available
422        isDataComplex = false;
423end
424
425% Open 2dseq file for reading
426fid = fopen(filename,'r',byteorder);
427if fid < 0
428        msg = sprintf('Could not open 2dseq file "%s" for reading.');
429        return
430end
431
432% Show waitbar
433if Dat.showWbar
434        wbh = aedes_calc_wait('Reading data from Bruker 2dseq file.');
435end
436
437% Read data
438if ~isDataComplex
439        % - Read magnitude data ---------------------------------
440        [data,count] = fread(fid,inf,precision);
441        fclose(fid);
442       
443        % Reshape and permute to correct size and orientation
444        data = reshape(data,[hdr.d3proc.IM_SIX,...
445                hdr.d3proc.IM_SIY,...
446                hdr.d3proc.IM_SIZ,...
447                hdr.d3proc.IM_SIT]);
448        data = permute(data,[2 1 3 4]);
449       
450        % Map integer data to single or double values
451        if ~strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
452                slope = hdr.reco.RECO_map_slope;
453                offset = hdr.reco.RECO_map_offset;
454                if length(hdr.reco.RECO_size) == 2
455                        nDim = 2;
456                elseif length(hdr.reco.RECO_size) == 3
457                        nDim = 3;
458                end
459                for ii=1:length(slope)
460                        if nDim == 2
461                                data(:,:,ii) = data(:,:,ii)/slope(ii)+offset(ii);
462                        else
463                                data(:,:,:,ii) = data(:,:,:,ii)/slope(ii)+offset(ii);
464                        end
465                end
466        end
467       
468        % The length of RECO_transposition should be same as NI.
469        NI = length(hdr.reco.RECO_transposition);
470        if NI~=hdr.d3proc.IM_SIZ
471                data = reshape(data,hdr.d3proc.IM_SIY,...
472                        hdr.d3proc.IM_SIX,NI,[]);
473        end
474       
475else
476        % - Read complex data -----------------------------------
477        [data,count] = fread(fid,inf,precision);
478        fclose(fid);
479       
480        % Reshape and permute to correct size and orientation
481        data = reshape(data,[hdr.d3proc.IM_SIX,...
482                hdr.d3proc.IM_SIY,...
483                hdr.d3proc.IM_SIZ,...
484                hdr.d3proc.IM_SIT]);
485        data = permute(data,[2 1 3 4]);
486       
487        % Map integer data to single or double values
488        if ~strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
489                slope = hdr.reco.RECO_map_slope;
490                offset = hdr.reco.RECO_map_offset;
491                if length(hdr.reco.RECO_size) == 2
492                        nDim = 2;
493                elseif length(hdr.reco.RECO_size) == 3
494                        nDim = 3;
495                end
496                for ii=1:length(slope)
497                        if nDim == 2
498                                data(:,:,ii) = data(:,:,ii)/slope(ii)+offset(ii);
499                        else
500                                data(:,:,:,ii) = data(:,:,:,ii)/slope(ii)+offset(ii);
501                        end
502                end
503        end
504       
505        data = reshape(data,hdr.reco.RECO_size(1),...
506                hdr.reco.RECO_size(2),...
507                hdr.reco.RECO_size(3),[]);
508       
509        kspace = complex(data(:,:,:,1),data(:,:,:,2));
510        data = [];
511end
512
513% Close waitbar
514if Dat.showWbar
515        close(wbh);
516end
517
518%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
519% Reconstruct k-space
520%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
521function [kspace,pe1_table,msg] = l_ReconstructKspace(kspace,hdr,Dat)
522
523msg = '';
524
525
526NI = hdr.acqp.NI; % Number of objects
527NSLICES = hdr.acqp.NSLICES; % Number of slices
528NR = hdr.acqp.NR; % Number of repetitions
529phase_factor = hdr.acqp.ACQ_phase_factor; % scans belonging to a single image
530nDims = hdr.acqp.ACQ_dim;
531im_size = hdr.acqp.ACQ_size;im_size(1)=im_size(1)/2;
532order = hdr.acqp.ACQ_obj_order;
533
534% Get phase table
535usePeTable = true;
536if isfield(hdr.acqp,'ACQ_spatial_phase_1')
537        [s,pe1_table] = sort(hdr.acqp.ACQ_spatial_phase_1);
538elseif isfield(hdr.method,'PVM_EncSteps1')
539        pe1_table = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
540else
541        pe1_table = [];
542        usePeTable = false;
543end
544
545% Reshape data so that all echoes are in correct planes
546kspace = reshape(kspace,im_size(1),...
547        phase_factor,NI,im_size(2)/phase_factor,NR);
548kspace = permute(kspace,[1 2 4 3 5]);
549kspace = reshape(kspace,im_size(1),im_size(2),NI,NR);
550
551% Handle EPI data
552if strncmpi(hdr.acqp.PULPROG,'EPI',3)
553        epi_matrix_size = hdr.method.PVM_Matrix;
554        kspace = reshape(kspace,epi_matrix_size(1),epi_matrix_size(2),NSLICES,NR);
555end
556
557% Sort echoes
558if usePeTable
559        kspace = kspace(:,pe1_table,:,:);
560end
561
562% Sort object order
563kspace(:,:,order+1,:) = kspace;
564
565% Permute to correct orientation
566kspace = flipdim(flipdim(kspace,1),2);
567
568%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
569% Check that the required reconstruction parameters are found in acqp
570%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
571function [ok,msg] = l_CheckAcqpParams(hdr)
572
573ok = false;
574msg = '';
575
576% Minimal list of ACQP parameters for data reconstruction
577minParamList = {...
578        'ACQ_dim',...
579        'ACQ_size',...
580        'NI',...
581        'NR',...
582        'BYTORDA',...
583        'ACQ_word_size',...
584        'ACQ_scan_size',...
585        'ACQ_scan_shift',...
586        'ACQ_phase_factor',...
587        'ACQ_rare_factor',...
588        'ACQ_phase_encoding_mode',...
589        'ACQ_phase_enc_start',...
590        'ACQ_spatial_size_1',...
591        'ACQ_spatial_size_2',...
592        'ACQ_spatial_phase_1',...
593        'ACQ_spatial_phase_2',...
594        'ACQ_obj_order',...
595        'ACQ_mod',...
596        'DSPFVS',...
597        'DSPFIRM',...
598        'DECIM'};
599
600param = fieldnames(hdr.acqp);
601
602
603
604%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
605% Fourier transform data
606%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
607function [data,msg] = l_doFFT(kspace,hdr,Dat)
608
609data = [];
610msg = '';
611
612% Check dimensions
613if hdr.acqp.ACQ_dim == 1
614  % 1D-image
615  AcqType = 1;
616elseif hdr.acqp.ACQ_dim == 2
617  % 2D-image
618  AcqType = 2;
619elseif hdr.acqp.ACQ_dim == 3
620  % 3D-image
621  AcqType = 3;
622else
623  AcqType = 4;
624end
625
626if AcqType==1
627  data = abs(fftshift(fft(kspace,[],1),1));
628elseif AcqType==2
629  data = abs(fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2));
630elseif AcqType==3
631  data = abs(fftshift(fftshift(fftshift(fft(fft(fft(kspace,[],1),[],2),[],3),1),2),3));
632end
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
Note: See TracBrowser for help on using the repository browser.

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