Skip to content
Imaging Evidence Path

How a voxel gets an (x,y,z) position

Tonight · ~30 min · read + a numeric example · energy: medium · setup: none

So the array has no location. How does a voxel acquire one? This lesson builds the affine matrix — the linear map from integer voxel index (i, j, k) to physical patient coordinate (x, y, z) — and works it through real numbers. Once you can write the affine, slice ordering, resampling and CT↔SEG comparison all become arithmetic.

The anchor: where is voxel (0,0,0)?

The patient-coordinate position of a voxel is determined by three things you already met: origin, spacing, and direction. The first voxel (0,0,0) sits at the origin. To reach any other voxel you step along the three axes, each step scaled by that axis’s spacing and pointed along that axis’s direction. String those steps together and you have the affine.

The geometry, in DICOM’s words

This is the durable reference for pixel index ↔ patient coordinate mapping. It follows the DICOM Image Plane Module (PS3.3 C.7.6.2.1.1).

DICOM pixel indices use (i, j) with:

  • i = column index (0 … Columns−1), increasing left→right along a row.
  • j = row index (0 … Rows−1), increasing top→bottom along a column.

ImageOrientationPatient (0020,0037) = six direction cosines:

  • first triplet [rx, ry, rz] = the row direction = direction of increasing i;
  • second triplet [cx, cy, cz] = the column direction = direction of increasing j.

PixelSpacing (0028,0030) = [row spacing, column spacing]:

  • first value = row spacing = distance between adjacent rows = step when row index j increases by 1;
  • second value = column spacing = distance between adjacent columns = step when column index i increases by 1.

The cross-coupling is the trap (we devote the next lesson to it): the row direction cosine pairs with the column spacing, and vice versa.

The in-plane equation

For a single frame, the patient coordinate of the pixel at (i, j) is:

P(i, j) = P0
        + i · (column spacing) · [rx, ry, rz]      # move along a row (i increases)
        + j · (row spacing)    · [cx, cy, cz]      # move along a column (j increases)

where P0 = ImagePositionPatient (0020,0032) of the first (top-left) pixel of this frame. The slice normal adds the third axis (next lesson but one).

A worked numeric example

Take a simple axis-aligned axial slice (the common case), so the row points along +x and the column along +y:

ImageOrientationPatient = [1, 0, 0,  0, 1, 0]      # row=+x, col=+y
PixelSpacing            = [0.8, 0.8]               # [row spacing, col spacing] mm
ImagePositionPatient    = [-200, -150, 42.0]       # P0, LPS mm

Where is the pixel at column i = 100, row j = 200?

P(100, 200) = [-200, -150, 42.0]
            + 100 · 0.8 · [1, 0, 0]     =  [+80,    0,   0]
            + 200 · 0.8 · [0, 1, 0]     =  [   0, +160,  0]
            -------------------------------
            = [-120, +10, 42.0]  mm (LPS)

Notice the z-component is unchanged (42.0) — moving within a slice never changes z. To change z you move to a different slice, which the slice normal handles (lesson 4).

Bridge

Read that result as: “voxel index (100, 200) on this slice is the physical point 120 mm to the right of and 10 mm posterior of the origin, at height 42 mm.” The array index became a millimetre position in the patient. That conversion is the affine doing its one job. Every “is this mask in the right place?” question reduces to applying this map and comparing coordinates.

The affine as a matrix

Pack origin, spacing and direction into one 4×4 matrix A (homogeneous coordinates) and the map becomes a single multiply:

[ x ]     [   D·diag(s)   |  origin ]   [ i ]
[ y ]  =  [               |         ] · [ j ]
[ z ]     [               |         ]   [ k ]
[ 1 ]     [   0   0   0   |    1    ]   [ 1 ]

where D is the 3×3 direction-cosine matrix (columns = axis directions) and diag(s) scales each axis by its spacing. This is exactly what TRACE-CT’s build_affine() constructs from IOP/IPP/spacing, and what SimpleITK carries internally as GetOrigin()/GetSpacing()/GetDirection().

What can go wrong even with the right affine

  • Oblique volumes. If the gantry was tilted, the direction matrix is not identity and the axes are not aligned with +x/+y/+z. Apply the full matrix; do not assume “axial means identity”.
  • Per-slice origin changes. Each slice has its own ImagePositionPatient. The in-plane affine is shared across a regular series, but P0 differs per slice — that is precisely the information lesson 4 uses to order slices.
  • LPS vs RAS. DICOM and SimpleITK are LPS; Slicer is RAS. Convert at the boundary or you mirror the patient.
Stop and think — then reveal

Using the example above, what is the patient coordinate of the origin voxel (i=0, j=0)? And by how many millimetres does the position change when i increases by 1?

(0,0) is just P0 = [-200, -150, 42.0] — the origin is literally the position of the first pixel. Increasing i by 1 adds 1 · (column spacing) · row direction = 1 · 0.8 · [1,0,0] = [0.8, 0, 0], i.e. 0.8 mm along +x. That 0.8 mm is the column spacing. So “spacing” is not an abstract metadata field — it is the physical size of one voxel, and you have just derived it from the geometry.

What to retain

  1. The affine maps voxel index (i,j,k) → patient (x,y,z) via origin + spacing × direction.
  2. In-plane: P(i,j) = P0 + i·colSpacing·rowDir + j·rowSpacing·colDir. The cross-coupling (row direction with column spacing) is the trap of lesson 3.
  3. The origin is the physical position of voxel (0,0,0); spacing is the physical size of one voxel; direction points the axes. Spacing is derivable from the geometry, not just a tag to copy.
  4. Pack the three into a 4×4 affine and the whole map is one matrix multiply — what SimpleITK carries and what build_affine() constructs.

Next: the cross-coupling trap in full — the DICOM row/column trap — and how NumPy and SimpleITK index the same voxel in opposite orders.