-
Notifications
You must be signed in to change notification settings - Fork 19
/
cifti_write.m
193 lines (190 loc) · 10.1 KB
/
cifti_write.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
function cifti_write(cifti, filename, varargin)
%function cifti_write(cifti, filename, option pairs...)
% Write a cifti file.
%
% Specifying "..., 'keepmetadata', true" leaves the file-level metadata as-is:
% if false or not specified, the 'Provenance' metadata value is moved to
% 'ParentProvenance', provenance keys are given generic values, and all other
% file-level metadata is removed.
%
% Specifying "..., 'disableprovenance', true" removes all file-level metadata
% before writing. If the 'keepmetadata' option is true, this option has no
% effect, but a warning is issued if both options are true.
%
% Example usage:
%
% >> cifti = cifti_read('91282_Greyordinates.dscalar.nii');
% >> cifti.cdata = outdata;
% >> cifti.diminfo{2} = cifti_diminfo_make_scalars(size(outdata, 2));
% >> cifti_write(cifti, 'ciftiout.dscalar.nii');
libversion = '2.2.2';
options = myargparse(varargin, {'stacklevel', 'disableprovenance', 'keepmetadata'}); %stacklevel is an implementation detail, don't add to help
if isempty(options.stacklevel) %stacklevel is so that so it doesn't get "ciftisave" all the time
options.stacklevel = 2;
end
options.keepmetadata = argtobool(options.keepmetadata, 'keepmetadata');
options.disableprovenance = argtobool(options.disableprovenance, 'disableprovenance');
if options.keepmetadata && options.disableprovenance
warning('both "keepmetadata" and "disableprovenance" are true, ignoring "disableprovenance"');
end
%dims_m = size(cifti.cdata); %NO: matlab can't represent [x, y, 1] matrices, always changes them to [x, y] on the fly
dims_m = [];
for i = 1:length(cifti.diminfo)
cifti.diminfo{i}.length = cifti_diminfo_length(cifti.diminfo{i}); %set diminfo .length from diminfo contents so it is synchronized
dims_m(i) = cifti.diminfo{i}.length; %#ok<AGROW>
end
sanity_check_cdata(cifti);
dims_c = dims_m([2 1 3:length(dims_m)]); %ciftiopen convention, first matlab index is down
if ~options.keepmetadata
if ~options.disableprovenance
stack = dbstack;
if options.stacklevel > length(stack)
newprov = 'written from matlab/octave prompt or script';
else
newprov = ['written from function ' stack(options.stacklevel).name ' in file ' stack(options.stacklevel).file];
end
prov = cifti_metadata_get(cifti.metadata, 'Provenance');
cifti.metadata = struct('key', {'Provenance', 'ProgramProvenance'}, 'value', {newprov, ['cifti_write.m ' libversion]});
if ~isempty(prov)
cifti.metadata = cifti_metadata_set(cifti.metadata, 'ParentProvenance', prov);
end
else
cifti.metadata = struct('key', {}, 'value', {});
end
end
xmlbytes = cifti_write_xml(cifti, true);
header = make_nifti2_hdr();
extension = struct('ecode', 32, 'edata', xmlbytes); %header writing function will pad the extensions with nulls
header.extensions = extension; %don't need concatenation for only one nifti extension
header.datatype = 16;
header.bitpix = 32;
header.dim(6:(6 + length(dims_c) - 1)) = dims_c;
header.dim(1) = length(dims_c) + 4;
[header.intent_code, header.intent_name] = cifti_intent_code(cifti, filename);
[fid, header, cleanupObj] = write_nifti2_hdr(header, filename); %#ok<ASGLU> %header writing also computes vox_offset for us
%the fseek probably isn't needed during writing, but to be safe
if(fseek(fid, header.vox_offset, 'bof') ~= 0)
error(['failed to seek to start data writing in file ' filename]);
end
%we need to swap the first 2 dims, and 'permute' effectively makes a copy of its input, so write large files in chunks instead
%FIXME: if we allow setting nifti scale/intercept, that needs to be added to this code
max_elems = 128 * 1024 * 1024 / 4; %assuming float32, use only 128MiB extra memory when writing (or the size of a row, if that manages to be larger)
if numel(cifti.cdata) <= max_elems
%file is small, use simple whole-array writing code
fwrite_excepting(fid, permute(cifti.cdata, [2 1 3:length(size(cifti.cdata))]), 'float32');
else
max_rows = max(1, floor(max_elems / dims_c(1)));
if max_rows < dims_c(2) %less than a plane at a time, don't write cross-plane
num_passes = ceil(dims_c(2) / max_rows); %even out the passes
chunk_rows = ceil(prod(dims_c(2:end)) / num_passes);
total_planes = prod(dims_c(3:end));
for plane = 1:total_planes
for chunkstart = 1:chunk_rows:dims_c(2)
%since it is part of one plane, technically matlab will return a 2D matrix, and we could just use transpose...
fwrite_excepting(fid, permute(cifti.cdata(chunkstart:min(dims_c(2), chunkstart + chunk_rows - 1), :, plane), [2 1 3]), 'float32');
end
end
else
max_planes = max(1, floor(max_rows / dims_c(2))); %just in case the division does something dumb
total_planes = prod(dims_c(3:end)); %flatten all dimensions 3+
num_passes = ceil(total_planes / max_planes);
chunk_planes = ceil(total_planes / num_passes);
for chunkstart = 1:chunk_planes:total_planes
fwrite_excepting(fid, permute(cifti.cdata(:, :, chunkstart:min(total_planes, chunkstart + chunk_planes - 1)), [2 1 3]), 'float32');
end
end
end
end
function [code, string] = cifti_intent_code(cifti, filename)
code = 3000;
string = 'ConnUnknown';
expectext = '';
explain = '';
numdims = length(cifti.diminfo); %this will be at least 2, we call this after checking dims (and after writing xml, so the types are all good)
switch cifti.diminfo{1}.type %NOTE: column
case 'dense'
switch cifti.diminfo{2}.type
case 'dense'
code = 3001; string = 'ConnDense'; expectext = '.dconn.nii'; explain = 'dense by dense';
case 'series'
code = 3002; string = 'ConnDenseSeries'; expectext = '.dtseries.nii'; explain = 'series by dense'; %order by cifti convention rather than matlab?
case 'scalars'
code = 3006; string = 'ConnDenseScalar'; expectext = '.dscalar.nii'; explain = 'scalars by dense';
case 'labels'
code = 3007; string = 'ConnDenseLabel'; expectext = '.dlabel.nii'; explain = 'labels by dense';
case 'parcels'
code = 3010; string = 'ConnDenseParcel'; expectext = '.dpconn.nii'; explain = 'parcels by dense';
end
case 'parcels'
switch numdims
case 3
if strcmp(cifti.diminfo{2}.type, 'parcels')
switch cifti.diminfo{3}.type
case 'series'
code = 3011; string = 'ConnPPSr'; expectext = '.pconnseries.nii'; explain = 'parcels by parcels by series';
case 'scalars'
code = 3012; string = 'ConnPPSc'; expectext = '.pconnscalar.nii'; explain = 'parcels by parcels by scalar';
end
end
case 2
switch cifti.diminfo{2}.type
case 'parcels'
code = 3003; string = 'ConnParcels'; expectext = '.pconn.nii'; explain = 'parcels by parcels';
case 'series' %these two are max length to have a null terminator in the field
code = 3004; string = 'ConnParcelSries'; expectext = '.ptseries.nii'; explain = 'series by parcels';
case 'scalars'
code = 3008; string = 'ConnParcelScalr'; expectext = '.pscalar.nii'; explain = 'scalars by parcels';
case 'dense'
code = 3009; string = 'ConnParcelDense'; expectext = '.pdconn.nii'; explain = 'dense by parcels';
end
end
end
if isempty(expectext)
periods = find(filename == '.', 2, 'last');
if length(periods) < 2 || length(filename) < 4 || ~myendswith(filename, '.nii')
warning(['cifti file with nonstandard mapping combination "' filename '" should be saved ending in .<something>.nii']);
else
problem = true;
switch filename(periods(1):end)
case '.dconn.nii'
case '.dtseries.nii'
case '.dscalar.nii'
case '.dlabel.nii'
case '.dpconn.nii'
case '.pconnseries.nii'
case '.pconnscalar.nii'
case '.pconn.nii'
case '.ptseries.nii'
case '.pscalar.nii'
case '.pdconn.nii'
case '.dfan.nii'
case '.fiberTemp.nii'
otherwise
problem = false;
end
if problem
warning(['cifti file with nonstandard mapping combination "' filename '" should Not be saved using an already-used cifti extension, please choose a different, reasonable cifti extension of the form .<something>.nii']);
end
end
else
if ~myendswith(filename, expectext)
if ~strcmp(expectext, '.dscalar.nii')
warning([explain ' cifti file "' filename '" should be saved ending in ' expectext]);
else
if ~(myendswith(filename, '.dfan.nii') || myendswith(filename, '.fiberTEMP.nii'))
warning([explain ' cifti file "' filename '" should be saved ending in ' expectext]);
end
end
end
end
end
function output = argtobool(input, argname)
switch input
case {0, false, '', '0', 'no', 'false'} %empty string defaults to false
output = false;
case {1, true, '1', 'yes', 'true'}
output = true;
otherwise
error(['unrecognized value for option "' argname '", please use 0/1, true/false, yes/no']);
end
end