Aligning a DICOM to a NIFTI, pixelbypixel
I have struggled at work to find how to align pixels. These pixels are particuarly vexing: they came from a DICOM file or a NIFTI file, which are both formats with awful documentation. If this helps a single person deal with this issue, then this post has reached its goal.
1 The problem
We encountered the following situation at work. We had two files:
 a DICOM file with a MRI scan of the subjects brain.
 a NIFTI file with the segmentation of the brain into regions.
These are essentially 3D arrays of numbers, with a little bit of metadata.
The issue is that the DICOM file and the NIFTI file can have arbitrary orientations with respect to oneanother:
the axes might not be in the same order. For example, maybe we have:
== 400, 300, 200 dicom.shape == 300, 400, 200 nifti.shape
the pixels of each axis might not be in the same order: maybe the first pixel in one format is the last pixel of the other format.
The problem is thus quite simple to state: there are 48 possibilities^{1}. Let’s just open the docs, read the correct header keys, figure out how the arrays are organized, code the correct reordering and call it a day. Right? Right?
Well, that’s certainly a plan. However, No plan survives first contact with the enemy^{2}, and I can state with absolute certainty that the “documentation” of the DICOM and NIFTI standards (and ancillary software) is definitely hostile. Our initial hopes for a quick adventure were thus immediately blown to bits.
Still, I pushed through and, after many bloopers, I have found the solution. I hope it can help somebody else in the future to tackle these f…antastic file formats.
2 The solution
2.1 A bit of geometry
To understand the solution, it is important to understand projective geometry. Do not sweat, we don’t need to understand everything^{3}. We just need to know that, if we want to transform coordinates between two reference frames, then that can be reframed as a matrix multiplication.
A change of reference frame from \(R_1\) to \(R_2\) can be represented by a \(4,4\) (or sometimes \(3,3\)) matrix \(M_{R_2 \leftarrow R_1}\).
Computing the inverse of the matrix gives the matrix for the opposite change of reference frame:
\[ \left(M_{R_2 \leftarrow R_1}\right)^{1} = M_{R_1 \leftarrow R_2} \]
Computing the matrix product \(M_{R_3 \leftarrow R_2} M_{R_2 \leftarrow R_1}\) gives the matrix for change of reference frame from \(R_1\) to \(R_3\).
That’s all we need to know, so feel free to skip to the next section: Section 2.2, or read my detailed explanations below.
For example:
 let \(x,y,z\) denote the coordinates in the MRI room, with respect to the earth: \(x\) is the southnorth axis, \(y\) is the westeast axis, \(z\) is the downup axis. Let the origin be the middle of the door into the room.
 let \(a,b,c\) denote the coordinates in the patient space. \(a\) is the leftright axis of the patient, \(b\) is the backtofront axis, \(c\) is the feettohead axis. The origin point is the middle of the head of the patient.
Then, there exists a matrix \(M\) of shape \(4, 4\) which can be used to translate from \(x,y,z\) coordinates to \(a,b,c\):
\[ \begin{pmatrix} a \\ b \\ c \\ 1 \end{pmatrix} = M \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix} \]
Wait, why do we have constant coordinates \(1\) here? Why are the vectors 4D instead of 3? It’s needed so that we can also represent the change of origin using \(M\). If we are representing transformations between reference frames with the same origin, we can work with a \(3, 3\) matrix instead. However, I wanted to present the 4D case, because it is what is described in the DICOM and NIFTI docs.
Inverting the matrix reverses the direction of the change of variables:
\[\begin{align} \begin{pmatrix} a \\ b \\ c \\ 1 \end{pmatrix} = M \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix} \\ M^{1} \begin{pmatrix} a \\ b \\ c \\ 1 \end{pmatrix} = \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix} \end{align}\]
Similary, if we had three frames of reference, then applying sequentially a change of reference from \(R_1\) to \(R_2\) then a second from \(R_2\) to \(R_3\) would give an overall change from \(R_1\) to \(R_3\). The same property holds for the associated matrices:
\[ M_{R_3 \leftarrow R_2} M_{R_2 \leftarrow R_1} = M_{R_3 \leftarrow R_1} \]
2.2 Finding the matrices
So now we know that we need to go looking for matrices.
For the NIFTI file format, this is immediate: the matrix is encoded as a header key^{4}. If using nibabel, it is immediately accessible.
Carefully reading the docs specifies that this affine matrix specifies the transformation between the voxel indices \(i,j,k\) and the RAS patientspace (the acronym gives the order and direction of the axes: Right then Anterior (i.e. backtofront) then Superior (foottohead)).
Surely, the dicom format must be similarly simple. Nope. However, since I know have digested the docs, here are the steps:
First, there is a transposition between what the DICOM format calls voxels and how the array is organized on disk^{5}. This is explained here.
Then, the
ImageOrientationPatient
header specifies the first two columns of the matrix.= np.array(dicom.ImageOrientationPatient[:3]) a = np.array(dicom.ImageOrientationPatient[3:]) b 1= np.zeros((4, 4)) matrix 3, 3] = 1 matrix[ 23, 0] = b matrix[:3, 0] = a matrix[:
 1
 Initializing the matrix and specifying the fourth column.
 2

Note the change of order with respect to
ImageOrientationPatient
. This is due to the transposition.
Finally, by considering the change of position between two different dicom slices, we can find the third column.
= (np.array(dicom2.ImagePositionPatient)  np.array(dicom.ImagePositionPatient)) / ( slice_diff  dicom.InstanceNumber dicom2.InstanceNumber )1= slice_diff / np.sum(slice_diff**2) ** 0.5 c
 1
 Normalizing the vector so that it has norm 1.
This gives us the matrix to transform from the DICOM array coordinates \(i,j,k\) to the LPS patientspace. This is almost the same as the RAS space used by the NIFTI format: the first two axes are just pointing in the opposite direction.
Overall, we are now able to combine:
 The matrix from DICOM coordinates to LPS,
 The matrix from LPS to RAS,
 The matrix from NIFTI coordinates to RAS,
in order to find the matrix corresponding to the change of variable we want:
\[ M_{\text{DICOM} \leftarrow \text{NIFTI}} = \left[ M_{\text{RAS} \leftarrow \text{LPS}} M_{\text{LPS} \leftarrow \text{DICOM}} \right]^{1} M_{\text{RAS} \leftarrow \text{NIFTI}} \]
And that’s it. We can now analyze the matrix to find how it swaps axes around and reorders them, and apply that transformation to the NIFTI array:
= ...
matrix_dicom_from_nifti = nib.orientations.io_orientation(matrix_dicom_from_nifti)
ornt = nib.orientations.apply_orientation(nifti, ornt) reoriented_nifti
Honestly, the only way you’ve made it this far is if you are yourself trying to deal with this exact problem. If so, then best of luck and bon courage. You will need both.
3 Bloopers
I can’t resist but tell you about all of the hilarious moments along the way were the DICOM spec blew up in my face.
The DICOM format does not define an ordering of DICOM slices. This means that different tools could choose different orderings. If you struggle with a mismatch of direction along the axis over which the slices are gathered, then the underlying issue might be that your tools do not align along this axis in the same way. NB: the difference in ordering in our case manifests in less than 5% of cases, so that was very tricky to debug.
After a little bit of (unsuccesfully) poking around and trying to read the docs, I asked a colleague about the issue. We quickly went from: “Oh, but all dicoms have the same orientation.” to “No wait, there are two.” to “Well technically it’s maybe three.”. We then found a fourth one later in our tests.
The nibabel documentation, and the NIFTI docs I have found all repeat: “the NIFTI format uses a RAS reference frame” throughout. Imagine my surprise when I discovered the existence of the
affine
header.In my first calculation of the shift between two DICOM slices, I initially naively thought that the slices would be in order and that the first file (with name
slice_00000
) would actually be the first slice. Ha. Ha. Ha. They are in a random order instead.
4 References
Footnotes
\(3!=6\) possibilities for the order of the axes; \(2^3=8\) possibilities for the order of each axis.↩︎
Seriously, don’t worry if you don’t quite get it: this projective geometry stuff is a bit crazy.↩︎
Technically three but nibabel automatically chooses the most appropriate one.↩︎
Possibly due to the differences between Fortran and C array layouts on disk?↩︎