source: aedes_readbruker.m @ 174

Last change on this file since 174 was 174, checked in by tjniskan, 8 years ago
  • Added first alpha version for Bruker data support.

NOTE: Only reconstructed 2dseq files work at this time. This is the first alpha version so caution is adviced...

M aedes_getfilefilter.m
A aedes_readbruker.m
M aedes_data_read.m
M aedes_revision.m
M aedes_getdataformat.m

File size: 11.0 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,msg]=aedes_readbruker(filename,'PropertyName1',value1,'PropertyName2',value2,...)
8%        [DATA,msg]=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 (data block headers, not stored by default)
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 HDR, FTDATA,
30%        KSPACE, and PROCPAR fields) if an error has occurred during
31%        reading.
32%       
33%        The first input argument is either a string containing full path to
34%        the FID-file or the header structure HDR. Only the data file header
35%        can be read by giving a string 'header' as the second input
36%        argument.
37%
38%        By default the k-space is not returned, i.e. the field KSPACE is
39%        empty. The returned data can be adjusted by using the 'return'
40%        property and values 1, 2, or 3 (see below for more information).
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,msg]=aedes_readvnmr(filename)             % Read image data
56%        from 'filename'
57%
58% See also:
59%        AEDES_READBRUKERRAW, AEDES_READBRUKERRECO, AEDES
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.precision = 'double';
81
82if nargin == 0 || isempty(filename)
83        [fname,fpath,findex] = uigetfile({'fid;2dseq','Bruker file formats (fid, 2dseq)';...
84                '*.*','All Files (*.*)'},'Select Bruker file formats');
85        if isequal(fname,0)
86                % Canceled
87                return
88        end
89        filename = [fpath,fname];
90end
91
92% Parse filename
93[fp,fn,fe] = fileparts(filename);
94if strcmpi(fn,'fid')
95        isRawData = true;
96elseif strcmpi(fn,'2dseq')
97        isRawData = false;
98else
99        error('Only Bruker FID and 2DSEQ files are supported.');
100end
101
102% Check varargin length
103if rem(length(varargin),2)==1
104        error('Arguments may only include filename and param/value pairs.')
105end
106
107% Parse varargin
108for ii=1:2:length(varargin)
109       
110end
111
112if isRawData
113        [hdr,msg] = l_ReadHeaderFid(filename);
114        if isempty(hdr)
115                error(msg);
116        end
117        [data,kspace,msg] = l_ReadDataFid(filename,hdr,Dat);
118        DATA.DataFormat = 'bruker_raw';
119else
120        [hdr,msg] = l_ReadHeader2dseq(filename);
121        if isempty(hdr)
122                error(msg);
123        end
124        [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat);
125        DATA.DataFormat = 'bruker_reco';
126end
127if isempty(data)
128        DATA = [];
129        error(msg);
130end
131
132DATA.HDR.FileHeader = hdr;
133DATA.HDR.BlockHeader = [];
134DATA.HDR.fname = fn;
135DATA.HDR.fpath = [fp,filesep];
136DATA.FTDATA = data;
137DATA.KSPACE = kspace;
138DATA.PROCPAR = [];
139DATA.PHASETABLE = [];
140
141
142% - Subfunctions -
143
144%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145% Read header files for raw data
146%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147function [hdr,msg] = l_ReadHeaderFid(filename)
148
149hdr = [];
150msg = '';
151
152% Parse filename
153[fp,fn,fe] = fileparts(filename);
154
155% Read Header files -------------------------
156
157% Read ACQP file
158acqp_file = [fp,filesep,'acqp'];
159if exist(acqp_file,'file')~=2
160        msg = sprintf('Cannot find ACQP file in "%s".',fp);
161        return
162end
163try
164        hdr.acqp = aedes_readjcamp(acqp_file);
165catch
166        hdr = [];
167        msg = sprintf('Error while reading "%s". The error was "%s".',...
168                acqp_file,lasterr);
169        return
170end
171
172% Read Method file
173method_file = [fp,filesep,'method'];
174if exist(method_file,'file')~=2
175        msg = sprintf('Cannot find METHOD file in "%s".',fp);
176        return
177end
178try
179        hdr.method = aedes_readjcamp(method_file);
180catch
181        hdr = [];
182        msg = sprintf('Error while reading "%s". The error was "%s".',...
183                method_file,lasterr);
184        return
185end
186
187
188
189%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
190% Read header files for reconstructed 2dseq data
191%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
192function [hdr,msg] = l_ReadHeader2dseq(filename)
193
194hdr = [];
195msg = '';
196
197% Parse filename
198[fp,fn,fe] = fileparts(filename);
199
200% Read Header files -------------------------
201
202% Read d3proc file
203d3proc_file = [fp,filesep,'d3proc'];
204if exist(d3proc_file,'file')~=2
205        hdr = [];
206        msg = sprintf('Cannot find D3PROC file in "%s".',fp);
207        return
208end
209try
210        hdr.d3proc = aedes_readjcamp(d3proc_file);
211catch
212        hdr = [];
213        msg = sprintf('Error while reading "%s". The error was "%s".',...
214                d3proc_file,lasterr);
215        return
216end
217
218% Read procs file
219procs_file = [fp,filesep,'procs'];
220if exist(procs_file,'file')~=2
221        hdr = [];
222        msg = sprintf('Cannot find PROCS file in "%s".',fp);
223        return
224end
225try
226        hdr.procs = aedes_readjcamp(procs_file);
227catch
228        hdr = [];
229        msg = sprintf('Error while reading "%s". The error was "%s".',...
230                procs_file,lasterr);
231        return
232end
233
234% Read reco file
235reco_file = [fp,filesep,'reco'];
236if exist(reco_file,'file')~=2
237        hdr = [];
238        msg = sprintf('Cannot find RECO file in "%s".',fp);
239        return
240end
241try
242        hdr.reco = aedes_readjcamp(reco_file);
243catch
244        hdr = [];
245        msg = sprintf('Error while reading "%s". The error was "%s".',...
246                reco_file,lasterr);
247        return
248end
249
250
251
252return
253
254%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
255% Read raw Bruker data (fid files)
256%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
257function [data,kspace,msg] = l_ReadDataFid(filename,hdr,Dat)
258
259data = [];
260kspace = [];
261msg = '';
262
263% Get word type
264if strcmpi(hdr.acqp.GO_raw_data_format,'GO_32BIT_SGN_INT')
265        precision = ['int32=>',Dat.precision];
266elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_SGN_INT')
267        precision = ['int16=>',Dat.precision];
268elseif strcmpi(hdr.acqp.GO_raw_data_format,'GO_16BIT_FLOAT')
269        precision = ['single=>',Dat.precision];
270else
271        msg = sprintf('Unknown word type "%s".',hdr.acqp.GO_raw_data_format);
272        return
273end
274
275% Get byteorder
276if strcmpi(hdr.acqp.BYTORDA,'little')
277        byteorder = 'ieee-le';
278else
279        byteorder = 'ieee-be';
280end
281
282% Get block size
283if strcmpi(hdr.acqp.GO_block_size,'continuous')
284        fixed_blocksize = false;
285elseif strcmpi(hdr.acqp.GO_block_size,'Standard_KBlock_format')
286        fixed_blocksize = true;
287else
288        msg = sprintf('Unknown block size identifier "%s".',hdr.acqp.GO_block_size);
289        return
290end
291
292% Open fid file for reading
293fid = fopen(filename,'r',byteorder);
294if fid < 0
295        hdr.acqp = [];
296        hdr.method = [];
297        msg = sprintf('Could not open file "%s" for reading.',filename);
298        return
299end
300
301% Read fid file
302tmp = fread(fid,inf,precision);
303kspace = complex(tmp(1:2:end),tmp(2:2:end));
304kspace = reshape(kspace,[hdr.acqp.ACQ_size(1)/2 hdr.acqp.ACQ_size(2) hdr.acqp.NI]);
305
306% Reconstruct k-space
307[kspace,msg] = l_ReconstructKspace(kspace,hdr,Dat);
308if isempty(kspace)
309        kspace = [];
310        return
311end
312
313% Fourier transform k-space
314[data,msg] = l_doFFT(kspace,hdr,Dat);
315if isempty(data)
316        kspace = [];
317        data = [];
318        return
319end
320
321
322%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
323% Read reconstructed Bruker data (2dseq files)
324%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
325function [data,kspace,msg] = l_ReadData2dseq(filename,hdr,Dat)
326
327data = [];
328kspace = [];
329msg = '';
330
331% Get word type
332if strcmpi(hdr.reco.RECO_wordtype,'_8BIT_UNSGN_INT')
333        precision = ['uint8=>',Dat.precision];
334elseif strcmpi(hdr.reco.RECO_wordtype,'_16BIT_SGN_INT')
335        precision = ['int16=>',Dat.precision];
336elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_SGN_INT')
337        precision = ['int32=>',Dat.precision];
338elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
339        precision = ['single=>',Dat.precision];
340else
341        msg = sprintf('Unknown word type "%s".',hdr.reco.RECO_wordtype);
342        return
343end
344
345% Get byteorder
346if strcmpi(hdr.reco.RECO_byte_order,'littleEndian')
347        byteorder = 'ieee-le';
348else
349        byteorder = 'ieee-be';
350end
351
352% Get image type
353if strcmpi(hdr.reco.RECO_image_type,'COMPLEX_IMAGE')
354        isDataComplex = true;
355elseif strcmpi(hdr.reco.RECO_image_type,'MAGNITUDE_IMAGE')
356        isDataComplex = false;
357else
358        msg = sprintf('Unknown image type "%s".',hdr.reco.RECO_image_type);
359        return
360end
361
362% Open 2dseq file for reading
363fid = fopen(filename,'r',byteorder);
364if fid < 0
365        msg = sprintf('Could not open 2dseq file "%s" for reading.');
366        return
367end
368
369% Read data
370if ~isDataComplex
371        % - Read magnitude data ---------------------------------
372        [data,count] = fread(fid,inf,precision);
373        fclose(fid);
374       
375        % Reshape to correct size
376        data = reshape(data,[hdr.d3proc.IM_SIX,...
377                hdr.d3proc.IM_SIY,...
378                hdr.d3proc.IM_SIZ,...
379                hdr.d3proc.IM_SIT]);
380        data = permute(data,[2 1 3 4]);
381       
382else
383        % - Read complex data -----------------------------------
384        msg = 'Reading complex data from 2dseq has not been implemented yet.';
385        return
386end
387
388
389
390%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
391% Reconstruct k-space
392%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
393function [kspace,msg] = l_ReconstructKspace(kspace,hdr,Dat)
394
395msg = '';
396
397% Get phase table
398pe_table = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
399
400% Sort echoes
401kspace(:,pe_table,:,:) = kspace;
402
403%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
404% Fourier transform data
405%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
406function [data,msg] = l_doFFT(kspace,hdr,Dat)
407
408data = [];
409msg = '';
410
411% Check dimensions
412if hdr.acqp.ACQ_dim == 1
413  % 1D-image
414  AcqType = 1;
415elseif hdr.acqp.ACQ_dim == 2
416  % 2D-image
417  AcqType = 2;
418elseif hdr.acqp.ACQ_dim == 3
419  % 3D-image
420  AcqType = 3;
421else
422  AcqType = 4;
423end
424
425if AcqType==1
426  data = abs(fftshift(ifft(kspace,[],1),1));
427elseif AcqType==2
428  data = abs(fftshift(fftshift(ifft(ifft(kspace,[],1),[],2),1),2));
429elseif AcqType==3
430  data = abs(fftshift(fftshift(fftshift(ifft(ifft(ifft(kspace,[],1),[],2),[],3),1),2),3));
431end
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
Note: See TracBrowser for help on using the repository browser.

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