-
Notifications
You must be signed in to change notification settings - Fork 0
/
nhdr_nrrd_read.m
905 lines (792 loc) · 43.6 KB
/
nhdr_nrrd_read.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
% Read NHDR/NRRD files
%
% headerInfo = nhdr_nrrd_read(nrrdFileName, bReadData)
%
% Inputs:
% * nrrdFileName (char): path to NHDR/NRRD file, either a detached header,
% a detached header pointing to detached data files or a NRRD standalone
% file with header and data included
% * bReadData (bool): set to true to read the data and store it in
% headerInfo.data, set to false to just import the header without the data
%
% Outputs:
% * headerInfo (struct): contains all the fields specified in the file
% header. If data was read it is contained in the 'data' field. That
% structure can be fed to the nhdr_nrrd_write module as is and produce
% valid NHDR/NRRD files.
%
% Format definition:
% http://teem.sourceforge.net/nrrd/format.html
%
% A few supported NRRD features:
% - detached headers with all variants of 'data file:'
% - raw, txt/text/ascii, gz/gzip encodings
% - definition of space and orientation
% - handling of diffusion-weighted MRI data with '<key>:=<value>' lines
% 'modality:=DWMRI', 'DWMRI_b-value:=' and 'DWMRI_gradient_0000:=',
% 'DWMRI_gradient_0001:=', 'DWMRI_gradient_0002:=', etc.
% (see https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:Nrrd_format)
%
%
% Other features:
% - exits cleanly upon error (no file should be left open)
% - informative error messages whenever possible
% - warnings printed to console
% - NHDR/NRRD writer fills out missing fields by inspecting the data
% whenever possible
%
%
% Unsupported features:
%
% In general, any header field that we were unable to parse is reported
% in a message printed to the console. Specific examples of
% unsupported features include:
%
% - reading data along more than 1 dimension or along a dimension other
% than the slowest (last) axis specified by the optional <subdim>, as in
% 'data file: <format> <min> <max> <step> [<subdim>]' or 'data file: LIST
% [<subdim>]'
% - storing all the comments found in headers
% - hex encoding
% - bz2/bzip2 encoding
% - byte skip; can only accept -1 with raw encoding, 0 for all other
% encodings
% - line skip: can only accept 0
% - checking that field specifications of the form '<field>: <desc>'
% appear no more than once in the NRRD header (unlike '<key>:=<value>'
% lines)
%
% Date: July 2018
% Author: Gaetan Rensonnet
%
%
%
% Notes :
%
% This is intended as a more comprehensive and accurate implementation of
% the NRRD format specification than most of the Matlab scripts that have
% been proposed so far. We try to fail graciously when an unsupported
% feature of the NRRD format is encountered. One of our main contributions
% is to propose a writer module which is compatible with the read module,
% in that the output of one can be given as an argument to the other to
% read or produce equivalent NHDR/NRRD files. This is still a version with
% much room for improvement.
%
% The following contributions inspired parts of our Matlab read/write
% modules:
%
% 1. The body of the writer module was pretty much written from scratch but
% the general structure of the reader's main body is based on the Matlab
% functions maReadNrrdHeader.m and maReadNrrdData.m by marc@bwh.harvard.edu
% and kquintus@bwh.harvard.edu (unpublished). Many additions were made and
% a few bugs fixed.
%
% 2. Jeff Mather's NRRD reader
% (http://nl.mathworks.com/matlabcentral/fileexchange/34653-nrrd-format-file-reader)
% and
% http://jeffmatherphotography.com/dispatches/2012/02/writing-a-file-reader-in-matlab/)
% provided the auxiliary functions:
% - adjustEndian
% - getDatatype, which we renamed nrrd_getMatlabDataType and now throws
% a gracious error if it encounters 'block'-type data,
% - readData was adapted to include a cross-platform fix to delete
% temporary files when using gzip encoding. David Feng's fix used a
% Windows-specific command to delete temporary files
% (https://nl.mathworks.com/matlabcentral/fileexchange/50830-nrrd-format-file-reader).
%
% 3. mdcacio's nrrdWriter
% (https://nl.mathworks.com/matlabcentral/fileexchange/48621-nrrdwriter-filename--matrix--pixelspacing--origin--encoding-)
% provided the auxiliary functions:
% - writeData(): we got rid of the 'unexpected end of input stream when
% attempting to gunzip the file' error when using gzip encoding, which we
% later found had been fixed by Quan Chen independenlty and in a similar
% manner.
% - setDatatype(), which we renamed get_nrrd_datatype() and is the
% reciprocal of getDatatype
%
function headerInfo = nhdr_nrrd_read(nrrdFileName, bReadData)
[mainFpath,mainFileName,mainFext] = fileparts(nrrdFileName);
headerInfo = struct();
% default header:
headerInfo.content = mainFileName; % default value, overwritten if content field is set
headerInfo.data = [];
fidr = fopen(nrrdFileName, 'r');
if (fidr == -1)
error('ABORT: %s does not exist.\n', nrrdFileName);
end
try
if ~(strcmpi(mainFext, '.nhdr') || strcmpi(mainFext, '.nrrd'))
warning('%s looks like a %s file, not a nhdr or nrrd file.\n', nrrdFileName, mainFext );
end
% Magic line
cs = fgetl(fidr);
assert(numel(cs) >= 8 && isequal(cs(1:4), 'NRRD'),...
'Bad signature. First line should be magic line of type "NRRD000X" with X an integer between 1 and 5.'); % FIXME should throw an error if a bad integer is provided
nrrd_version = sscanf(cs(5:end), '%d');
if nrrd_version > 5
error('This reader only supports versions of the NRRD file format up to 5. Detected %d.', nrrd_version)
end
% Always optional: defining orientation
define_orientation = 0; % internal-use variable
% Parse header
while ~feof(fidr)
cs = fgetl(fidr); % content string
if isempty(cs)
% End of header
break;
end
if foundKeyword( 'CONTENT:', cs )
headerInfo.content = strtrim( cs( length('CONTENT:')+1:end ) );
elseif foundKeyword('TYPE:', cs )
headerInfo.type = strtrim( cs( length('TYPE:')+1:end ) );
elseif foundKeyword('ENDIAN:', cs )
headerInfo.endian = strtrim( cs( length('ENDIAN:')+1:end ) );
elseif foundKeyword('ENCODING:', cs )
headerInfo.encoding = strtrim( cs( length('ENCODING:')+1:end ) );
elseif foundKeyword('DIMENSION:', cs )
headerInfo.dimension = sscanf( cs( length('DIMENSION:')+1:end ), '%i' );
elseif foundKeyword('SIZES:', cs )
iSizes = sscanf( cs(length('SIZES:')+1:end), '%i' );
headerInfo.sizes = iSizes(:)'; % store as row vector
elseif foundKeyword('KINDS:', cs )
headerInfo.kinds = extractStringList( cs(length('KINDS:')+1:end) ); % bug fixed with extractStringList where 2 entries are present
% FIXME: check that axis sizes match each kind according to nrrd standard
elseif foundKeyword('SPACE:', cs )
% Starts defining orientation (either space or space dimension,
% not both)
define_orientation = 1;
if isfield(headerInfo, 'spacedimension')
fprintf(['WARNING nhdr_nrrd_read %s:\n ''space'' field specifier will ' ...
'be checked for consistency but will be ignored afterwards ' ...
'because ''space dimension'' was specified before.\n'], fopen(fidr));
end
tmp_space = strtrim( cs( length('SPACE:')+1:end ) );
tmp_spacedimension = nrrd_getSpaceDimensions(tmp_space);
if tmp_spacedimension <= 0
error('%s: unrecognized ''space'' descriptor ''%s''.', fopen(fidr), tmp_space)
end
if isfield(headerInfo, 'spacedimension')
% internal_spacedimension already set
if internal_spacedimension ~= tmp_spacedimension
error(['%s: ''space'' field specifier implies a spatial ' ...
'(world) dimension equal to %d, which differs from ' ...
'the ''space dimension'' field specifier set to %d.'],...
fopen(fidr), tmp_spacedimension, internal_spacedimension)
end
% if no inconsistencies found, just ignore space field
else
% Set space info for the first time:
headerInfo.space = tmp_space;
internal_spacedimension = tmp_spacedimension; % for internal use only
end
elseif foundKeyword('SPACE DIMENSION:', cs)
% Starts defining orientation (either space or space dimension,
% not both)
define_orientation = 1;
if isfield(headerInfo, 'space')
fprintf(['WARNING nhdr_nrrd_read %s:\n ''space dimension'' field specifier ' ...
' will be checked for consistency but will be ignored afterwards ' ...
'because ''space'' was specified before.\n'],...
fopen(fidr));
end
tmp_spacedimension = sscanf( cs( length('SPACE DIMENSION:')+1:end), '%i' );
if numel(tmp_spacedimension) ~= 1
error(['%s: ''space dimension'' should be specified as one' ...
' integer number. Found %d element(s) instead.'],...
fopen(fidr), numel(tmp_spacedimension))
end
if tmp_spacedimension <= 0
error('%s: ''space dimension'' should be specified as a strictly positive integer (found %d).',...
fopen(fidr), tmp_spacedimension)
end
if isfield(headerInfo, 'space')
if tmp_spacedimension ~= internal_spacedimension
error(['%s: ''space dimension'' field specifier set to %d, ' ...
'which differs from the space (world) dimension implied by the ' ...
'''space'' field specifier which is equal to %d.'],...
fopen(fidr), tmp_spacedimension, internal_spacedimension)
end
% if no inconsistencies found, just ignore space dimension
% field
else
% Set space info for the first time:
headerInfo.spacedimension = tmp_spacedimension;
internal_spacedimension = tmp_spacedimension; % for internal use only
end
elseif foundKeyword('SPACE DIRECTIONS:', cs )
% Required if orientation defined but must come after space or
% space dimension
% space directions: <vector[0]> <vector[1]> ... <vector[dim-1]>
% The format of the <vector> is as follows. The vector is
% delimited by "(" and ")", and the individual components are
% comma-separated.
if ~define_orientation
error('%s: field specifier ''space directions'' cannot be set before ''space'' or ''space dimension''.',fopen(fidr))
end
space_dir_tmp = strtrim(cs(length('SPACE DIRECTIONS:')+1:end)); % remove leading and trailing spaces
spacedir_vecs = strsplit(space_dir_tmp); % cell array of strings after split at {' ','\f','\n','\r','\t','\v'}
SD_data = zeros(internal_spacedimension, internal_spacedimension);
% Check each vector: either none or (f1,...,f_spaceDim) with fi
% a floating-point number
cnt_space_vectors = 0;
for i = 1:length(spacedir_vecs)
none_start_index = strfind(lower(spacedir_vecs{i}), 'none');
if ~isempty(none_start_index)
% Axis-specific entry contains substring none
if ~strcmpi(spacedir_vecs{i}, 'none')
fprintf(['WARNING nhdr_nrrd_read: detected %s instead of ' ...
'expected none for axis %d of the per-axis field specifications "space directions:".\n' ...
' There should be no quotation marks, parentheses or any other characters, just plain none.\n'],...
spacedir_vecs{i}, i);
% Clean none vector specification
spacedir_vecs{i} = 'none';
end
else
% Axis-specific entry is a numerical vector
cnt_space_vectors = cnt_space_vectors + 1;
if cnt_space_vectors > internal_spacedimension
error(['%s:\n ''space directions'' field specifier: ' ...
'number of space vectors detected exceeds space (world)' ...
' dimension, which is equal to %d.'],...
fopen(fidr), internal_spacedimension)
end
btw_parentheses = spacedir_vecs{i}(1) == '(' ...
&& spacedir_vecs{i}(end) == ')';
axis_vector = regexprep(spacedir_vecs{i}, '[()]', ''); % stripped off all parens
if ~btw_parentheses
fprintf(['WARNING nhdr_nrrd_read: vector should be delimited ' ...
'by parentheses for axis %d of the per-axis field ' ...
'specifications "space directions:".\n' ...
' At least one missing parenthesis in ''%s''.\n'],...
i, spacedir_vecs{i})
end
% Clean up by forcing single enclosing parentheses:
spacedir_vecs{i} = ['(', axis_vector, ')'];
% Check vector and extract numerical data
vector_entries = strsplit(axis_vector, ',');
if length(vector_entries) ~= internal_spacedimension
error(['%s:\n vector for data axis %d (space axis %d) of the ' ...
'per-axis field specifications "space directions:" should ' ...
'contain %d entries corresponding to the space (or world) dimension ' ...
'specified in the "space" or "space dimension" field specification.' ...
' Found %d here.'],...
fopen(fidr), i, cnt_space_vectors, internal_spacedimension, ...
length(vector_entries))
end
for j = 1:length(vector_entries)
vector_entry = sscanf(vector_entries{j}, '%f');
if isempty(vector_entry)
error(['%s\n in field specification "space directions:", ' ...
'vector for data axis %d (space axis %d) too short. ' ...
'Detected %d entries instead of expected %d corresponding to ' ...
'space (world) dimension.'],...
fopen(fidr), i, cnt_space_vectors, j-1, internal_spacedimension)
end
SD_data(j, cnt_space_vectors) = vector_entry;
end
end
end
% Store array of cleaned up strings
headerInfo.spacedirections = spacedir_vecs; % cell array of strings, ideally of the form {'(1,0,0)' '(0,2,0.1)' 'none' '(0.1,0.1,1.1)'} if internal_spacedimension==3
% Extract numerical data more leniently:
SD_data_chk = strrep(space_dir_tmp, 'none', '' ); % ignore "none" entries
SD_data_chk = extractNumbersWithout(SD_data_chk, {'(',')',',', '"', ''''} ); % detects numbers up to first non numerical entry
if numel(SD_data_chk) ~= (internal_spacedimension)^2
error(['Expected ''space directions'' to specify a %d-by-%d matrix ' ...
'(%d elements in total) in agreement with world space dimension.' ...
' Found %d element(s) instead.\n'],...
internal_spacedimension, internal_spacedimension,...
(internal_spacedimension)^2, numel(SD_data_chk));
end
% Sanity check
if ~isequal(SD_data_chk(:), SD_data(:))
error(['%s:\n ''space directions'' field specifier: couldn''t' ...
' read space vectors. Please refer to the NRRD format definition.'],...
fopen(fidr))
end
headerInfo.spacedirections_matrix = SD_data;
% Correctness of field specification is checked below after
% whole header is parsed because we need to be sure that the
% "dimension" basic field specification was set
elseif foundKeyword('SPACE UNITS:', cs )
% Always optional, must come after space or space dimension
if define_orientation ~= 1
error('Field specification ''space units'' cannot be specified before ''space'' or ''space dimension''.')
end
space_units_tmp = strrep( cs(length('SPACE UNITS:')+1:end), 'none', ''); % ignore none entries
% FIXME: ideally, check that the sum of elements including none and
% " " matches headerInfo.dimension, i.e. the total dimension, as
% specified in the standard. Standard a bit unclear: should
% unknown units be specified with "???", "none" or "" ? (quotes
% seem to be required).
headerInfo.spaceunits = extract_spaceunits_list( space_units_tmp );
if length(headerInfo.spaceunits) ~= internal_spacedimension
error(['Expected ''space units'' to contain %d elements enclosed in double quotes' ...
' to match the ''space'' or ''space dimension'' field but found the following ' ...
'%d element(s):\n%s'],...
internal_spacedimension, length(headerInfo.spaceunits), sprintf('%s\t',headerInfo.spaceunits{:}))
end
elseif foundKeyword('SPACE ORIGIN:', cs )
% Always optional, must come after space or space dimension
assert(define_orientation==1, ...
sprintf(['%s: field ''space origin'' cannot be specified' ...
' before ''space'' or ''space dimension''.'], ...
fopen(fidr)));
iSO = extractNumbersWithout( cs(length('SPACE ORIGIN:')+1:end), {'(',')',','} );
assert(numel(iSO) == internal_spacedimension,...
sprintf(['%s: expected ''space origin'' to specify a ' ...
'%d-element vector to match the ''space'' or ' ...
'''space dimension'' field but found %d ' ...
'element(s).'],...
fopen(fidr),internal_spacedimension, numel(iSO)) );
headerInfo.spaceorigin = iSO(:);
elseif foundKeyword('MEASUREMENT FRAME:', cs )
% Always optional, must come after space or space dimension
assert(define_orientation==1,...
sprintf(['%s: field ''measurement frame'' cannot be ' ...
'specified before ''space'' or ''space dimension''.'],...
fopen(fidr)));
measframe_str = strrep( cs(length('MEASUREMENT FRAME:')+1:end), 'none', '');
iMF = extractNumbersWithout( measframe_str, {'(',')',','} ); % fails if non number entries are not previously removed
assert(numel(iMF) == (internal_spacedimension)^2,...
sprintf(['%s: expected ''measurement frame'' to specify a ' ...
'%d-by-%d matrix (%d total elements) but found %d ' ...
'element(s).'],...
fopen(fidr), internal_spacedimension, ...
internal_spacedimension, (internal_spacedimension)^2,...
numel(iMF)));
headerInfo.measurementframe = reshape(iMF(:), [internal_spacedimension, internal_spacedimension]);
% FIXME: this will gladly accept '1 0 0 0 1 0 0 0 1' instead of
% '(1,0,0) (0,1,0) (0,0,1)' for instance. But it should have
% length spacedimension according to standard so this is not too bad.
elseif foundKeyword('THICKNESSES:', cs )
sThicknesses = extractStringList( cs(length('THICKNESSES:')+1:end) ); % fixed bug with extractStringList where 2 entries are present
iThicknesses = [];
lenThicknesses = length( sThicknesses );
for iI=1:lenThicknesses
iThicknesses = [iThicknesses, str2double(sThicknesses{iI}) ];
end
headerInfo.thicknesses = iThicknesses;
elseif foundKeyword('CENTERINGS:', cs )
headerInfo.centerings = extractStringList( cs(length('CENTERINGS:')+1:end ) ); % fixed bug with extractStringList where 2 entries are present
elseif foundKeyword('LINE SKIP:', cs) || foundKeyword('LINESKIP', cs)
if foundKeyword('LINE SKIP:', cs)
headerInfo.lineskip = sscanf( cs( length('LINE SKIP:')+1:end ), '%d' );
else
headerInfo.lineskip = sscanf( cs( length('LINESKIP:')+1:end ), '%d' );
end
assert(headerInfo.lineskip >= 0,...
sprintf(['Field lineskip or line skip should be greater' ...
' than or equal to zero, detected %d.'], ...
headerInfo.lineskip));
elseif foundKeyword('BYTE SKIP:', cs) || foundKeyword('BYTESKIP:', cs)
if foundKeyword('BYTE SKIP:', cs)
headerInfo.byteskip = sscanf( cs( length('BYTE SKIP:')+1:end ), '%d' );
else
headerInfo.byteskip = sscanf( cs( length('BYTESKIP:')+1:end ), '%d' );
end
assert(headerInfo.byteskip >= -1, ...
sprintf(['Field byteskip or byte skip can only take ' ...
'non-negative integer values or -1, detected %d.'], ...
headerInfo.byteskip));
elseif foundKeyword('MODALITY', cs )
headerInfo.modality = strtrim( extractKeyValueString( cs(length('MODALITY')+1:end ) ) );
elseif foundKeyword('DWMRI_B-VALUE', cs )
headerInfo.bvalue = str2double( extractKeyValueString( cs(length('DWMRI_B-VALUE')+1:end ) ) );
elseif foundKeyword('DWMRI_GRADIENT_', cs )
[iGNr, dwiGradient] = extractGradient(cs(length('DWMRI_GRADIENT_')+1:end ));
headerInfo.gradients(iGNr+1,:) = dwiGradient;
% FIXME: make it more general than numbering limited to 4 digits?
elseif foundKeyword('DATA FILE:', cs ) || foundKeyword('DATAFILE:', cs)
% This tells us it is a detached header and ends it. 3 possibilities here
field_value = strtrim( cs(length('DATA FILE:')+1:end) );
if foundKeyword('DATAFILE:', cs)
field_value = strtrim( cs(length('DATAFILE:')+1:end) );
end
[filelist, LIST_mode, subdim] = extract_datafiles(field_value);
headerInfo.datafiles = filelist;
% In LIST mode, filelist is empty and is filled by reading the
% rest of the header
if LIST_mode
datafile_cnt = 0;
while ~feof(fidr)
cs = fgetl(fidr);
if isempty(cs)
break;
end
datafile_cnt = datafile_cnt + 1;
headerInfo.datafiles{datafile_cnt} = cs;
end
end
% We could technically break the loop here bc data file/datafile is
% supposed to close a header; however I think it does no harm to
% keep parsing (in LIST_mode, end of file reached anyways)
else
% see if we are dealing with a comment
csTmp = strtrim( cs );
if csTmp(1)~='#' && ~strcmp(cs(1:4),'NRRD')
fprintf('WARNING nhdr_nrrd_read: Could not parse input line: ''%s'' \n', cs );
% REVIEW: is it better to blindly write it to the output
% structure in order to be able to write the exact same file
% later on ?
end
end
end
% Check for required fields
% REVIEW: should this be checked only when the data is read? People
% might be interested in just parsing the header no matter what is
% in it...)
assert(isfield(headerInfo, 'sizes'), 'Missing required ''sizes'' field in header');
assert(isfield(headerInfo, 'dimension'), 'Missing required ''dimension'' field in header');
assert(isfield(headerInfo, 'encoding'), 'Missing required ''encoding'' field in header');
assert(isfield(headerInfo, 'type'), 'Missing required ''type'' field in header');
% Check other cross-field dependencies, etc.
% TODO: assert lengths of all
% fields specified, check that kinds and axis sizes match, etc.
if isfield(headerInfo, 'byteskip') && headerInfo.byteskip == -1
assert( strcmpi(headerInfo.encoding, 'raw'), ...
sprintf('byte skip value of -1 is only valid with raw encoding. See definition of NRRD File Format for more information.\n'));
end
matlabdatatype = nrrd_getMatlabDataType(headerInfo.type);
% endian field required if data defined on 2 bytes or more
if ~(any(strcmpi(headerInfo.encoding,{'txt', 'text', 'ascii'})) || any(strcmpi(matlabdatatype, {'int8', 'uint8'})))
assert(isfield(headerInfo, 'endian'), 'Missing required ''endian'' field in header');
else
if ~isfield(headerInfo,'endian')
headerInfo.endian = 'little'; % useless but harmless, just for code compatibility because endian field may be accessed later on while reading data
end
end
% Space and orientation information
if define_orientation
assert(isfield(headerInfo, 'spacedirections'), ...
['Missing field ''space directions'', required if either' ...
' ''space'' or ''space dimension'' is set.']);
if length(headerInfo.spacedirections) ~= headerInfo.dimension
fprintf(['WARNING nhdr_nrrd_read %s:\n Unexpected format found for ''space directions'' specifier.\n',...
' Expected none entries and vectors delimited by parentheses such as (1.2,0.2,0.25), ', ...
'the number of entries being equal to the dimension field specification.\n',...
' See definition of NRRD File Format for more information.\n'], fopen(fidr));
end
end
% TODO: add support for positive line skips
if isfield(headerInfo, 'lineskip') && headerInfo.lineskip > 0
assert(isfield(headerInfo, 'byteskip') && headerInfo.byteskip == -1, ...
sprintf(['lineskip option is currently not supported and can ' ...
'only be set to zero unless raw encoding is used and ' ...
'byte skip is set to -1 (which cancels the effect of ' ...
'lineskip altogether).']));
end
% TODO: add support for positive byte skips (see line skip)
if isfield(headerInfo, 'byteskip') && ~strcmpi(headerInfo.encoding,'raw')
assert( headerInfo.byteskip == 0, ...
sprintf('byte skip option with non raw encoding is currently not supported and can only be set to zero.\n'));
end
if isfield(headerInfo, 'byteskip') && strcmpi(headerInfo.encoding, 'raw')
assert( headerInfo.byteskip == -1, ...
sprintf('non-negative byte skip values with raw encoding are currently not supported; byte skip can only be set to -1.\n'));
end
catch me
% Clean up before raising error
fclose(fidr);
rethrow(me);
end
% Read the data. Detect if data is in the file or in a detached data file
% and check what type of detached file it is if need be.
if bReadData
N_data_tot = prod(headerInfo.sizes);
if isfield(headerInfo, 'datafiles')
% This is a detached header file
% We no longer need to read the header (closing it before reading
% all the detached data files is safer)
fclose(fidr);
% TODO: we currently only read slices along the slowest axis (i.e.,
% last coordinate)
if ~isempty(subdim) && subdim~=headerInfo.dimension
error(['(detached header): reading data from slices along axis other than the last' ...
' (i.e. slowest) one, is currently not supported.\n' ...
'Last argument [<subdim>] in ''data file'' or ''datafile'' field ' ...
'should be removed or set equal to %d, which is the detected' ...
' ''dimension'' field value, for now.'],...
headerInfo.dimension);
end
% Read data chunk by chunk from detached data files
N_data_files = length(headerInfo.datafiles);
assert(mod(N_data_tot, N_data_files)==0, ...
sprintf(['Number of detected data files (%d) does not divide total' ...
' number of values contained in data %d obtained from prod(sizes=[%s]).\n'],...
N_data_files, N_data_tot, sprintf('%d ',headerInfo.sizes)));
N_data_per_file = N_data_tot/N_data_files;
headerInfo.data = zeros(headerInfo.sizes, matlabdatatype); % specify right type of zeros, otherwise double by default
for i = 1:N_data_files
% Check type of detached data file
[~,fname_data,ext_data] = fileparts(headerInfo.datafiles{i}); % redundant because done above
data_ind = (i-1)*N_data_per_file+1:i*N_data_per_file;
if strcmpi(ext_data,'.nhdr')
error(['datafile %di/%d: nhdr file should not be used as ' ...
'detached data file.'], ...
i, length(headerInfo.datafiles));
elseif strcmpi(ext_data, '.nrrd')
% Detached NRRD file
bRead_Detached_Data = true;
tmp_struct = nhdr_nrrd_read(fullfile(mainFpath, ...
[fname_data, ext_data]), ...
bRead_Detached_Data); % recursive call
% TODO: check the rest of the structure for
% inconsistencies with metadata from header
assert(N_data_files==1 || headerInfo.dimension == tmp_struct.dimension +1, ...
sprintf(['Detached header %s: the number of dimensions in detached ' ...
'nrrd data file %d/%d (%s) should be one fewer than that of' ...
' the header file.\nDetected %d instead of %d.\nDifferent ' ...
'datafile dimensions as specified by the nrrd standard are' ...
' not supported as of now.\n'],...
fopen(fidr), i, N_data_files, headerInfo.datafiles{i},...
tmp_struct.dimension, headerInfo.dimension-1 ));
headerInfo.data(data_ind) = tmp_struct.data(:);
% Store detached data header for last detached file seen,
% assuming that all are the same. Do not store data though
% as it would be redundant.
tmp_struct = rmfield(tmp_struct, 'data');
headerInfo.detached_header = tmp_struct;
else
% e.g., detached .raw file
fid_data = fopen( fullfile(mainFpath, [fname_data, ext_data]), 'r');
if( fid_data < 1 )
error(['While reading detached header file %s:\ndetached ' ...
'data file number %d/%d (%s) could not be opened.'], ...
nrrdFileName, i, N_data_files, headerInfo.datafiles{i});
end
try
tmp_data = readData(fid_data, N_data_per_file, headerInfo.encoding, matlabdatatype);
fclose(fid_data);
catch me_detached
fclose(fid_data);
rethrow(me_detached);
end
tmp_data = adjustEndian(tmp_data, headerInfo.endian);
headerInfo.data(data_ind) = tmp_data(:);
end
end
else
% This is a NRRD standalone file: read data directly from it
try
headerInfo.data = readData(fidr, N_data_tot, headerInfo.encoding, matlabdatatype);
fclose(fidr);
catch me_detached
fclose(fidr);
rethrow(me_detached);
end
headerInfo.data = adjustEndian(headerInfo.data, headerInfo.endian);
headerInfo.data = reshape(headerInfo.data, headerInfo.sizes(:)'); % data into expected form. Transpose required by reshape function bc size vectors need to be row vectors
end
else
fclose(fidr);
end
end
% -------------------------------------------------------------------%
% ====================================================================
% --- Auxiliary functions -------------------------------------------%
% ====================================================================
function [iGNr, dwiGradient] = extractGradient( st )
% first get the gradient number
iGNr = str2num( st(1:4) ); % FIX numbering limited to 4 digits (synchronize with writer module)
% find where the assignment is
assgnLoc = strfind( st, ':=' );
if ( isempty(assgnLoc) )
dwiGradient = [];
return;
else
dwiGradient = sscanf( st(assgnLoc+2:end), '%f' );
end
end
% Return part of string after :=
function kvs = extractKeyValueString( st )
assgnLoc = strfind( st, ':=' );
if ( isempty(assgnLoc) )
kvs = [];
return;
else
kvs = st(assgnLoc(1)+2:end);
end
end
% Turn space-separated list into cell array of strings.
function sl = extractStringList( strList )
sl = strsplit(strtrim(strList)); % old Matlab file exchange version had a strange bug with lists of length 2
end
% Store in an array the list of numbers separated by the tokens listed in
% the withoutTokens cell array
function iNrs = extractNumbersWithout( inputString, withoutTokens )
auxStr = inputString;
for iI=1:length( withoutTokens )
auxStr = strrep( auxStr, withoutTokens{iI}, ' ' );
end
iNrs = sscanf( auxStr, '%f' );
end
% Return true if the string keyWord is the beginning of the content string
% cs (ignoring case)
function fk = foundKeyword( keyWord, cs )
lenKeyword = length( keyWord );
fk = (lenKeyword <= length(cs)) && strcmpi( cs(1:lenKeyword), keyWord);
end
% Return cell array of strings with space units in the form {mm, km, m}
% from input of the form "mm " "mm" " m" where undesired blank spaces may
% have crept in. The double quotes are removed in the process but will be
% added by the nhdr/nrrd writer function.
function su_ca = extract_spaceunits_list( fieldValue )
fv_trimmed = strtrim( fieldValue );
su_ca = strsplit(fv_trimmed, '"'); % units are delimited by double quotes
su_ca = su_ca(~ ( strcmp(su_ca, '') | strcmp(su_ca, ' ') ) ); % remove empty or blank space strings
for i = 1:length(su_ca)
su_ca{i} = strtrim( su_ca{i} );
end
end
% Extract data file list into a cell array from the datafile or data file
% field in a NHDR file, stripped off its leading and trailing spaces;
% LIST_mode is set to one if a list of files closes the header;
% subdim is empty by default and may be set through either of these field
% specifications:
% data file: <format> <min> <max> <step> [<subdim>]
% data file: LIST [<subdim>]
function [filelist, LIST_mode, subdim] = extract_datafiles( field_string)
field_string = strtrim(field_string);
filelist = {};
LIST_mode = 0;
subdim = [];
if length(field_string)>= 4 && strcmpi(field_string(1:4), 'LIST')
% LIST mode
subdim = sscanf(field_string(5:end), '%d'); % assert 1<=dimensions
LIST_mode = 1;
return;
else
% single detached data file or multiple files written in concise form, non LIST mode
str_lst = strsplit(field_string); % without delimiters specified, splits at any sequence in the set {' ','\f','\n','\r','\t','\v'}
if length(str_lst) == 1
% Single detached data file:
filelist{1} = str_lst{1};
% subdim does not make sense here since all of the data is
% contained in the deatached data file
elseif any(length(str_lst)==[0, 2, 3]) || length(str_lst)>5
% Invalid cases (2, 3 or striclty more than 5 entries)
error(['error in ''data list'' (or ''data list'') field, in non' ...
' ''LIST'' mode:\nexpected ''<filename>'' or ''<format> ' ...
'<min> <max> <step> [<subdim>]'' but found instead ''%s'' ' ...
'(%d elements instead of 1, 4 or 5).'],...
sprintf('%s ',str_lst{:}), length(str_lst));
else
% Multiple detached data files with filenames including integral
% values generated by min:step:max
str_format = str_lst{1};
id_min = sscanf(str_lst{2}, '%d');
id_max = sscanf(str_lst{3}, '%d');
step = sscanf(str_lst{4}, '%d');
% check format and indexing values
assert(step~=0,...
sprintf(['detached data files with names specified by ' ...
'sprintf()-like format: step should be strictly positive' ...
' or negative, not zero.']));
if step < 0
assert(id_min >= id_max, ...
sprintf(['detached data files with names specified by ' ...
'sprintf()-like format: when step is <0, min ' ...
'should be larger than or equal to max, here we ' ...
'found min=%d < max=%d.'],id_min, id_max));
else
assert(id_min <= id_max,...
sprintf(['detached data files with names specified by ' ...
'sprintf()-like format: : when step is >0, min ' ...
'should be smaller than or equal to max, here we ' ...
'found min=%d > max=%d.'],id_min, id_max));
end
% populate data file list
fileIDs = id_min:step:id_max;
filelist = cell(length(fileIDs), 1);
for i = 1:length(fileIDs) % does not necessarily incude <max>, it cannot be larger than <max>
expanded_fname = sprintf(str_format, fileIDs(i));
filelist{i} = expanded_fname;
end
if length(str_lst) == 5
subdim = sscanf(str_lst{5}, '%d');
end
end
end
end
% Read data in a column vector from open file with pointer fidIn at the
% beginning of the data to read. The other arguments specify the number of
% elements to read (determined in main body of reader function above), data
% file encoding (specified in NHDR/NRRD header) and the Matlab type the
% data should be converted to (determined in main body of reader function
% above).
function data = readData(fidIn, Nelems, meta_encoding, matlabdatatype)
switch (meta_encoding)
case {'raw'}
% This implicitely assumes that byteskip was set to -1 in raw
% encoding
data_bytes = (Nelems)*sizeOf(matlabdatatype); % bytes of useful info to return in data array
fseek(fidIn,0,'eof'); % set position pointer to end of file
tot_bytes_file = ftell(fidIn); % current location of the position pointer (I believe as a byte offset from beginning of file)
fseek(fidIn,tot_bytes_file-data_bytes,'bof'); % make sure position pointer is right before useful data starts
% Just read binary file (original version of the code)
data = fread(fidIn, inf, [matlabdatatype '=>' matlabdatatype]);
case {'gzip', 'gz'}
tmp = fread(fidIn, Inf, 'uint8=>uint8'); % byte per byte
tmpBase = tempname(pwd);
tmpFile = [tmpBase '.gz'];
fidTmp = fopen(tmpFile, 'wb'); % this creates tmpFile
assert(fidTmp > 3, 'Could not open temporary file for GZIP decompression')
try
fwrite(fidTmp, tmp, 'uint8');
catch me
fclose(fidTmp);
delete(tmpFile);
rethrow(me);
end
fclose(fidTmp);
try
gunzip(tmpFile); % this creates tmpBase (tmpFile without the .gz)
catch me
delete(tmpFile);
rethrow(me);
end
delete(tmpFile);
fidTmp = fopen(tmpBase, 'rb');
% cleaner = onCleanup(@() fclose(fidTmp));
try
data = readData(fidTmp, Nelems, 'raw', matlabdatatype); % recursive call
catch me
fclose(fidTmp);
delete(tmpBase);
rethrow(me);
end
fclose(fidTmp);
delete(tmpBase);
case {'txt', 'text', 'ascii'}
data = fscanf(fidIn, '%f');
data = cast(data, matlabdatatype);
otherwise
assert(false, 'Unsupported encoding')
end
% Check number of read elements (sometimes there is an offest of about 1)
assert(Nelems == numel(data),...
sprintf(['Error reading binary content of %s: detected %d elements' ...
' instead of %d announced in header.'],...
fopen(fidIn), numel(data), Nelems));
end
% Switch byte ordering according to endianness specified in nhdr/nrrd file.
% The metadata should just contain a field 'endian' set to 'little' or
% 'big'.
function data = adjustEndian(data, meta_endian)
[~,~,endian] = computer();
needToSwap = (isequal(endian, 'B') && isequal(lower(meta_endian), 'little')) || ...
(isequal(endian, 'L') && isequal(lower(meta_endian), 'big'));
if (needToSwap)
data = swapbytes(data);
end
end
% Size of Matlab's numeric classes in bytes
% See: https://stackoverflow.com/questions/16104423/size-of-a-numeric-class
function numBytes = sizeOf(dataClass)
% create temporary variable of data type specified
eval(['var = ' dataClass '(0);']);
% Use the functional form of whos, and get the number of bytes
W = whos('var');
numBytes = W.bytes;
end