MVloader

MVloader is meant to be a tiny helper to load and save medical volumetric data (therefore MV) aka. image volumes, i.e. three-dimensional medical images (DICOM, NIfTI, or NRRD), in Python 3.5+. It is also meant to simplify dealing with their different anatomical world coordinate systems.

Table of Contents

• Installing MVloader
• Motivation: Voxel Indices, World Coordinates, and Patient Anatomy
• The Volume Class
• Loading Images
• Saving Images
• Acknowledgments
• History

Installing MVloader

MVloader can be installed using pip or pip3, depending on your system, from its github repository (requiring git being installed on your system in either case):

pip install git+https://github.com/spezold/mvloader.git
pip3 install git+https://github.com/spezold/mvloader.git

This should also work inside conda.

MVloader's only remaining dependencies are the underlying image libraries (pydicom, nibabel, pynrrd) and NumPy, which should all be resolved automatically during installation.

Motivation: Voxel Indices, World Coordinates, and Patient Anatomy

When dealing with medical image volumes, one must realize that they live in two different worlds: their voxel space and a world coordinate system, the latter of which has an attached anatomical meaning.

What does a voxel index stand for …

The voxel space very unsurprisingly tells us what image intensity is stored in what sampling position of the volume: Say, we have a three-dimensional array voxel_data_array that contains our image volume:

import numpy as np
voxel_data_array = np.linspace(0, 1, num=1000).reshape(10, 10, 10)
i, j, k = 0, 0, 0
print(voxel_data_array[i, j, k])
# 0.0
i, j, k = 9, 9, 9
print(voxel_data_array[i, j, k])
# 1.0

The call of print(voxel_data_array[i, j, k]) simply tells us we have a value of zero in the [0, 0, 0] corner of the image cube, a value of 1 in the [9, 9, 9] corner, and so on.

… in terms of the physical world?

What the voxel index [i, j, k] does not tell us, is, where in the imaged patient (or healthy subject) we found this value. This, however, may be crucial for medical applications.

Medical image formats therefore provide a mapping from voxel indices to a patient-based world coordinate system: Via rotation, scaling, and translation, we may map from voxel indices to patient coordinates. Using homogeneous coordinates, we can store this mapping in a 4×4 matrix M:

r_11, r_12, r_13, r_21, r_22, r_23, r_31, r_32, r_33 = ...  # rotation
s_i, s_j, s_k = ...  # scaling (world units per voxel)
t_x, t_y, t_z = ...  # translation (world units)

M = [[r_11 * s_i, r_12 * s_j, r_13 * s_k, t_x],
     [r_21 * s_i, r_22 * s_j, r_23 * s_k, t_y],
     [r_31 * s_i, r_32 * s_j, r_33 * s_k, t_z],
     [         0,          0,          0,   1]]
M = np.asarray(M)

We may now use M to find for each voxel index [i, j, k] its respective position [x, y, z] in world coordinates:

homogeneous = lambda c3d: np.r_[c3d, 1]  # append 1 to 3D coordinate
x, y, z = M[:3] @ homogeneous([i, j, k])

Note that, as mentioned, we had to use homogeneous coordinates for the transformation, which explains why we append 1 for the voxel index's coordinate; however, by using only the first three rows of the transformation matrix M, our resulting world coordinate [x, y, z] consists of only three values, as one would expect.

… in terms of the patient?

What so far remains open, is, what does [x, y, z] stand for? We may well have a value in millimeters (or whatever world units are assumed/defined) by now, but we still do not know what the value means in terms of the imaged patient's anatomy. For this reason, medical image formats define the world coordinate system's axes relative to the patient's body axes:

• one world axis points along the patient's left–right axis and therefore its value increases when moving from the patient's left side to their right side – or vice versa,
• one world axis points along the patient's anterior–posterior axis and therefore its value increases when moving from the patient's front to their back – or vice versa,
• one world axis points along the patient's superior–inferior axis and therefore its value increases when moving from the patient's head to their feet – or vice versa.

This may be encapsulated in a definition like "left-posterior-superior (LPS)" or "right-anterior-superior (RAS)". In the "LPS" case, this means "x increases to the left (L), y increases to the back (P as in posterior), z increases to the head (S as in superior)"; in the "RAS" case, this means "x increases to the right (R), y increases to the front (A as in anterior), z increases to the head (S as in superior)". Such a definition of mapping from world coordinate system axes to the patient's anatomy is either implicitly assumed by a particular image format (for example, DICOM uses LPS, NIfTI uses RAS) or explicity stored in the image volume's meta information (NRRD defines the space field in its header for that purpose). Notice we still don't know where we are absolutely positioned within the patient (the origin of the coordinate system is usually defined with respect to some point in the imaging modality, as far as I know), but at least we now have a relative understanding of what a voxel index means with respect to the patient's anatomy.

How can we simplify this in practice?

A remaining open issue is a more practical one: what if we want to display or process an image volume in a certain anatomical orientation? Say, we want to display axial slices of the patient (i.e. slices perpendicular to their superior–inferior axis) or apply a certain image filter along its left–right axis – do we always have to consult the mapping M from voxel indices to world coordinates (or rather, its inverse) in order to find the necessary voxel indices?

Theoretically, indeed we have to do this: notice that a definition like "LPS" does not tell us that the first voxel index increases from the patient's right to their left. For example, if our transformation matrix M looks something like:

M = np.asarray([[0, 0, 1, 0],
                [0, 1, 0, 0],
                [1, 0, 0, 0],
                [0, 0, 0, 1]])

then voxel index i will actually increase with the z coordinate of the world coordinate system (which, remember, in the case of RAS and LPS means moving towards the patient's head):

i, j, k = 1, 0, 0
x, y, z = M[:3] @ homogeneous([i, j, k])
print("x={}, y={}, z={}".format(x, y, z))
# x=0, y=0, z=1

But couldn't we rearrange the voxel array, aligning it with the world coordinate system so that we can be sure increasing voxel index i indeed always means moving to the patient's right side (for an RAS world) or left side (for an LPS world)? We can – precisely, if the rotational part of M contains zeros, ones, and minus ones only; approximately, if it contains arbitrary rotations. And that is where the Volume class comes into play.

The Volume Class

MVloader represents all loaded image volumes as instances of the Volume class.

Volume serves as an abstraction layer to handle transformations from voxel indices to arbitrary anatomical world coordinate systems. In particular

1. the user may choose each Volume's anatomical world coordinate system independent of the underlying file format's world coordinate system, and
2. Volume provides a voxel representation (called aligned_data) with the voxel axes aligned as closely as possible to the user's choice of anatomical world coordinate system. This representation simplifies visualizing volumes in almost correct anatomical orientation (see above) and processing data differently and consistently along different body axes.

Voxel data representations

Each Volume instance provides two representations of the image volume's voxel data as three-dimensional NumPy arrays:

• src_data provides the voxels in the same order as they have been returned by the underlying image library. Specifically, if temporal and/or multiple data dimensions (i.e. vector data) are present, the spatial dimensions will also remain along the same axes in which they have been returned by the underlying image library.
• aligned_data provides the voxels with their axes aligned as closely as possible to the anatomical world coordinate system that has been chosen by the user. If temporal and/or multiple data dimensions (i.e. vector data) are present, the spatial dimensions are brought to the front, the others brought to the back, with the latter keeping their original order. For example, if the user chooses an "LPS" anatomical world coordinate system (meaning that the first coordinate axis will point to the patient's left side, the second axis will point to their back, and the third will point towards their head), then aligned_data[1, 0, 0] will lie to the left of aligned_data[0, 0, 0], aligned_data[0, 1, 0] will lie closer to the patient's back, and aligned_data[0, 0, 1] will lie closer to their head.

Transformation matrix representations

Each Volume instance provides three 4×4 transformation matrices to map from voxel indices to anatomical world coordinates:

• src_transformation maps from src_data's voxel indices to the world coordinate system that has been assumed by the underlying image format (which is provided v