The DICOM row/column trap
Tonight · ~25 min · read · energy: medium · setup: none
You wrote arr[x, y] to index column x, row y of a CT — and the anatomy came
out transposed. Or you read PixelSpacing, assumed the first number was the x
spacing, and your resampling stretched the image the wrong way. Both are the same
bug: DICOM’s row/column naming crosses itself, and it crosses again against
NumPy. This lesson is the one place to nail down the conventions so you never
transpose by accident again.
The anchor: the same voxel, three naming systems
One physical voxel is described by three different tools using three different index orders. They all mean the same voxel, but they index it differently:
DICOM pixel (i, j) = (column i, row j)SimpleITK index (x, y, z) = (column i, row j, slice k) # NOT NumPy orderNumPy array [k, j, i] = (slice k, row j, column i)NumPy and SimpleITK index the same voxel in opposite orders. That is the bite.
The cross-coupling, spelled out
The confusion has two layers. Layer 1 — the names cross:
ImageOrientationPatientfirst triplet = row direction = direction of increasing column indexi.PixelSpacingfirst value = row spacing = step when row indexjincreases by 1.
So the row direction cosine pairs with the column spacing (it is the step along
i), and the column direction cosine pairs with the row spacing (the step
along j). The word “row” appears in both but refers to different axes. This is
not pedantry — getting it backwards resamples or registers along the wrong axis.
Layer 2 — the index orders cross between libraries:
A DICOM pixel at (i = column, j = row) lives in a NumPy array (slices-first,
image order) as numpy[k, j, i], i.e. numpy[slice, row, column]. So arr[0] is
slice 0, and within a slice arr[k, j, i] indexes row j, column i — the
transpose of the naive arr[x, y].
flowchart TD
D["DICOM: pixel (i=col, j=row)"]
N["NumPy: arr[k, j, i]<br/>= slice, row, column"]
S["SimpleITK: GetPixel(x=i, y=j, z=k)<br/>= column, row, slice"]
D -- "same voxel" --> N
D -- "same voxel" --> S
N -. "opposite index order" .-> S
The one-line rule
When in doubt, hold this in your head:
DICOM
iis column,jis row;PixelSpacingis[row, col]; NumPy is[k, j, i]; SimpleITK is(x=i, y=j, z=k).
And the operational version: always print GetSpacing()/GetOrigin()/GetDirection()
rather than assuming, and convert explicitly with sitk.GetArrayFromImage /
GetImageFromArray at the NumPy↔SimpleITK boundary.
A concrete check you can run
import SimpleITK as sitkimg = sitk.ReadImage("series_dir") # or build from DICOMarr = sitk.GetArrayFromImage(img) # -> NumPy [k, j, i]
print(img.GetSize()) # (x=cols, y=rows, z=slices) -- SimpleITK orderprint(arr.shape) # (slices, rows, cols) -- NumPy order (reversed)print(img.GetSpacing()) # (sx, sy, sz) in mm, LPSNote the reversal: img.GetSize() is (x, y, z) while arr.shape is (z, y, x).
The voxel at SimpleITK (x=10, y=20, z=5) is arr[5, 20, 10]. If you forget the
reversal you have silently transposed the volume.
Recognising orientations from the cosines
The two IOP triplets tell you which way the row/column axes point in the patient. For a standard axial series the row points left→right and the column points anterior→posterior, giving the axial IOP we used in lesson 2:
ImageOrientationPatient = [1, 0, 0, 0, 1, 0] # row=+x(L), col=+y(P), axial planeA coronal or sagittal series has different cosines; an oblique series (acquired at a tilt) has non-axis cosines. Never assume “axial” — read the cosines, because the direction matrix is part of the geometry that decides whether two objects share space (lesson 5).
Stop and think — then reveal
You read PixelSpacing = [0.8, 0.7] and need the step in millimetres when the
column index i increases by 1. Which value do you use — 0.8 or 0.7?
0.7 — the second value (column spacing). The column index i increases along
a row, so the step is the column spacing, which is the second element of
PixelSpacing. The first value (0.8, row spacing) is the step when the row index
j increases. Mixing these up stretches or compresses your resample along the
wrong axis — a silent error that looks fine until you overlay a mask.
What to retain
- The same voxel has three index orders: DICOM
(i=col, j=row), SimpleITK(x=i, y=j, z=k), NumPy[k, j, i]. NumPy and SimpleITK are reversed. PixelSpacingis[row spacing, col spacing]; the row direction cosine pairs with the column spacing. The naming crosses itself.img.GetSize()is(x,y,z);arr.shapeis(z,y,x). Convert at the boundary and print the geometry rather than assuming.- Never assume “axial” — read
ImageOrientationPatient; oblique volumes have non-identity direction matrices.
Next: with index/physical space sorted, the question that makes or breaks a 3D volume — how to order slices without trusting filenames.