-
Notifications
You must be signed in to change notification settings - Fork 260
json header extension
A draft specification of the JSON header for Nibabel.
JSON, as y'all know, encodes strings, numbers, objects and arrays, An object is like a Python dict, with strings as keys, and an array is like a Python list.
In what follows, I will build dicts and lists corresponding to the objects and
arrays of the JSON header. In each case, the json.dumps of the given Python
object gives the corresponding JSON string.
I'll use the term field to refer to a (key, value) pair from a Python dict / JSON object.
We specify image axes by name in the header, and give the correspondence of the
names to the image array axes by the order of the names. This is the
axis_names field at the top level of the header.
If the user transposes or otherwise reorders the axes of the data array, the
header should change only in the ordering of the axis names in
axis_names. Call this the "axis transpose" principle.
The JSON header should make sense as a key, value pair store for DICOM fields using a standard way of selecting DICOM fields -- the simple DICOM principle.
-
JSON-LD - provides a way of using json that can be mapped into the Resource Description Framework (RDF). It is highly recommended to take a look at the RDF Primer to get a sense of why we might want to use JSON-LD/RDF, but essentially it boils down to a couple points:
-
JSON keys are turned into URIs
-
URIs can dereference to a Web URL with additional documentation, such as a definition, a pretty label (e.g.,
nipy_header_versionhas_label"NIPY Header Version"), etc. -
The URI link to documentation makes the meaning of your JSON keys explicit, in a machine readable way (i.e., the json key becomes a "resource" on the Web that avoids name clashes)
-
JSON-LD/RDF has a full query language called SPARQL and a python library called RDFLib that acts as a parser, serializer, database, and query engine.
-
In the example below, the
@contextsection provides the namespace prefixdcmas a placeholder for the URLhttp://neurolex.org/wiki/Category:, thusdcm:Echo_Timedereferences to http://neurolex.org/wiki/Category:Echo_Time where additional documentation is provided:{ "@context": { "dcm": "http://neurolex.org/wiki/Category:#" }, "dcm:Echo_Time": 45, "dcm:Repetition_Time": 2, }
-
>>> hdr = dict(nipy_header_version='1.0')We chose the name "nipy_header_version" in the hope that this would not often occur in an unrelated JSON file.
- First version will be "1.0".
- Versioning will use Semantic Versioning of form
major.minor[.pico[-extra]]wheremajor,minor,picoare all integers,extramay be a string, and bothpicoandextraare optional. Header versions with the samemajorvalue are forwards compatible -- that is, a reader that can read a header with a particular major version should be able to read any header with that major version. Specifically, any changes to the header format within major version number should allow older readers of that major version to read the header correctly, but can expand on the information in the header, so that older readers can safely ignore new information in the header. - All fields other than
nipy_header_versionare optional.hdrabove is therefore the minimal valid header.
The base level header will usually also have image metadata fields giving information about the whole image. A field is an "image metadata field" if it is defined at the top level of the header. For example:
>>> hdr = dict(nipy_header_version='1.0',
... Manufacturer="SIEMENS")All image metadata fields are optional.
As for all keys in this standard, IM (Image Metadata) keys are case sensitive. IM keys that begin with a capital letter must be from the DICOM data dictionary standard short names (DICOM keyword). Call these "DICOM IM keys". This is to conform to the simple DICOM principle.
Keys beginning with "extended" will be read and written, but not further processed by a header reader / writer. If you want to put extra fields into the header that are outside this standard you could use a dict / object of form:
>>> hdr = dict(nipy_header_version='1.0',
... extended=dict(my_field1=0.1, my_field2='a string'))or:
>>> hdr = dict(nipy_header_version='1.0',
... extended_mysoft=dict(mysoft_one='expensive', mysoft_two=1000))Values for DICOM IM keys are constrained by the DICOM standard. This standard constrains values for ("nipy_header_version", "axis_names", "axis_metadata"). Other values have no constraint.
Questions:
-
Should all DICOM values be allowed?
-
Should DICOM values be allowed at this level that in fact refer to a particular axis, and therefore might go in the
axis_metadataelements? -
How should we relate the DICOM standard values to JSON? For example, how should we store dates and times? One option would be to use the new DICOM JSON encoding for DICOM values, but omitting the tag and value representation (VR). For example, the DICOM JSON spec has:
"00080070": { "vr": "LO", "Value": [ "SIEMENS" ] },but we might prefer:
"Manufacturer": "SIEMENS"
Using the DICOM data dictionary we can reconstruct the necessary tag and VR, so our version is lossless if the DICOM keyword exists in the DICOM data dictionary. Of course this may well not be true for private tags, or if the keyword comes from a DICOM dictionary that is later than the one we are using to look up the keyword. For the latter, we could make sure we're always using the latest dictionary. For the private tags, we might want to recode these in any case, maybe using our own dictionary. Maybe it is unlikely we will want to reconstruct the private tags of a DICOM file from the JSON. Comments welcome.
axis_names is a list of strings corresponding to the axes of the image data
to which the header refers.
>>> hdr = dict(nipy_header_version='1.0',
... axis_names=["frequency", "phase", "slice", "time"])The names must be valid Python identifiers (should not begin with a digit, nor contain spaces etc).
There must be the same number of names as axes in the image to which the header refers. For example, the header above is valid for a 4D image but invalid for a 3D or 5D image.
The names appear in fastest-slowest order in which the image data is stored on
disk. The first name in axis_names corresponds to the axis over which
the data on disk varies fastest, and the last corresponds to the axis over which
the data varies slowest.
For a Nifti image, nibabel (and nipy) will create an image where the axes have
this same fastest to slowest ordering in memory. For example, let's say the
read image is called img. img has shape (4, 5, 6, 10), and a 2-byte
datatype such as int16. In the case of the nifti default fastest-slowest ordered
array, the distance in memory between img[0, 0, 0, 0] and img[1, 0, 0,
0] is 2 bytes, and the distance between img[0, 0, 0, 0] and img[0, 0, 0,
1] is 4 * 5 * 6 * 2 = 240 bytes. The names in axis_names will then refer
to the first, second, third and fourth axes respectively. In the example above,
"frequency" is the first axis and "time" is the last.
axis_names is optional only if axis_metadata is empty or absent.
Otherwise, the set() of axis_names must be a superset of the union of
all axis names specified in the applies_to fields of axis_metadata
elements.
axis_metadata is a list of axis metadata elements.
Each axis metadata element in the axis_metadata list gives data that
applies to a particular axis, or combination of axes. axis_metadata can
be empty:
>>> hdr['axis_metadata'] = []We prefer you delete this section if it is empty, to avoid clutter, but hey, mi casa, su casa.
An axis metadata element must contain a field applies_to, with a value that
is a list that contains one or more values from axis_names. From the above
example, the following would be valid axis metadata elements:
>>> hdr = dict(nipy_header_version='1.0',
... axis_names = ["frequency", "phase", "slice", "time"],
... axis_metadata = [
... dict(applies_to = ['time']),
... dict(applies_to = ['slice']),
... dict(applies_to = ['slice', 'time']),
... ])The elements will usually also have axis metadata fields. For example:
>>> element = dict(applies_to=['time'],
... TimeSliceVector=[0, 2, 4])As for image metadata keys, keys that begin with a capital letter are DICOM standard keywords.
A single axis name for applies_to specifies that any axis metadata values in
the element apply to the named axis.
In this case, axis metadata values may be:
- a scalar. The value applies to every point along the corresponding image axis.
- a vector of length N (where N is the length of the corresponding image axis).
Value
$v_i$ in the vector$v$ corresponds to the image slice at point$i$ on the corresponding axis. - an array of shape (N, ...) or (1, ...) where "..." can be any further shape. The (N, ...) case gives N vectors or arrays with one (vector, array) corresponding to each point in the image axis, The (1, ...) case gives a single vector or array corresponding to every point on the image axis.
More than one axis name for applies_to specifies that any values in the
element apply to the combination of the given axes.
In the case of more than one axis for applies_to, the axis metadata values
apply to the Cartesian product of the image axis values. For example, if the
values of applies_to == ['slice', 'time'], and the slice and time axes
in the array are lengths (6, 10) respectively, then the values apply to all
combinations of the 6 possible values for slice indices and the 10 possible
values for the time indices (ie apply to all 6x10=60 values). The axis metadata
values in this case can be:
- a scalar. The value applies to every combination of (slice, time)
- an array of shape (S, T) (where S is the length of the slice axis and T is the
length of the time axis). Value
$a_{i,j}$ in the array$a$ corresponds to the image slice at point$i$ on the slice axis and$j$ on the time axis. - an array of shape (S, T, ...) or (1, 1, ...) where "..." can be any further shape. The (S, T, ...) case gives N vectors or arrays with one vector / array corresponding to each combination of slice, time points in the image, The (1, 1, ...) case gives a single vector / array corresponding to every (slice, time) point in the image.
- In the spirit of numpy array broadcasting, we also allow a value array for (slice, time) to be of shape (S, 1, ...) or (1, T, ...) where the arrays give (respectively): one vector / array for each point of the slice axis, but applying to every value of the time axis, or; one vector / array for each point of the time axis but applying to every point of the slice axis. Of course these may better go as axis metadata for the slice and time axes respectively. See :ref:`q_vector` for an example where this kind of specification can be useful.
In general, for a given value applies_to, we can take the corresponding axis
lengths:
>>> shape_of_image = [4, 5, 6, 10]
>>> image_names = ['frequency', 'phase', 'slice', 'time']
>>> applies_to = ['slice', 'time']
>>> axis_indices = [image_names.index(name) for name in applies_to]
>>> axis_lengths = [shape_of_image[i] for i in axis_indices]
>>> axis_lengths
[6, 10]The axis metadata value can therefore be of shape:
- () (a scalar) (a scalar value for every combination of points)
-
axis_lengths(a scalar value for each combination of points) -
[1 for i in axis_lengths] + any_other_list(a single array or vector for every combination of points, where the shape of the array or vector is given byany_other_list) -
axis_lengths + any_other_list(an array or vector corresponding to each combination of points, where the shape of the array or vector is given byany_other_list) -
axis_lengths_or_1 + any_other_list-- whereaxis_lengths_or_1is a vector that, at each positioni, contains either 1 oraxis_length[i]. The array or vector implied by theany_other_listindices applies to every point whereaxis_lengths_or_1[i] == 1for axisi, and each point along axisiwhereaxis_lengths_or_1[i] == axis_lengths[i]. Obviously the case of[1 for i in axis_lengths] + any_other_listis a specific case in this category.
We define an axis metadata field q_vector which gives the q vector
corresponding to the diffusion gradients applied.
The q_vector should apply to (applies_to) four axes, of which three are
the 3D spatial dimensions, and the fourth corresponds to image volume.
For example:
>>> import numpy as np
>>> element = dict(applies_to=['frequency', 'phase', 'slice', 'time'],
... q_vector = [[[[
... [0, 0, 0],
... [1000, 0, 0],
... [0, 1000, 0],
... [0, 0, 1000],
... [0, 0, 0],
... [1000, 0, 0],
... [0, 1000, 0],
... [0, 0, 1000],
... [0, 0, 0],
... [1000, 0, 0]
... ]]]])
>>> np.array(element['q_vector']).shape
(1, 1, 1, 10, 3)We specify the spatial axes in applies_to because of the axis transpose
principle -- we need to know which axes the vector directions correspond to, to
be robust to the case of an image transpose.
An individual (3,) vector is the unit vector expressing the direction of the gradient, multiplied by the scalar b value of the gradient. In the example, there are three b == 0 scans (corresponding to volumes 0, 4, 8), with the rest having b value of 1000.
The first value corresponds to the direction along the first named image axis ('frequency'), the second value to direction along the second named axis ('phase'), and the third to direction along the 'slice' axis.
Note that the q_vector is always specified in the axes of the image. This is
the same convention as FSL uses for its bvals and bvecs files.
This gives a list of times of acquisition of each slice in milliseconds relative to the beginning of the volume.
Applies to an image axis representing slices.
We use "slice" as the axis name here, but any name is valid.
>>> element = dict(applies_to=['slice'],
... slice_acquisition_times=[0, 20, 40, 60, 80, 100])NIfTI 1 and 2 can encode some slice acquisition times using a somewhat complicated scheme, but they cannot - for example - encode multi-slice acquisitions, and NIfTI slice time encoding is rarely set.
I propose that slice_acquisition_times always override the NIfTI
specification when this value is specified in JSON.
This is a list of times of acquisition of each volume in milliseconds relative to the beginning of the acquisition of the run.
The field applies to an image axis representing volumes.
These values can be useful for recording runs with missing or otherwise not-continuous time data.
We use "time" as the axis name, but any name is valid.
>>> element = dict(applies_to=['time'],
... volume_acquisition_times=[0, 120, 240, 480, 600])Milliseconds is a reasonable choice for units because a Python integer can decode / encode any integer number in the JSON correctly, a signed 32-bit int can encode to around 6000 hours, and a 32-bit float can encode to 23 hours without loss of precision.
So far we are allowing any axis to be a slice or volume axis, but it might be nice to check. One way of doing this is:
>>> element = dict(applies_to=['mytime'],
... axis_meanings=["volume", "time"],
... volume_acquisition_times=[0, 120, 240, 480, 600])
>>> element = dict(applies_to=['myslice'],
... axis_meanings=["slice"],
... slice_acquisition_times=[0, 20, 40, 60, 80, 100])In this case we can assert that slice_acquisition_times applies to an axis
with meanings that include "slice" and that volume_acquisition_times
applies to an axis with meaning "volume". For example:
>>> # Should raise an error on reading full JSON
>>> element = dict(applies_to=['myslice'],
... axis_meanings=["volume"],
... slice_acquisition_times=[0, 20, 40, 60, 80, 100])This might also help for the situation where there is more than one frequency axis:
>>> hdr = dict(nipy_header_version='1.0',
... axis_names = ["frequency1", "frequency2", "slice", "time"],
... axis_metadata = [
... dict(applies_to = ["frequency1"],
... axis_meanings = ["frequency"]),
... dict(applies_to = ["frequency2"],
... axis_meanings = ["frequency"]),
... dict(applies_to = ['slice'],
... axis_meanings = ["slice"]),
... dict(applies_to = ['time'],
... axis_meanings = ["time", "volume"]),
... ])We can also check that space axes really are space axes:
>>> hdr = dict(nipy_header_version='1.0',
... axis_names = ["frequency", "phase", "slice", "time"],
... axis_metadata = [
... dict(applies_to = ["frequency"],
... axis_meanings = ["frequency", "space"]),
... dict(applies_to = ["phase"],
... axis_meanings = ["phase", "space"]),
... dict(applies_to = ["slice"],
... axis_meanings = ["slice", "space"]),
... dict(applies_to = ["time"],
... axis_meanings = ["time", "volume"]),
... dict(applies_to=["frequency", "phase", "slice", "time"],
... q_vector = [[[[
... [0, 0, 0],
... [1000, 0, 0],
... ]]]])
... ])For the q_vector field, we can check that all of the applies_to axes
("frequency", "phase", "slice", "time").
When doing motion correction on a 4D image, we calculate the required affine transformation from |--| say |--| the second image to the first image; the third image to the first image; etc. If there are N volumes in the 4D image, we would need to store N-1 affine transformations. If we have registered to the mean volume of the volume series instead of one of the volumes in the volume series, then we need to store all N transforms.
We often want to store this set of required transformations with the image,
but NIfTI does not allow us to do that. SPM therefore stores these transforms
in a separate MATLAB-format .mat file. We currently don't read these
transformations because we have no API in nibabel to present or store multiple
affines.
Assume the 4D volume has T time points (volumes).
There are two ways we could implement the multi-affines. The first would be to have (1 x 1 x 1 x T x 3 x 4) array of affines, with one for each volume / time point. This would use the same general idea as the q_vector field:
>>> element = dict(applies_to=['frequency', 'phase', 'slice', 'time'],
... multi_affine = [[[[
... [[ 2.86, -0.7 , 0.83, -80.01],
... [ 0.71, 2.91, 0.01, -114.59],
... [ -0.54, 0.13, 4.42, -54.34]],
... [[ 2.87, -0.38, 1.19, -92.77],
... [ 0.31, 2.97, 0.45, -110.87],
... [ -0.82, -0.2 , 4.32, -33.89]],
... [[ 2.97, -0.39, 0.31, -78.95],
... [ 0.33, 2.9 , 1.06, -116.99],
... [ -0.29, -0.68, 4.36, -36.41]],
... [[ 2.93, -0.5 , 0.61, -78.02],
... [ 0.4 , 2.9 , 0.99, -118.9 ],
... [ -0.5 , -0.59, 4.35, -33.61]],
... [[ 2.95, -0.44, 0.49, -77.86],
... [ 0.3 , 2.78, 1.62, -125.83],
... [ -0.46, -1.03, 4.17, -21.66]]
... ]]]])
>>> np.array(element['multi_affine']).shape
(1, 1, 1, 5, 3, 4)Another option would be to partially follow NRRD in giving the column vectors from the affine to the axis to which they apply, and split the translation into a separate offset vector:
>>> hdr = dict(nipy_header_version='1.0',
... axis_names = ["frequency", "phase", "slice", "time"],
... axis_metadata = [
... dict(applies_to=['frequency', 'time'],
... output_vector = [[
... [ 2.86, 0.71, -0.54],
... [ 2.87, 0.31, -0.82],
... [ 2.97, 0.33, -0.29],
... [ 2.93, 0.4 , -0.5 ],
... [ 2.95, 0.3 , -0.46],
... ]]),
... dict(applies_to=['phase', 'time'],
... output_vector = [[
... [ -0.7 , 2.91, 0.13],
... [ -0.38, 2.97, -0.2 ],
... [ -0.39, 2.9 , -0.68],
... [ -0.5 , 2.9 , -0.59],
... [ -0.44, 2.78, -1.03],
... ]]),
... dict(applies_to=['slice', 'time'],
... output_vector = [[
... [ 0.83, 0.01, 4.42],
... [ 1.19, 0.45, 4.32],
... [ 0.31, 1.06, 4.36],
... [ 0.61, 0.99, 4.35],
... [ 0.49, 1.62, 4.17],
... ]]),
... dict(applies_to=['time'],
... output_offset = [
... [ -80.01, -114.59, -54.34],
... [ -92.77, -110.87, -33.89],
... [ -78.95, -116.99, -36.41],
... [ -78.02, -118.9, -33.61],
... [ -77.86, -125.83, -21.66],
... ])],
... )
>>> np.array(hdr['axis_metadata'][0]['output_vector']).shape
(1, 5, 3)
>>> np.array(hdr['axis_metadata'][1]['output_vector']).shape
(1, 5, 3)
>>> np.array(hdr['axis_metadata'][2]['output_vector']).shape
(1, 5, 3)
>>> np.array(hdr['axis_metadata'][3]['output_offset']).shape
(5, 3)