1 | function [DATA,msg_out] = aedes_readvnmr(filename,varargin) |
---|
2 | % AEDES_READVNMR - Read VNMR (Varian) FID-files |
---|
3 | % |
---|
4 | % |
---|
5 | % Synopsis: |
---|
6 | % [DATA,msg]=aedes_readvnmr(filename,'PropertyName1',value1,'PropertyName2',value2,...) |
---|
7 | % [DATA,msg]=aedes_readvnmr(filename,'header') |
---|
8 | % [DATA,msg]=aedes_readvnmr(filename,output_filename) |
---|
9 | % |
---|
10 | % Description: |
---|
11 | % The function reads VNMR (Varian) FID-file and return 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 'vnmr') |
---|
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_READFID 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) |
---|
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. The error message is returned in the second output |
---|
32 | % argument msg. If AEDES_READVNMR is called with only one output argument |
---|
33 | % (i.e. without MSG), the error message is echoed to the workspace. |
---|
34 | % |
---|
35 | % The first input argument is either a string containing full path to |
---|
36 | % the FID-file or the header structure HDR. Only the data file header |
---|
37 | % can be read by giving a string 'header' as the second input |
---|
38 | % argument. |
---|
39 | % |
---|
40 | % By default the k-space is not returned, i.e. the field KSPACE is |
---|
41 | % empty. The returned data can be adjusted by using the 'return' |
---|
42 | % property and values 1, 2, or 3 (see below for more information). |
---|
43 | % |
---|
44 | % The supported property-value pairs in AEDES_READVNMR (property strings |
---|
45 | % are not case sensitive): |
---|
46 | % |
---|
47 | % Property: Value: Description: |
---|
48 | % ********* ****** ************ |
---|
49 | % 'Return' : [ {1} | 2 | 3 | 4 ] % 1=return only ftdata (default) |
---|
50 | % % 2=return only k-space |
---|
51 | % % 3=return both ftdata & kspace |
---|
52 | % % 4=return raw kspace |
---|
53 | % |
---|
54 | % 'FastRead' : [{'off'} | 'on' ] % Enables an experimental |
---|
55 | % % method for very fast reading |
---|
56 | % % of fid-files. This is not as |
---|
57 | % % robust as the older |
---|
58 | % % method and can consume a lot |
---|
59 | % % of memory. |
---|
60 | % |
---|
61 | % 'SumOfSquares' : [ {1} | 2 | 3 ] % 1=Return only the |
---|
62 | % % sum-of-squares image |
---|
63 | % % for multireceiver |
---|
64 | % % data (default). |
---|
65 | % % 2=Return both SoQ and |
---|
66 | % % individual receiver data |
---|
67 | % % 3=Return only individual |
---|
68 | % % receiver data |
---|
69 | % % NOTE: This property has |
---|
70 | % % no effect on single |
---|
71 | % % receiver data. |
---|
72 | % |
---|
73 | % 'FileChunkSize' : [ megabytes ] % Chunk size in megabytes for |
---|
74 | % % reading fid-files |
---|
75 | % % (default=500 MB). |
---|
76 | % % Lowering the chunk size |
---|
77 | % % helps to conserve memory |
---|
78 | % % when reading large files, |
---|
79 | % % but is slower. This |
---|
80 | % % property has no effect if |
---|
81 | % % FastRead='off'. |
---|
82 | % |
---|
83 | % |
---|
84 | % 'OutputFile' : filename % Writes the slices into |
---|
85 | % % NIfTI-format files |
---|
86 | % % using the given filename. |
---|
87 | % |
---|
88 | % 'DCcorrection': [ {'off'} | 'on' ] % 'on'=use DC correction |
---|
89 | % % 'off'=don't use DC correction |
---|
90 | % % (default) |
---|
91 | % |
---|
92 | % 'Procpar' : (procpar-structure) % Input procpar |
---|
93 | % % structure. If this |
---|
94 | % % property is not set the |
---|
95 | % % procpar structure |
---|
96 | % % will be read using |
---|
97 | % % AEDES_READPROCPAR |
---|
98 | % |
---|
99 | % 'SortPSS' : [ 'off' | {'on'} ] % Sort slices in 2D stacks |
---|
100 | % % using pss. Turn this 'off' |
---|
101 | % % if interleaved slice |
---|
102 | % % order is to be preserved. |
---|
103 | % % The original pss is saved |
---|
104 | % % in procpar.pss_orig |
---|
105 | % |
---|
106 | % 'ZeroPadding' : ['off'|'on'|{'auto'}] % 'off' = force |
---|
107 | % % zeropadding off |
---|
108 | % % 'on' = force |
---|
109 | % % zeropadding on (force |
---|
110 | % % images to be square) |
---|
111 | % % 'auto' = autoselect |
---|
112 | % % using relative FOV |
---|
113 | % % dimensions (default) |
---|
114 | % |
---|
115 | % 'PhaseTable' : (custom phase table) % User-specified |
---|
116 | % % line order for |
---|
117 | % % k-space. The phase table |
---|
118 | % % is obtained from the file |
---|
119 | % % specified in |
---|
120 | % % procpar.petable if |
---|
121 | % % necessary. |
---|
122 | % |
---|
123 | % 'sorting' : [ 'off' | {'on'} ] % 'off' =disable k-space |
---|
124 | % % sorting, i.e. ignore the |
---|
125 | % % phase table and/or arrays. |
---|
126 | % % 'on' =sort k-space using |
---|
127 | % % phase table and/or array |
---|
128 | % % information if necessary |
---|
129 | % % (default) |
---|
130 | % |
---|
131 | % 'wbar' : [ {'on'} | 'off' ] % 'on' = show waitbar |
---|
132 | % % 'off' = don't show waitbar |
---|
133 | % |
---|
134 | % 'Precision' : [{'single'}|'double'] % Precision of the |
---|
135 | % % outputted data. |
---|
136 | % % 'single' = 32-bit float |
---|
137 | % % 'double' = 64-bit float |
---|
138 | % |
---|
139 | % 'OrientImages': [ {'on'} | 'off' ] % Orient FTDATA after |
---|
140 | % % reading the data using |
---|
141 | % % PROCPAR.orient property. |
---|
142 | % |
---|
143 | % 'RemoveEPIphaseIm' : [{'on'}|'off'] % Remove the phase image |
---|
144 | % % from EPI data. This |
---|
145 | % % option only affects EPI |
---|
146 | % % images. |
---|
147 | % |
---|
148 | % 'EPIBlockSize' : [integer value] % Block size (number of |
---|
149 | % % volumes) used for |
---|
150 | % % processing multireceiver |
---|
151 | % % EPI data. Default=300 |
---|
152 | % |
---|
153 | % 'EPIPhasedArrayData' : ['on'|{'off'}] % Return data from |
---|
154 | % % individual receivers from |
---|
155 | % % phased array EPI files. |
---|
156 | % |
---|
157 | % 'EPI_PC' : [{'auto'}|'pointwise'| % Phase correction for EPI. |
---|
158 | % 'triple'|'off'] % 'auto' = Choose correction |
---|
159 | % % based on procpar. |
---|
160 | % % 'pointwise' = Pointwise |
---|
161 | % % single reference. |
---|
162 | % % 'triple' = Use triple |
---|
163 | % % reference correction |
---|
164 | % % 'off' = Do not perform |
---|
165 | % % phase correction. |
---|
166 | % |
---|
167 | % |
---|
168 | % Examples: |
---|
169 | % [DATA,msg]=aedes_readvnmr(filename) % Read image data from 'filename' |
---|
170 | % [DATA,msg]=aedes_readvnmr(filename,'header') % Read only data file header |
---|
171 | % |
---|
172 | % [DATA,msg]=aedes_readvnmr(filename,'return',1) % Return only image data (default) |
---|
173 | % [DATA,msg]=aedes_readvnmr(filename,'return',2) % Return only k-space |
---|
174 | % [DATA,msg]=aedes_readvnmr(filename,'return',3) % Return both image data and k-space |
---|
175 | % |
---|
176 | % % Read first data file header and then image data and k-space |
---|
177 | % [DATA,msg]=aedes_readvnmr(filename,'header') |
---|
178 | % [DATA,msg]=aedes_readvnmr(DATA.HDR,'return',3) |
---|
179 | % |
---|
180 | % % Read VNMR-data with default options and save slices in NIfTI |
---|
181 | % % format |
---|
182 | % DATA=aedes_readvnmr('example.fid','saved_slices.nii'); |
---|
183 | % |
---|
184 | % % If you want to use non-default options and also write |
---|
185 | % % NIfTI-files, use the property/value formalism, for example |
---|
186 | % DATA=aedes_readvnmr('example.fid','ZeroPadding','off',... |
---|
187 | % 'OutputFile','saved_slices.nii'); |
---|
188 | % |
---|
189 | % See also: |
---|
190 | % AEDES_READFIDPREFS, AEDES_READPROCPAR, AEDES_READPHASETABLE, |
---|
191 | % AEDES_DATA_READ, AEDES_WRITE_NIFTI, AEDES |
---|
192 | |
---|
193 | % This function is a part of Aedes - A graphical tool for analyzing |
---|
194 | % medical images |
---|
195 | % |
---|
196 | % Copyright (C) 2010 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uef.fi> |
---|
197 | % |
---|
198 | % Department of Physics, Department of Neurobiology |
---|
199 | % University of Kuopio, FINLAND |
---|
200 | % |
---|
201 | % This program may be used under the terms of the GNU General Public |
---|
202 | % License version 2.0 as published by the Free Software Foundation |
---|
203 | % and appearing in the file LICENSE.TXT included in the packaging of |
---|
204 | % this program. |
---|
205 | % |
---|
206 | % This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE |
---|
207 | % WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. |
---|
208 | |
---|
209 | |
---|
210 | if nargout==2 |
---|
211 | Dat.ShowErrors = false; |
---|
212 | msg_out=''; |
---|
213 | else |
---|
214 | Dat.ShowErrors = true; |
---|
215 | end |
---|
216 | |
---|
217 | %% --------------------- |
---|
218 | % Defaults |
---|
219 | Dat.ReturnKSpace = false; |
---|
220 | Dat.ReturnFTData = true; |
---|
221 | Dat.ReturnRawKspace = false; |
---|
222 | Dat.DCcorrection = false; |
---|
223 | Dat.ZeroPadding = 2; |
---|
224 | Dat.Sorting = true; |
---|
225 | Dat.UsePhaseTable = true; |
---|
226 | Dat.FastDataRead = true; |
---|
227 | Dat.precision = 'single'; |
---|
228 | Dat.FileChunkSize = 500; % Chunk size (in MB) for FastRead |
---|
229 | |
---|
230 | %% Other defaults |
---|
231 | Dat.ShowWaitbar = true; |
---|
232 | procpar=[]; |
---|
233 | Dat.phasetable=[]; |
---|
234 | Dat.OutputFile = false; |
---|
235 | Dat.ReorientEPI = false; |
---|
236 | Dat.RemoveEPIphaseIm = true; |
---|
237 | Dat.EPIBlockSize = 300; |
---|
238 | Dat.EPIPhasedArrayData = false; |
---|
239 | Dat.EPI_PC = 'auto'; |
---|
240 | Dat.OrientImages = true; |
---|
241 | Dat.SumOfSquares = 1; |
---|
242 | Dat.UseCustomRecon = true; |
---|
243 | Dat.SortPSS = true; |
---|
244 | |
---|
245 | % Default custom reconstruction file directory |
---|
246 | tmp = which('aedes'); |
---|
247 | if ~isempty(tmp) |
---|
248 | tmp_path=fileparts(tmp); |
---|
249 | recon_dir = [tmp_path,filesep,'vnmr_recon',filesep]; |
---|
250 | else |
---|
251 | % If all else fails, look in the current directory |
---|
252 | recon_dir = [pwd,filesep,'vnmr_recon',filesep]; |
---|
253 | end |
---|
254 | |
---|
255 | %% Set data format label |
---|
256 | DATA.DataFormat = 'vnmr'; |
---|
257 | |
---|
258 | %% Parse input arguments |
---|
259 | if nargin==0 || isempty(filename) |
---|
260 | |
---|
261 | %% Get default directory |
---|
262 | try |
---|
263 | defaultdir = getpref('Aedes','GetDataFileDir'); |
---|
264 | catch |
---|
265 | defaultdir = [pwd,filesep]; |
---|
266 | end |
---|
267 | |
---|
268 | %% If no input arguments are given, ask for a file |
---|
269 | [f_name, f_path, f_index] = uigetfile( ... |
---|
270 | {'fid;FID','Varian FID-files (fid)'; ... |
---|
271 | '*.*', 'All Files (*.*)'}, ... |
---|
272 | 'Select VNMR (Varian) FID file',defaultdir); |
---|
273 | if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed |
---|
274 | DATA=[]; |
---|
275 | msg_out='Canceled'; |
---|
276 | return |
---|
277 | end |
---|
278 | ReadHdr = true; |
---|
279 | ReadData = true; |
---|
280 | filename=[f_path,f_name]; |
---|
281 | |
---|
282 | %% Set default directory to preferences |
---|
283 | setpref('Aedes','GetDataFileDir',f_path) |
---|
284 | |
---|
285 | end |
---|
286 | |
---|
287 | if nargin==1 |
---|
288 | if isstruct(filename) |
---|
289 | hdr = filename; |
---|
290 | filename = [hdr.fpath,hdr.fname]; |
---|
291 | ReadHdr = false; |
---|
292 | ReadData = true; |
---|
293 | elseif ischar(filename) |
---|
294 | ReadHdr = true; |
---|
295 | ReadData = true; |
---|
296 | end |
---|
297 | elseif nargin==2 |
---|
298 | if strcmpi(varargin{1},'header') |
---|
299 | ReadHdr = true; |
---|
300 | ReadData = false; |
---|
301 | elseif ischar(varargin{1}) |
---|
302 | ReadHdr = true; |
---|
303 | ReadData = true; |
---|
304 | Dat.OutputFile = varargin{1}; |
---|
305 | else |
---|
306 | if ~Dat.ShowErrors |
---|
307 | DATA=[]; |
---|
308 | msg_out=sprintf('Invalid second input argument "%s".',varargin{1}); |
---|
309 | return |
---|
310 | else |
---|
311 | error('Invalid second input argument "%s".',varargin{1}) |
---|
312 | end |
---|
313 | end |
---|
314 | else |
---|
315 | if isstruct(filename) |
---|
316 | hdr = filename; |
---|
317 | filename = [hdr.fpath,hdr.fname]; |
---|
318 | ReadHdr = false; |
---|
319 | ReadData = true; |
---|
320 | elseif isempty(filename) |
---|
321 | [f_name, f_path, f_index] = uigetfile( ... |
---|
322 | {'fid;FID','Varian FID-files (fid)'; ... |
---|
323 | '*.*', 'All Files (*.*)'}, ... |
---|
324 | 'Select VNMR (Varian) FID file'); |
---|
325 | if ( all(f_name==0) || all(f_path==0) ) % Cancel is pressed |
---|
326 | DATA=[]; |
---|
327 | msg_out='Canceled'; |
---|
328 | return |
---|
329 | end |
---|
330 | ReadHdr = true; |
---|
331 | ReadData = true; |
---|
332 | filename=[f_path,f_name]; |
---|
333 | else |
---|
334 | ReadHdr = true; |
---|
335 | ReadData = true; |
---|
336 | end |
---|
337 | |
---|
338 | for ii=1:2:length(varargin) |
---|
339 | switch lower(varargin{ii}) |
---|
340 | case 'return' |
---|
341 | if length(varargin)<2 |
---|
342 | DATA=[]; |
---|
343 | if ~Dat.ShowErrors |
---|
344 | msg_out='"Return" value not specified!'; |
---|
345 | else |
---|
346 | error('"Return" value not specified!') |
---|
347 | end |
---|
348 | return |
---|
349 | else |
---|
350 | if varargin{ii+1}==1 |
---|
351 | Dat.ReturnKSpace = false; |
---|
352 | Dat.ReturnFTData = true; |
---|
353 | elseif varargin{ii+1}==2 |
---|
354 | Dat.ReturnKSpace = true; |
---|
355 | Dat.ReturnFTData = false; |
---|
356 | elseif varargin{ii+1}==3 |
---|
357 | Dat.ReturnKSpace = true; |
---|
358 | Dat.ReturnFTData = true; |
---|
359 | elseif varargin{ii+1}==4 |
---|
360 | Dat.ReturnRawKspace = true; |
---|
361 | Dat.ReturnKSpace = true; |
---|
362 | Dat.ReturnFTData = false; |
---|
363 | end |
---|
364 | end |
---|
365 | |
---|
366 | case {'dc','dccorrection','dccorr'} |
---|
367 | if strcmpi(varargin{ii+1},'on') |
---|
368 | Dat.DCcorrection = true; |
---|
369 | else |
---|
370 | Dat.DCcorrection = false; |
---|
371 | end |
---|
372 | |
---|
373 | case 'procpar' |
---|
374 | procpar=varargin{ii+1}; |
---|
375 | |
---|
376 | case 'zeropadding' |
---|
377 | if ischar(varargin{ii+1}) |
---|
378 | if strcmpi(varargin{ii+1},'on') |
---|
379 | Dat.ZeroPadding = 1; % on |
---|
380 | elseif strcmpi(varargin{ii+1},'off') |
---|
381 | Dat.ZeroPadding = 0; % off |
---|
382 | else |
---|
383 | Dat.ZeroPadding = 2; % auto |
---|
384 | end |
---|
385 | else |
---|
386 | % Undocumented |
---|
387 | Dat.ZeroPadding = 3; % Custom |
---|
388 | Dat.CustomZeroPadding = varargin{ii+1}; |
---|
389 | end |
---|
390 | |
---|
391 | case {'sumofsquares','soq'} |
---|
392 | Dat.SumOfSquares = varargin{ii+1}; |
---|
393 | |
---|
394 | case 'phasetable' |
---|
395 | Dat.phasetable = varargin{ii+1}; |
---|
396 | |
---|
397 | case 'sorting' |
---|
398 | if strcmpi(varargin{ii+1},'on') |
---|
399 | Dat.UsePhaseTable = true; |
---|
400 | Dat.Sorting = true; |
---|
401 | else |
---|
402 | Dat.UsePhaseTable = false; |
---|
403 | Dat.Sorting = false; |
---|
404 | end |
---|
405 | |
---|
406 | case 'fastread' |
---|
407 | if strcmpi(varargin{ii+1},'on') |
---|
408 | Dat.FastDataRead = true; |
---|
409 | else |
---|
410 | Dat.FastDataRead = false; |
---|
411 | end |
---|
412 | |
---|
413 | case 'filechunksize' |
---|
414 | Dat.FileChunkSize = round(varargin{ii+1}); |
---|
415 | |
---|
416 | case 'wbar' |
---|
417 | if strcmpi(varargin{ii+1},'on') |
---|
418 | Dat.ShowWaitbar = 1; |
---|
419 | else |
---|
420 | Dat.ShowWaitbar = 0; |
---|
421 | end |
---|
422 | |
---|
423 | case 'precision' |
---|
424 | if strcmpi(varargin{ii+1},'single') |
---|
425 | Dat.precision = 'single'; |
---|
426 | elseif strcmpi(varargin{ii+1},'double') |
---|
427 | Dat.precision = 'double'; |
---|
428 | else |
---|
429 | warning('Unsupported precision "%s". Using single...',varargin{ii+1}); |
---|
430 | Dat.precision = 'single'; |
---|
431 | end |
---|
432 | |
---|
433 | case 'outputfile' |
---|
434 | Dat.OutputFile = varargin{ii+1}; |
---|
435 | |
---|
436 | case 'reorientepi' |
---|
437 | if strcmpi(varargin{ii+1},'on') |
---|
438 | Dat.ReorientEPI = true; |
---|
439 | end |
---|
440 | |
---|
441 | case 'removeepiphaseim' |
---|
442 | if strcmpi(varargin{ii+1},'on') |
---|
443 | Dat.RemoveEPIphaseIm = true; |
---|
444 | end |
---|
445 | case 'epiblocksize' |
---|
446 | Dat.EPIBlockSize = round(varargin{ii+1}); |
---|
447 | |
---|
448 | case 'epiphasedarraydata' |
---|
449 | if strcmpi(varargin{ii+1},'on') |
---|
450 | Dat.EPIPhasedArrayData = true; |
---|
451 | else |
---|
452 | Dat.EPIPhasedArrayData = false; |
---|
453 | end |
---|
454 | case 'epi_pc' |
---|
455 | |
---|
456 | if any(strcmpi(varargin{ii+1},... |
---|
457 | {'off','auto','triple','pointwise'})) |
---|
458 | Dat.EPI_PC = varargin{ii+1}; |
---|
459 | else |
---|
460 | if ~Dat.ShowErrors |
---|
461 | msg_out=sprintf('Unknown EPI phase correction type "%s".',varargin{ii+1}); |
---|
462 | else |
---|
463 | error('Unknown EPI phase correction type "%s".',varargin{ii+1}) |
---|
464 | end |
---|
465 | return |
---|
466 | end |
---|
467 | case 'orientimages' |
---|
468 | if strcmpi(varargin{ii+1},'off') |
---|
469 | Dat.OrientImages = false; |
---|
470 | end |
---|
471 | otherwise |
---|
472 | DATA=[]; |
---|
473 | if ~Dat.ShowErrors |
---|
474 | msg_out = ['Unknown property "' varargin{ii} '"']; |
---|
475 | else |
---|
476 | error(['Unknown property "' varargin{ii} '"']) |
---|
477 | end |
---|
478 | return |
---|
479 | end |
---|
480 | end |
---|
481 | end |
---|
482 | |
---|
483 | |
---|
484 | |
---|
485 | % Parse filename |
---|
486 | [fpath,fname,fext]=fileparts(filename); |
---|
487 | if ~strcmpi(fname,'fid') |
---|
488 | if isempty(fname) |
---|
489 | fpath = [fpath,filesep]; |
---|
490 | else |
---|
491 | if isempty(fpath) |
---|
492 | fpath = [pwd,filesep,fname,fext,filesep]; |
---|
493 | else |
---|
494 | fpath = [fpath,filesep,fname,fext,filesep]; |
---|
495 | end |
---|
496 | end |
---|
497 | fname = 'fid'; |
---|
498 | else |
---|
499 | fpath = [fpath,filesep]; |
---|
500 | end |
---|
501 | |
---|
502 | %% Read procpar if not given as input argument |
---|
503 | if isempty(procpar) || nargin~=2 |
---|
504 | [procpar,proc_msg]=aedes_readprocpar([fpath,'procpar']); |
---|
505 | if isempty(procpar) |
---|
506 | DATA=[]; |
---|
507 | if ~Dat.ShowErrors |
---|
508 | msg_out=proc_msg; |
---|
509 | else |
---|
510 | error(proc_msg) |
---|
511 | end |
---|
512 | return |
---|
513 | end |
---|
514 | end |
---|
515 | |
---|
516 | % Don't use DC correction if number of transients is greater that 1 |
---|
517 | if procpar.nt>1 |
---|
518 | Dat.DCcorrection = false; |
---|
519 | end |
---|
520 | |
---|
521 | %% Read phasetable if necessary |
---|
522 | if Dat.Sorting |
---|
523 | |
---|
524 | % Look in preferences for tablib-directory |
---|
525 | if isfield(procpar,'petable') && ~strcmpi(procpar.petable{1},'n') |
---|
526 | try |
---|
527 | tabpath = getpref('Aedes','TabLibDir'); |
---|
528 | catch |
---|
529 | % If the table path was not found in preferences, try to use aedes |
---|
530 | % directory |
---|
531 | tmp_path = which('aedes'); |
---|
532 | if ~isempty(tmp_path) |
---|
533 | aedes_path=fileparts(tmp_path); |
---|
534 | tabpath = [aedes_path,filesep,'tablib',filesep]; |
---|
535 | else |
---|
536 | % If all else fails, look in the current directory |
---|
537 | tabpath = [pwd,filesep,'tablib',filesep]; |
---|
538 | end |
---|
539 | end |
---|
540 | [phasetable,msg] = aedes_readphasetable([tabpath,procpar.petable{1}]); |
---|
541 | else |
---|
542 | phasetable = []; |
---|
543 | end |
---|
544 | |
---|
545 | % If petable cannot be found, check if it is in procpar... |
---|
546 | if isempty(phasetable) && isfield(procpar,'pe_table') |
---|
547 | phasetable = {'t1',procpar.pe_table(:)}; |
---|
548 | elseif isempty(phasetable) && isfield(procpar,'pelist') && ... |
---|
549 | ~isempty(procpar.pelist) && isnumeric(procpar.pelist) && ... |
---|
550 | length(procpar.pelist)>1 |
---|
551 | phasetable = {'t1',reshape(procpar.pelist,procpar.etl,[]).'}; |
---|
552 | end |
---|
553 | |
---|
554 | if ~isempty(phasetable) |
---|
555 | Dat.phasetable = phasetable{1,2}; |
---|
556 | else |
---|
557 | Dat.phasetable = []; |
---|
558 | end |
---|
559 | |
---|
560 | %% Abort and throw error if phasetable cannot be read and the FID-file |
---|
561 | % has not been sorted |
---|
562 | % if isempty(phasetable) && ~isfield(procpar,'flash_converted') |
---|
563 | % DATA=[]; |
---|
564 | % if ~Dat.ShowErrors |
---|
565 | % msg_out={['Could not find the required phase table "' ... |
---|
566 | % procpar.petable{1} '" in the following folders'],... |
---|
567 | % ['"' tabpath '".']}; |
---|
568 | % else |
---|
569 | % error('Could not find the required phase table "%s" in %s',... |
---|
570 | % procpar.petable{1},tabpath) |
---|
571 | % end |
---|
572 | % return |
---|
573 | % elseif ( isempty(phasetable) && isfield(procpar,'flash_converted') ) || ... |
---|
574 | % ~Dat.UsePhaseTable |
---|
575 | % %% If the FID-file has been sorted, the reading can continue but |
---|
576 | % % throw a warning. |
---|
577 | % fprintf(1,'Warning: aedes_readfid: Could not find phasetable "%s" in "%s"!\n',procpar.petable{1},tabpath) |
---|
578 | % Dat.phasetable=[]; |
---|
579 | % else |
---|
580 | % Dat.phasetable = phasetable{1,2}; |
---|
581 | % end |
---|
582 | end |
---|
583 | |
---|
584 | % Convert phasetable indices to MATLAB indices |
---|
585 | if ~isempty(Dat.phasetable) |
---|
586 | Dat.phasetable=Dat.phasetable+(-min(min(Dat.phasetable))+1); |
---|
587 | else |
---|
588 | Dat.UsePhaseTable = false; |
---|
589 | end |
---|
590 | |
---|
591 | %% Open FID-file |
---|
592 | [file_fid,msg]=fopen([fpath,fname],'r','ieee-be'); |
---|
593 | if file_fid < 0 |
---|
594 | DATA=[]; |
---|
595 | if ~Dat.ShowErrors |
---|
596 | msg_out=['Could not open file "' filename '" for reading.']; |
---|
597 | else |
---|
598 | error('Could not open file "%s" for reading.',filename) |
---|
599 | end |
---|
600 | return |
---|
601 | end |
---|
602 | |
---|
603 | % Read only header |
---|
604 | if ReadHdr && ~ReadData |
---|
605 | [hdr,msg_out]=l_ReadDataFileHeader(file_fid); |
---|
606 | if ~isempty(msg_out) |
---|
607 | DATA=[]; |
---|
608 | fclose(file_fid); |
---|
609 | if Dat.ShowErrors |
---|
610 | error(msg_out) |
---|
611 | end |
---|
612 | return |
---|
613 | else |
---|
614 | DATA.HDR.FileHeader = hdr.FileHeader; |
---|
615 | DATA.FTDATA=[]; |
---|
616 | DATA.KSPACE=[]; |
---|
617 | DATA.PROCPAR=[]; |
---|
618 | DATA.PHASETABLE=[]; |
---|
619 | end |
---|
620 | elseif ~ReadHdr && ReadData % Header structure given as input argument |
---|
621 | % Read only data. |
---|
622 | [hdr,data,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat,procpar); |
---|
623 | if ~isempty(msg_out) |
---|
624 | DATA=[]; |
---|
625 | fclose(file_fid); |
---|
626 | if Dat.ShowErrors |
---|
627 | error(msg_out) |
---|
628 | end |
---|
629 | return |
---|
630 | else |
---|
631 | DATA.HDR.FileHeader=hdr.FileHeader; |
---|
632 | DATA.HDR.BlockHeaders = hdr.BlockHeaders; |
---|
633 | DATA.FTDATA=data; |
---|
634 | DATA.KSPACE=kspace; |
---|
635 | DATA.PROCPAR=procpar; |
---|
636 | DATA.PHASETABLE=Dat.phasetable; |
---|
637 | end |
---|
638 | elseif ReadHdr && ReadData % Read headers and data |
---|
639 | [hdr,msg_out]=l_ReadDataFileHeader(file_fid); |
---|
640 | [hdr,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat); |
---|
641 | if ~isempty(msg_out) |
---|
642 | DATA=[]; |
---|
643 | fclose(file_fid); |
---|
644 | if Dat.ShowErrors |
---|
645 | error(msg_out) |
---|
646 | end |
---|
647 | return |
---|
648 | else |
---|
649 | DATA.HDR.FileHeader=hdr.FileHeader; |
---|
650 | DATA.HDR.BlockHeaders = hdr.BlockHeaders; |
---|
651 | DATA.FTDATA=[]; |
---|
652 | DATA.KSPACE=[]; |
---|
653 | DATA.PROCPAR=procpar; |
---|
654 | DATA.PHASETABLE=Dat.phasetable; |
---|
655 | end |
---|
656 | end |
---|
657 | |
---|
658 | % Close file |
---|
659 | fclose(file_fid); |
---|
660 | |
---|
661 | |
---|
662 | % Set file name and path to the HDR structure |
---|
663 | DATA.HDR.fname=fname; |
---|
664 | DATA.HDR.fpath=fpath; |
---|
665 | |
---|
666 | % Set parameter values |
---|
667 | if ~Dat.DCcorrection |
---|
668 | DATA.HDR.Param.DCcorrection = 'off'; |
---|
669 | else |
---|
670 | DATA.HDR.Param.DCcorrection = 'on'; |
---|
671 | end |
---|
672 | |
---|
673 | % Reconstruct kspace ======================================= |
---|
674 | if Dat.ReturnRawKspace |
---|
675 | % If asked to return raw unsorted kspace, return immediately |
---|
676 | DATA.KSPACE=kspace; |
---|
677 | return |
---|
678 | else |
---|
679 | |
---|
680 | % Check if the sequence used has a custom reconstruct code |
---|
681 | recon_func=l_GetReconFunc(recon_dir); |
---|
682 | recon_func_ind = 0; |
---|
683 | if ~isempty(recon_func) |
---|
684 | % Get sequence name |
---|
685 | seqfil = procpar.seqfil; |
---|
686 | |
---|
687 | % Check if any of the custom reconstruction files would |
---|
688 | % like to do the reconstruction... |
---|
689 | for ii=1:length(recon_func) |
---|
690 | try |
---|
691 | list=recon_func{ii}(); |
---|
692 | catch |
---|
693 | warning('Error in custom reconsruction function "%s", skipping...',func2str(recon_func{ii})); |
---|
694 | continue |
---|
695 | end |
---|
696 | if any(strcmp(seqfil,list)) |
---|
697 | recon_func_ind = ii; |
---|
698 | break |
---|
699 | end |
---|
700 | end |
---|
701 | end |
---|
702 | |
---|
703 | if recon_func_ind==0 |
---|
704 | Dat.UseCustomRecon = false; |
---|
705 | end |
---|
706 | |
---|
707 | if Dat.UseCustomRecon |
---|
708 | if Dat.ShowWaitbar |
---|
709 | wbh=aedes_calc_wait(sprintf('%s\n%s',... |
---|
710 | ['Using custom function ',upper(func2str(recon_func{recon_func_ind}))],... |
---|
711 | ['to reconstruct sequence ',procpar.seqfil{1}])); |
---|
712 | drawnow |
---|
713 | end |
---|
714 | [kspace,data,msg_out]=recon_func{recon_func_ind}(kspace,Dat,procpar); |
---|
715 | if Dat.ShowWaitbar |
---|
716 | close(wbh) |
---|
717 | end |
---|
718 | if isempty(data) && Dat.ReturnFTData |
---|
719 | % Fourier transform data if not done in custom reconstruction code |
---|
720 | [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar); |
---|
721 | end |
---|
722 | else |
---|
723 | % Use the default reconstruction code |
---|
724 | [kspace,msg_out,procpar]=l_ReconstructKspace(kspace,Dat,procpar); |
---|
725 | DATA.PROCPAR = procpar; |
---|
726 | |
---|
727 | % Fourier transform data |
---|
728 | if Dat.ReturnFTData |
---|
729 | [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar); |
---|
730 | end |
---|
731 | end |
---|
732 | end |
---|
733 | if Dat.ReturnKSpace |
---|
734 | DATA.KSPACE = kspace; |
---|
735 | else |
---|
736 | DATA.KSPACE = []; |
---|
737 | end |
---|
738 | if Dat.ReturnFTData |
---|
739 | DATA.FTDATA = data; |
---|
740 | else |
---|
741 | DATA.FTDATA = []; |
---|
742 | end |
---|
743 | clear kspace data |
---|
744 | |
---|
745 | |
---|
746 | return |
---|
747 | |
---|
748 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
749 | % Read Data File (Main) Header |
---|
750 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
751 | function [hdr,msg_out]=l_ReadDataFileHeader(file_fid) |
---|
752 | |
---|
753 | msg_out=''; |
---|
754 | %% Read Data File Header |
---|
755 | try |
---|
756 | FH.nblocks = fread(file_fid,1,'long'); % Number of blocks in file |
---|
757 | FH.ntraces = fread(file_fid,1,'long'); % Number of traces per block |
---|
758 | FH.np = fread(file_fid,1,'long'); % Number of elements per trace |
---|
759 | FH.ebytes = fread(file_fid,1,'long'); % Number of bytes per element |
---|
760 | FH.tbytes = fread(file_fid,1,'long'); % Number of bytes per trace |
---|
761 | FH.bbytes = fread(file_fid,1,'long'); % Number of bytes per block |
---|
762 | FH.vers_id = fread(file_fid,1,'short'); % Software version, file_id status bits |
---|
763 | FH.status = fread(file_fid,1,'short'); % Status of whole file |
---|
764 | FH.nbheaders = fread(file_fid,1,'long'); % Number of block headers per block |
---|
765 | |
---|
766 | hdr.FileHeader = FH; |
---|
767 | |
---|
768 | %% Parse status bits |
---|
769 | hdr.FileHeader.status=[]; |
---|
770 | tmp_str = {'S_DATA',... % 0 = no data, 1 = data |
---|
771 | 'S_SPEC',... % 0 = FID, 1 = spectrum |
---|
772 | 'S_32',... % 0 = 16 bit, 1 = 32 bit integer |
---|
773 | 'S_FLOAT',... % 0 = integer, 1 = floating point |
---|
774 | 'S_COMPLEX',... % 0 = real, 1 = complex |
---|
775 | 'S_HYPERCOMPLEX',... % 1 = hypercomplex |
---|
776 | '',... % Unused bit |
---|
777 | 'S_ACQPAR',... % 0 = not Acqpar, 1 = Acqpar |
---|
778 | 'S_SECND',... % 0 = first FT, 1 = second FT |
---|
779 | 'S_TRANSF',... % 0 = regular, 1 = transposed |
---|
780 | '',... % Unused bit |
---|
781 | 'S_NP',... % 1 = np dimension is active |
---|
782 | 'S_NF',... % 1 = nf dimension is active |
---|
783 | 'S_NI',... % 1 = ni dimension is active |
---|
784 | 'S_NI2',... % 1 = ni2 dimension is active |
---|
785 | ''... % Unused bit |
---|
786 | }; |
---|
787 | status_bits = fliplr(double(dec2bin(FH.status,16))==49); |
---|
788 | for ii=1:length(tmp_str) |
---|
789 | if ~isempty(tmp_str{ii}) |
---|
790 | hdr.FileHeader.status.(tmp_str{ii}) = status_bits(ii); |
---|
791 | end |
---|
792 | end |
---|
793 | catch |
---|
794 | hdr=[]; |
---|
795 | msg_out=['Error occurred while reading file header from "' filename '".']; |
---|
796 | return |
---|
797 | end |
---|
798 | |
---|
799 | |
---|
800 | |
---|
801 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
802 | % Read Block Headers and Data |
---|
803 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
804 | function [hdr,kspace,msg_out]=l_ReadBlockData(file_fid,hdr,Dat) |
---|
805 | |
---|
806 | hdr.BlockHeaders = []; |
---|
807 | %hdr.HypercmplxBHeaders = []; |
---|
808 | kspace=[]; |
---|
809 | count = []; |
---|
810 | msg_out=''; |
---|
811 | |
---|
812 | |
---|
813 | BlockHeadStatusLabels = {... |
---|
814 | 'S_DATA',... % 0 = no data, 1 = data |
---|
815 | 'S_SPEC',... % 0 = FID, 1 = spectrum |
---|
816 | 'S_32',... % 0 = 16 bit, 1 = 32 bit integer |
---|
817 | 'S_FLOAT',... % 0 = integer, 1 = floating point |
---|
818 | 'S_COMPLEX',... % 0 = real, 1 = complex |
---|
819 | 'S_HYPERCOMPLEX',... % 1 = hypercomplex |
---|
820 | '',... % Unused bit |
---|
821 | 'MORE_BLOCKS',... % 0 = absent, 1 = present |
---|
822 | 'NP_CMPLX',... % 0 = real, 1 = complex |
---|
823 | 'NF_CMPLX',... % 0 = real, 1 = complex |
---|
824 | 'NI_CMPLX',... % 0 = real, 1 = complex |
---|
825 | 'NI2_CMPLX',... % 0 = real, 1 = complex |
---|
826 | '',... % Unused bit |
---|
827 | '',... % Unused bit |
---|
828 | '',... % Unuesd bit |
---|
829 | ''... % Unused bit |
---|
830 | }; |
---|
831 | |
---|
832 | BlockHeadModeLabels = {... |
---|
833 | 'NP_PHMODE',... % 1 = ph mode |
---|
834 | 'NP_AVMODE',... % 1 = av mode |
---|
835 | 'NP_PWRMODE',... % 1 = pwr mode |
---|
836 | '',... % Unused bit |
---|
837 | 'NF_PHMODE',... % 1 = ph mode |
---|
838 | 'NF_AVMODE',... % 1 = av mode |
---|
839 | 'NF_PWRMODE',... % 1 = pwr mode |
---|
840 | '',... % Unused bit |
---|
841 | 'NI_PHMODE',... % 1 = ph mode |
---|
842 | 'NI_AVMODE',... % 1 = av mode |
---|
843 | 'NI_PWRMODE',... % 1 = pwr mode |
---|
844 | '',... % Unused bit |
---|
845 | 'NI2_PHMODE',... % 1 = ph mode |
---|
846 | 'NI2_AVMODE',... % 1 = av mode |
---|
847 | 'NI2_PWRMODE',... % 1 = pwr mode |
---|
848 | ''... % Unused bit |
---|
849 | }; |
---|
850 | |
---|
851 | |
---|
852 | % The nbheaders -field is not correct in some cases. Thus, this field |
---|
853 | % cannot be trusted and the real nbheaders has to be calculated from the |
---|
854 | % byte counts. |
---|
855 | tmp_bytes=hdr.FileHeader.bbytes-hdr.FileHeader.tbytes*hdr.FileHeader.ntraces; |
---|
856 | header_bytes=28; |
---|
857 | if rem(tmp_bytes,header_bytes)~=0 |
---|
858 | msg_out = 'Block header byte count does not match with file header'; |
---|
859 | return |
---|
860 | else |
---|
861 | nbheaders = tmp_bytes/28; |
---|
862 | end |
---|
863 | |
---|
864 | % Allocate space for kspace |
---|
865 | % kspace = complex(zeros(hdr.FileHeader.np/2,... |
---|
866 | % hdr.FileHeader.ntraces*hdr.FileHeader.nblocks,... |
---|
867 | % Dat.precision)); |
---|
868 | |
---|
869 | kspace = complex(zeros(hdr.FileHeader.np/2,... |
---|
870 | hdr.FileHeader.ntraces*hdr.FileHeader.nblocks,Dat.precision)); |
---|
871 | |
---|
872 | |
---|
873 | %% - The older robust (but also slower) way for reading the data. |
---|
874 | %% When the blocksize is relatively small, this is also quite fast. |
---|
875 | if ~Dat.FastDataRead |
---|
876 | |
---|
877 | % Initialize waitbar |
---|
878 | if Dat.ShowWaitbar |
---|
879 | wb_h = aedes_wbar(0/hdr.FileHeader.nblocks,... |
---|
880 | {['Reading VNMR data.'],... |
---|
881 | ['(Processing data block ' ... |
---|
882 | num2str(0) '/' num2str(hdr.FileHeader.nblocks) ')']}); |
---|
883 | end |
---|
884 | |
---|
885 | %% Read data and block headers |
---|
886 | for ii=1:hdr.FileHeader.nblocks |
---|
887 | %% Update waitbar |
---|
888 | if Dat.ShowWaitbar |
---|
889 | aedes_wbar(ii/hdr.FileHeader.nblocks,... |
---|
890 | wb_h,... |
---|
891 | {['Reading VNMR data.'],... |
---|
892 | ['(Processing data block ' ... |
---|
893 | num2str(ii) '/' num2str(hdr.FileHeader.nblocks) ')']}) |
---|
894 | end |
---|
895 | |
---|
896 | %% Read block header and hypercomplex header |
---|
897 | for kk=1:nbheaders |
---|
898 | %% Read block header |
---|
899 | if kk==1 |
---|
900 | hdr.BlockHeaders.scale = fread(file_fid,1,'short'); % Scaling factor |
---|
901 | tmp_status = fread(file_fid,1,'short'); % Status of data in block |
---|
902 | hdr.BlockHeaders.status = []; |
---|
903 | hdr.BlockHeaders.index = fread(file_fid,1,'short'); % Block index |
---|
904 | tmp_mode = fread(file_fid,1,'short'); % Mode of data in block |
---|
905 | hdr.BlockHeaders.mode = []; |
---|
906 | hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID |
---|
907 | hdr.BlockHeaders.lpval = fread(file_fid,1,'float'); % f2 (2D-f1) left phase in phasefile |
---|
908 | hdr.BlockHeaders.rpval = fread(file_fid,1,'float'); % f2 (2D-f1) right phase in phasefile |
---|
909 | hdr.BlockHeaders.lvl = fread(file_fid,1,'float'); % level drift correction |
---|
910 | hdr.BlockHeaders.tlt = fread(file_fid,1,'float'); % tilt drift correction |
---|
911 | |
---|
912 | %% Parse status and mode bits |
---|
913 | status_bits = fliplr(double(dec2bin(tmp_status,16))==49); |
---|
914 | mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49); |
---|
915 | for tt=1:length(BlockHeadStatusLabels) |
---|
916 | if ~isempty(BlockHeadStatusLabels{tt}) |
---|
917 | hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt); |
---|
918 | end |
---|
919 | if ~isempty(BlockHeadModeLabels{tt}) |
---|
920 | hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt); |
---|
921 | end |
---|
922 | end |
---|
923 | |
---|
924 | |
---|
925 | else %% Read hypercomplex header |
---|
926 | fread(file_fid,1,'short'); % Spare |
---|
927 | hdr.BlockHeaders.HCHstatus = fread(file_fid,1,'short'); |
---|
928 | fread(file_fid,1,'short'); % Spare |
---|
929 | fread(file_fid,1,'short'); % Spare |
---|
930 | fread(file_fid,1,'long'); % Spare |
---|
931 | hdr.BlockHeaders.HCHlpval1 = fread(file_fid,1,'float'); |
---|
932 | hdr.BlockHeaders.HCHrpval1 = fread(file_fid,1,'float'); |
---|
933 | fread(file_fid,1,'float'); % Spare |
---|
934 | fread(file_fid,1,'float'); % Spare |
---|
935 | end |
---|
936 | end |
---|
937 | |
---|
938 | %% Check block index to be sure about the data type |
---|
939 | if hdr.BlockHeaders.index~=ii |
---|
940 | fprintf(1,'Warning: Index mismatch in "%s" block %d\n',fopen(file_fid),ii); |
---|
941 | |
---|
942 | % Use information from the file header instead of the BlockHeader if |
---|
943 | % there is a mismatch in blockheader index... |
---|
944 | useFileHeader = true; |
---|
945 | else |
---|
946 | useFileHeader = false; |
---|
947 | end |
---|
948 | |
---|
949 | %% Determine data precision |
---|
950 | if useFileHeader |
---|
951 | if hdr.FileHeader.status.S_FLOAT==1 |
---|
952 | prec_str = ['single=>',Dat.precision]; |
---|
953 | elseif hdr.FileHeader.status.S_32==1 ... |
---|
954 | && hdr.FileHeader.status.S_FLOAT==0 |
---|
955 | prec_str = ['int32=>',Dat.precision]; |
---|
956 | elseif hdr.FileHeader.status.S_32==0 ... |
---|
957 | && hdr.FileHeader.status.S_FLOAT==0 |
---|
958 | prec_str = ['int16=>',Dat.precision]; |
---|
959 | end |
---|
960 | |
---|
961 | else |
---|
962 | if hdr.BlockHeaders.status.S_FLOAT==1 |
---|
963 | prec_str = ['single=>',Dat.precision]; |
---|
964 | elseif hdr.BlockHeaders.status.S_32==1 ... |
---|
965 | && hdr.BlockHeaders.status.S_FLOAT==0 |
---|
966 | prec_str = ['int32=>',Dat.precision]; |
---|
967 | elseif hdr.BlockHeaders.status.S_32==0 ... |
---|
968 | && hdr.BlockHeaders.status.S_FLOAT==0 |
---|
969 | prec_str = ['int16=>',Dat.precision]; |
---|
970 | end |
---|
971 | end |
---|
972 | |
---|
973 | % Read k-space |
---|
974 | tmp=fread(file_fid,... |
---|
975 | [hdr.FileHeader.np,hdr.FileHeader.ntraces],... |
---|
976 | prec_str); |
---|
977 | |
---|
978 | |
---|
979 | % Store complex block and perform DC correction |
---|
980 | if ~Dat.DCcorrection |
---|
981 | complex_block = complex(tmp(1:2:end,:),tmp(2:2:end,:)); |
---|
982 | else |
---|
983 | complex_block = complex(tmp(1:2:end,:)-hdr.BlockHeaders.lvl,... |
---|
984 | tmp(2:2:end,:)-hdr.BlockHeaders.tlt); |
---|
985 | end |
---|
986 | |
---|
987 | inds = ((ii-1)*hdr.FileHeader.ntraces+1):(ii*hdr.FileHeader.ntraces); |
---|
988 | kspace(:,inds) = complex_block; |
---|
989 | |
---|
990 | % Do not save blockheaders by default. When reading data files with a lot of |
---|
991 | % blocks (e.g. over 1000) the processing of the DATA structure can be |
---|
992 | % slowed down considerably. If you for some reason want to save also the |
---|
993 | % block headers in the DATA structure comment out the code line below. |
---|
994 | hdr.BlockHeaders = []; |
---|
995 | end % for ii=1:hdr. |
---|
996 | |
---|
997 | if Dat.ShowWaitbar |
---|
998 | delete(wb_h) |
---|
999 | end |
---|
1000 | |
---|
1001 | else |
---|
1002 | %% ------------------------------------------------------------------- |
---|
1003 | %% Fast Method for reading data. This may fail with some datas and can |
---|
1004 | %% also require a relatively large amount of memory. This |
---|
1005 | %% method should be used for EPI datas that contain a large number |
---|
1006 | %% of block headers... |
---|
1007 | |
---|
1008 | % Check the size of the FID-file |
---|
1009 | d=dir(fopen(file_fid)); |
---|
1010 | file_sz = d.bytes/1024/1024; % File size in MB |
---|
1011 | if file_sz<Dat.FileChunkSize |
---|
1012 | nBlocks = 1; |
---|
1013 | else |
---|
1014 | nBlocks = ceil(file_sz/Dat.FileChunkSize); % Read data in 500 MB blocks |
---|
1015 | end |
---|
1016 | |
---|
1017 | % Initialize waitbar |
---|
1018 | if Dat.ShowWaitbar |
---|
1019 | if nBlocks==1 |
---|
1020 | wb_h = aedes_calc_wait(['Reading VNMR data.']); |
---|
1021 | else |
---|
1022 | wb_h = aedes_wbar(1/nBlocks,{['Reading VNMR data.'],... |
---|
1023 | sprintf('Reading block 0/%d',nBlocks)}); |
---|
1024 | end |
---|
1025 | end |
---|
1026 | |
---|
1027 | % The first block header is read and it is assumed that the values in |
---|
1028 | % the other block headers don't change. |
---|
1029 | hdr.BlockHeaders.scale = fread(file_fid,1,'short'); % Scaling factor |
---|
1030 | tmp_status = fread(file_fid,1,'short'); % Status of data in block |
---|
1031 | hdr.BlockHeaders.status = []; |
---|
1032 | hdr.BlockHeaders.index = fread(file_fid,1,'short'); % Block index |
---|
1033 | tmp_mode = fread(file_fid,1,'short'); % Mode of data in block |
---|
1034 | hdr.BlockHeaders.mode = []; |
---|
1035 | hdr.BlockHeaders.ctcount = fread(file_fid,1,'long'); % ct value for FID |
---|
1036 | hdr.BlockHeaders.lpval = fread(file_fid,1,'float'); % f2 (2D-f1) left phase in phasefile |
---|
1037 | hdr.BlockHeaders.rpval = fread(file_fid,1,'float'); % f2 (2D-f1) right phase in phasefile |
---|
1038 | hdr.BlockHeaders.lvl = fread(file_fid,1,'float'); % level drift correction |
---|
1039 | hdr.BlockHeaders.tlt = fread(file_fid,1,'float'); % tilt drift correction |
---|
1040 | |
---|
1041 | %% Parse status and mode bits |
---|
1042 | status_bits = fliplr(double(dec2bin(tmp_status,16))==49); |
---|
1043 | mode_bits = fliplr(double(dec2bin(tmp_mode,16))==49); |
---|
1044 | for tt=1:length(BlockHeadStatusLabels) |
---|
1045 | if ~isempty(BlockHeadStatusLabels{tt}) |
---|
1046 | hdr.BlockHeaders.status.(BlockHeadStatusLabels{tt}) = status_bits(tt); |
---|
1047 | end |
---|
1048 | if ~isempty(BlockHeadModeLabels{tt}) |
---|
1049 | hdr.BlockHeaders.mode.(BlockHeadModeLabels{tt}) = mode_bits(tt); |
---|
1050 | end |
---|
1051 | end |
---|
1052 | |
---|
1053 | %% Determine data precision |
---|
1054 | if hdr.BlockHeaders.status.S_FLOAT==1 |
---|
1055 | prec_str = ['single=>',Dat.precision]; |
---|
1056 | prec = 4; % Precision in bytes |
---|
1057 | elseif hdr.BlockHeaders.status.S_32==1 ... |
---|
1058 | && hdr.BlockHeaders.status.S_FLOAT==0 |
---|
1059 | prec_str = ['int32=>',Dat.precision]; |
---|
1060 | prec = 4; |
---|
1061 | elseif hdr.BlockHeaders.status.S_32==0 ... |
---|
1062 | && hdr.BlockHeaders.status.S_FLOAT==0 |
---|
1063 | prec_str = ['int16=>',Dat.precision]; |
---|
1064 | prec = 2; |
---|
1065 | end |
---|
1066 | |
---|
1067 | % Seek the file back to the beginning of the first block header |
---|
1068 | fseek(file_fid,-28,0); |
---|
1069 | |
---|
1070 | % Determine the number of values that will result from block header(s) |
---|
1071 | nVals = (nbheaders*28)/prec; |
---|
1072 | |
---|
1073 | nbh = floor(hdr.FileHeader.nblocks/nBlocks); |
---|
1074 | szh = nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces; |
---|
1075 | for ii=1:nBlocks |
---|
1076 | if nBlocks~=1 |
---|
1077 | aedes_wbar(ii/nBlocks,wb_h,{['Reading VNMR data.'],... |
---|
1078 | sprintf('Reading block %d/%d',ii,nBlocks)}); |
---|
1079 | end |
---|
1080 | |
---|
1081 | % Read the whole data including block headers etc... |
---|
1082 | if ii==nBlocks |
---|
1083 | tmp = fread(file_fid,inf,prec_str); |
---|
1084 | else |
---|
1085 | tmp = fread(file_fid,nbh*szh,prec_str); |
---|
1086 | end |
---|
1087 | tmp=reshape(tmp,nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces,[]); |
---|
1088 | tmp(1:nVals,:)=[]; |
---|
1089 | tmp=reshape(tmp,hdr.FileHeader.np,[]); |
---|
1090 | |
---|
1091 | if ii==nBlocks |
---|
1092 | inds = ((ii-1)*nbh*hdr.FileHeader.ntraces+1):size(kspace,2); |
---|
1093 | else |
---|
1094 | inds = ((ii-1)*nbh*hdr.FileHeader.ntraces+1):ii*nbh*hdr.FileHeader.ntraces; |
---|
1095 | end |
---|
1096 | |
---|
1097 | % Do DC-correction if necessary |
---|
1098 | if ~Dat.DCcorrection |
---|
1099 | kspace(:,inds)=complex(tmp(1:2:end,:,:),tmp(2:2:end,:,:)); |
---|
1100 | else |
---|
1101 | kspace(:,inds)=complex(tmp(1:2:end,:,:)-hdr.BlockHeaders.lvl,... |
---|
1102 | tmp(2:2:end,:,:)-hdr.BlockHeaders.tlt); |
---|
1103 | end |
---|
1104 | end |
---|
1105 | clear tmp |
---|
1106 | |
---|
1107 | if Dat.ShowWaitbar |
---|
1108 | delete(wb_h) |
---|
1109 | pause(0.1) |
---|
1110 | end |
---|
1111 | |
---|
1112 | end |
---|
1113 | |
---|
1114 | % ================================================================== |
---|
1115 | % Reconstruct k-space |
---|
1116 | % ================================================================== |
---|
1117 | function [kspace,msg_out,procpar]=l_ReconstructKspace(kspace,Dat,procpar) |
---|
1118 | |
---|
1119 | msg_out = ''; |
---|
1120 | |
---|
1121 | % Get the strategic parameters |
---|
1122 | seqcon = procpar.seqcon{1}; |
---|
1123 | seqfil = procpar.seqfil{1}; |
---|
1124 | np = procpar.np; |
---|
1125 | nv = procpar.nv; |
---|
1126 | nv2 = procpar.nv2; |
---|
1127 | nv3 = procpar.nv3; |
---|
1128 | ns = procpar.ns; |
---|
1129 | if isfield(procpar,'etl') |
---|
1130 | etl = procpar.etl; |
---|
1131 | else |
---|
1132 | etl = 1; |
---|
1133 | end |
---|
1134 | pss = procpar.pss; |
---|
1135 | nt = procpar.nt; |
---|
1136 | nf = procpar.nf; |
---|
1137 | ne = procpar.ne; |
---|
1138 | if isfield(procpar,'flash_converted') |
---|
1139 | % Don't try to sort already sorted data... |
---|
1140 | Dat.Sorting = false; |
---|
1141 | seqcon(2:3) = 'sc'; |
---|
1142 | end |
---|
1143 | if isfield(procpar,'nseg') |
---|
1144 | nseg = procpar.nseg; |
---|
1145 | else |
---|
1146 | nseg = 1; |
---|
1147 | end |
---|
1148 | |
---|
1149 | |
---|
1150 | % Check dimensions |
---|
1151 | if nv==0 |
---|
1152 | % 1D-image (specrum) |
---|
1153 | AcqType = 1; |
---|
1154 | elseif nv2==0 |
---|
1155 | % 2D-image |
---|
1156 | AcqType = 2; |
---|
1157 | elseif nv3==0 |
---|
1158 | % 3D-image |
---|
1159 | AcqType = 3; |
---|
1160 | else |
---|
1161 | AcqType = 4; |
---|
1162 | end |
---|
1163 | |
---|
1164 | % Check number of receivers |
---|
1165 | if isfield(procpar,'rcvrs') && length(procpar.rcvrs{1})>1 |
---|
1166 | nRcvrs = length(find(procpar.rcvrs{1}=='y')); |
---|
1167 | else |
---|
1168 | nRcvrs = 1; |
---|
1169 | end |
---|
1170 | |
---|
1171 | |
---|
1172 | % Check for arrayed acquisition |
---|
1173 | if isfield(procpar,'array') && ~isempty(procpar.array{1}) |
---|
1174 | %if length(procpar.array)==1 && ~iscell(procpar.array{1}) && ... |
---|
1175 | % strcmp(procpar.array{1},'pad') && all(procpar.pad==0) |
---|
1176 | % % Skip the array parsing if the array is a dummy "pad"-array... |
---|
1177 | % isAcqArrayed = false; |
---|
1178 | % ArrayLength = 1; |
---|
1179 | %else |
---|
1180 | isAcqArrayed = true; |
---|
1181 | ArrayLength = []; |
---|
1182 | |
---|
1183 | % Determine array length |
---|
1184 | for ii=1:length(procpar.array) |
---|
1185 | if iscell(procpar.array{ii}) |
---|
1186 | ArrayLength(ii) = length(procpar.(procpar.array{ii}{1})); |
---|
1187 | else |
---|
1188 | ArrayLength(ii) = length(procpar.(procpar.array{ii})); |
---|
1189 | end |
---|
1190 | end |
---|
1191 | ArrayLength = prod(ArrayLength); |
---|
1192 | %end |
---|
1193 | else |
---|
1194 | isAcqArrayed = false; |
---|
1195 | ArrayLength = 1; |
---|
1196 | end |
---|
1197 | |
---|
1198 | |
---|
1199 | % Reconstruct k-space ---------------------- |
---|
1200 | if AcqType==1 |
---|
1201 | % Reconstruct 1D data ... |
---|
1202 | elseif AcqType==2 |
---|
1203 | |
---|
1204 | % Reconstruct 2D data |
---|
1205 | if strcmpi(seqcon,'ncsnn') |
---|
1206 | kspace = reshape(kspace,[np/2,etl,ns,nRcvrs,ArrayLength,nv/etl]); |
---|
1207 | kspace = permute(kspace,[1 2 6 3 5 4]); |
---|
1208 | kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]); |
---|
1209 | elseif strcmpi(seqcon,'nscnn') |
---|
1210 | if isfield(procpar,'flash_converted') |
---|
1211 | kspace = reshape(kspace,[np/2,nRcvrs,ArrayLength,ns,nv]); |
---|
1212 | kspace = permute(kspace,[1 5 4 3 2]); |
---|
1213 | kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]); |
---|
1214 | else |
---|
1215 | kspace = reshape(kspace,[np/2,nRcvrs,nv,ArrayLength,ns]); |
---|
1216 | kspace = permute(kspace,[1 3 5 4 2]); |
---|
1217 | kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]); |
---|
1218 | end |
---|
1219 | elseif strcmpi(seqcon,'nccnn') |
---|
1220 | kspace = reshape(kspace,[np/2,etl,ns,nv/etl,ArrayLength,nRcvrs]); |
---|
1221 | kspace = permute(kspace,[1 2 4 3 5 6]); |
---|
1222 | kspace = reshape(kspace,[np/2,nv,ns*ArrayLength,1,nRcvrs]); |
---|
1223 | end |
---|
1224 | elseif AcqType==3 |
---|
1225 | % Reconstruct 3D data |
---|
1226 | if strcmpi(seqcon,'nscsn') |
---|
1227 | kspace = reshape(kspace,[np/2,etl,nv/etl,nRcvrs,ArrayLength*nv2]); |
---|
1228 | kspace = permute(kspace,[1 2 3 5 4]); |
---|
1229 | kspace = reshape(kspace,[np/2,nv,nv2,ArrayLength,nRcvrs]); |
---|
1230 | end |
---|
1231 | if strcmpi(seqcon,'ccccn') |
---|
1232 | kspace = reshape(kspace,[np/2,etl,nRcvrs,ne,nv/etl,ArrayLength*nv2]); |
---|
1233 | kspace = permute(kspace,[1 2 5 6 4 3]); |
---|
1234 | kspace = reshape(kspace,[np/2,nv,nv2,ne*ArrayLength,nRcvrs]); |
---|
1235 | end |
---|
1236 | else |
---|
1237 | % Reconstruct nD data ... to be written |
---|
1238 | end |
---|
1239 | |
---|
1240 | |
---|
1241 | % Sort data using phasetable ------------------ |
---|
1242 | if Dat.Sorting && ~isempty(Dat.phasetable) |
---|
1243 | Dat.phasetable = Dat.phasetable.'; |
---|
1244 | kspace(:,Dat.phasetable(:),:,:,:)=kspace; |
---|
1245 | end |
---|
1246 | |
---|
1247 | % Sort interleaved 2D data using pss |
---|
1248 | [sorted_pss,I_pss]=sort(procpar.pss); |
---|
1249 | if ~ismember(procpar.pss,sorted_pss,'rows') && Dat.SortPSS |
---|
1250 | kspace = kspace(:,:,I_pss,:,:); |
---|
1251 | procpar.pss_orig = procpar.pss; |
---|
1252 | procpar.pss = sorted_pss; |
---|
1253 | end |
---|
1254 | |
---|
1255 | % ================================================================== |
---|
1256 | % Fourier transform k-space |
---|
1257 | % ================================================================== |
---|
1258 | function [data,msg_out]=l_CalculateFFT(kspace,Dat,procpar) |
---|
1259 | |
---|
1260 | msg_out = ''; |
---|
1261 | data = []; |
---|
1262 | |
---|
1263 | |
---|
1264 | % Check dimensions |
---|
1265 | if procpar.nv==0 |
---|
1266 | % 1D-image |
---|
1267 | AcqType = 1; |
---|
1268 | elseif procpar.nv2==0 |
---|
1269 | % 2D-image |
---|
1270 | AcqType = 2; |
---|
1271 | elseif procpar.nv3==0 |
---|
1272 | % 3D-image |
---|
1273 | AcqType = 3; |
---|
1274 | else |
---|
1275 | AcqType = 4; |
---|
1276 | end |
---|
1277 | |
---|
1278 | % Check number of receivers |
---|
1279 | if isfield(procpar,'rcvrs') && length(procpar.rcvrs{1})>1 |
---|
1280 | nRcvrs = length(find(procpar.rcvrs{1}=='y')); |
---|
1281 | else |
---|
1282 | nRcvrs = 1; |
---|
1283 | end |
---|
1284 | |
---|
1285 | % If zeropadding is requested, calculate the padded size |
---|
1286 | if Dat.ZeroPadding~=0 |
---|
1287 | if Dat.ZeroPadding==1 |
---|
1288 | % Zeropad to square |
---|
1289 | if AcqType==1 |
---|
1290 | padSize = procpar.np/2; |
---|
1291 | elseif AcqType==2 |
---|
1292 | padSize = ones(1,2)*procpar.np/2; |
---|
1293 | padSize(3) = size(kspace,3); |
---|
1294 | else |
---|
1295 | padSize = ones(1,3)*procpar.np/2; |
---|
1296 | end |
---|
1297 | else |
---|
1298 | % Zeropadding is on "auto", i.e. zeropad to FOV |
---|
1299 | lpe = procpar.lpe; |
---|
1300 | lpe2 = procpar.lpe2; |
---|
1301 | lro = procpar.lro; |
---|
1302 | if AcqType==2 |
---|
1303 | % 2D data |
---|
1304 | padSize = [procpar.np/2 ... |
---|
1305 | procpar.np/2*(lpe/lro) ... |
---|
1306 | size(kspace,3)]; |
---|
1307 | elseif AcqType==3 && lpe2~=0 |
---|
1308 | % 3D data |
---|
1309 | padSize = [procpar.np/2 ... |
---|
1310 | procpar.np/2*(lpe/lro) ... |
---|
1311 | procpar.np/2*(lpe2/lro)]; |
---|
1312 | end |
---|
1313 | end |
---|
1314 | ks_sz = [size(kspace,1) ... |
---|
1315 | size(kspace,2) ... |
---|
1316 | size(kspace,3)]; |
---|
1317 | padSize = round(padSize); |
---|
1318 | if any(padSize>ks_sz) |
---|
1319 | %kspace(padSize) = 0; |
---|
1320 | kspace(padSize(1),padSize(2),padSize(3)) = 0; |
---|
1321 | end |
---|
1322 | else |
---|
1323 | padSize = [size(kspace,1) ... |
---|
1324 | size(kspace,2) ... |
---|
1325 | size(kspace,3)]; |
---|
1326 | end |
---|
1327 | |
---|
1328 | % Allocate space for Fourier transformed data |
---|
1329 | if nRcvrs>1 && any(Dat.SumOfSquares==[1 2]) |
---|
1330 | data_sz = [padSize,size(kspace,4),size(kspace,5)+1]; |
---|
1331 | data = zeros(data_sz,Dat.precision); |
---|
1332 | else |
---|
1333 | data = zeros(size(kspace),Dat.precision); |
---|
1334 | end |
---|
1335 | %data = []; |
---|
1336 | %if strcmpi(Dat.precision,'single') |
---|
1337 | % data = single(data); |
---|
1338 | %end |
---|
1339 | |
---|
1340 | % Fourier transform data |
---|
1341 | if nRcvrs>1 && any(Dat.SumOfSquares==[1 2]) |
---|
1342 | ind = [2:size(data,5)]; |
---|
1343 | else |
---|
1344 | ind = [1:size(data,5)]; |
---|
1345 | end |
---|
1346 | if AcqType==1 |
---|
1347 | data(:,:,:,:,ind) = abs(fftshift(ifft(kspace,[],1),1)); |
---|
1348 | elseif AcqType==2 |
---|
1349 | data(:,:,:,:,ind) = abs(fftshift(fftshift(ifft(ifft(kspace,[],1),[],2),1),2)); |
---|
1350 | elseif AcqType==3 |
---|
1351 | data(:,:,:,:,ind) = abs(fftshift(fftshift(fftshift(ifft(ifft(ifft(kspace,[],1),[],2),[],3),1),2),3)); |
---|
1352 | end |
---|
1353 | |
---|
1354 | % Calculate sum-of-squares image |
---|
1355 | if nRcvrs>1 && any(Dat.SumOfSquares==[1 2]) |
---|
1356 | % Calculate sum-of-squares |
---|
1357 | data(:,:,:,:,1) = sqrt(sum(data(:,:,:,:,ind).^2,5)); |
---|
1358 | data=abs(data); |
---|
1359 | if Dat.SumOfSquares==1 |
---|
1360 | % Remove individual receiver data |
---|
1361 | data=data(:,:,:,:,1); |
---|
1362 | end |
---|
1363 | end |
---|
1364 | |
---|
1365 | % ================================================================== |
---|
1366 | % Find custom functions for VNMR k-space reconstruction |
---|
1367 | % ================================================================== |
---|
1368 | function recon_func=l_GetReconFunc(recon_dir) |
---|
1369 | |
---|
1370 | recon_func = {}; |
---|
1371 | |
---|
1372 | if ~isdir(recon_dir) |
---|
1373 | return |
---|
1374 | end |
---|
1375 | |
---|
1376 | dir_struct=dir(recon_dir); |
---|
1377 | recon_files = {dir_struct(~[dir_struct(:).isdir]).name}; |
---|
1378 | if isempty(recon_files) |
---|
1379 | return |
---|
1380 | end |
---|
1381 | |
---|
1382 | % Remove files that don't have .m extension |
---|
1383 | ind = regexpi(recon_files,'\.m$'); |
---|
1384 | recon_files = {recon_files{~cellfun('isempty',ind)}}; |
---|
1385 | |
---|
1386 | currentDir = pwd; |
---|
1387 | try |
---|
1388 | cd(recon_dir); |
---|
1389 | for ii=1:length(recon_files) |
---|
1390 | recon_func{ii}=str2func(recon_files{ii}(1:end-2)); |
---|
1391 | end |
---|
1392 | cd(currentDir); |
---|
1393 | catch |
---|
1394 | cd(currentDir); |
---|
1395 | end |
---|
1396 | |
---|
1397 | |
---|
1398 | |
---|