source: aedes_readbruker.m @ 175

Last change on this file since 175 was 175, checked in by tjniskan, 8 years ago
  • Added scaling for Bruker data (reconstructed 2dseq -files).
  • Some progress in reconstructing Bruker raw data (FID-files).

M aedes_readbruker.m
M aedes_revision.m

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

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