source: aedes_readbruker.m @ 176

Last change on this file since 176 was 176, checked in by tjniskan, 8 years ago
  • EPI is now read into a 4D matrix from Bruker 2dseq files.
  • Some modifications to Bruker raw data reading. Still far from ready though...

M aedes_readbruker.m
M aedes_revision.m

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

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