The NumPy array has no physical location

A segmentation lands a few centimetres off the anatomy it was drawn for. The fix is rarely where people look first: the mask’s values are fine, the drawing was fine, and the code ran without an error. The failure is upstream of all of that — the mask was treated as if it lived in the patient, when all it ever carried was a position in an array. This section is about the single fact that prevents many errors in this class: a NumPy array is a box of numbers with no position in the patient.

The same fact decides quieter questions. You load a CT into NumPy and get an array of shape (100, 512, 512). A colleague loads a different scan and also gets (100, 512, 512). Are the two arrays comparable? Can you subtract one from the other? You cannot know from the arrays.

Same shape, different universe

Two arrays of identical shape can be:

• the same patient, same series, different reconstructions (different spacing);
• the same patient scanned twice (different positions in the scanner);
• different patients entirely;
• a 30 cm body CT and a 5 cm specimen.

The array carries none of that. It carries rank, shape, dtype, and values — nothing about millimetres, nothing about orientation, nothing about which patient or which way is up. To reason spatially you must attach a geometry to it.

The same index is a different physical point when origin, spacing or direction changes; shape alone supplies no patient-space alignment.

Three coordinate systems

A medical image implicitly lives in three coordinate systems, and confusing them is the root of most geometry bugs.

Array / voxel indices
(col i, row j, slice k)
integers, shape W,H,D
        |
        | affine matrix
        v
Physical / patient coords
(x, y, z) in mm, LPS
continuous, real numbers
        |
        | window/level, zoom
        v
Display coords
screen pixels
view-dependent

The first system is the array itself: pixel_array[k, j, i], the ordering read as (slice, row, column) in the index-to-patient convention this chapter uses, indexed by integers with shape (slices, rows, columns). Indices tell you nothing about physical size. Cropping or resampling the array changes its indices, and the index-to-physical map changes with them unless origin, spacing, direction, and size are correctly propagated or recomputed so that every intended physical location stays paired with the voxel that represents it. The second system is physical patient space, measured in continuous millimetres. DICOM defines the patient coordinate system as LPS: +x → Left, +y → Posterior, +z → Head — a shorthand for ordinary biped (human) data. ImagePositionPatient values are coordinates in LPS millimetres, while ImageOrientationPatient values are dimensionless direction cosines stating where the image’s row and column axes point in that frame, not millimetres. This is where geometry becomes comparable across objects. The third system — display coordinates, screen pixels, view-dependent — is an interaction interface: physical measurements are defined in patient/image physical space, and screen coordinates are typically how a user selects the patient-space points to measure. This chapter leaves that interface behind.

Two boundary traps follow directly. 3D Slicer works internally in RAS (Right-Anterior-Superior), converting LPS↔RAS at its boundary; omitting that conversion, or applying it twice, reflects or misplaces coordinates along x and y. NumPy has no coordinate system at all — only array axes. That is precisely why a bare ndarray is geometrically meaningless until you attach origin, spacing, and direction. SimpleITK stays in LPS, matching DICOM, which is why it serves as the bridge between arrays and physical space.

The transform connecting the first system to the second is the affine matrix, built explicitly in the next section.

Geometry is four things

A complete regular image grid is described by four quantities. Origin, spacing, and direction are defined in patient space; size is the extent of the index space they map from:

• Origin — the physical coordinate of voxel index (0,0,0).
• Spacing — the centre-to-centre step between adjacent samples: (sx, sy, sz) in mm along the three axes. It is the distance between neighbouring sample points, not universally the physical support — the tissue volume — a voxel covers.
• Direction — a 3×3 matrix whose columns are the direction cosines of the volume’s three axes in patient space. The identity matrix means the volume’s axes align with the selected physical axes, not necessarily the scanner itself; non-identity may include obliquity (tilted gantry), a permutation of the axes, or a reflection.
• Size — the number of samples along each axis: (slices, rows, columns), bounding which indices exist.

These four are the entire vocabulary of medical-image geometry. Origin answers where, spacing answers how far apart, direction answers which way, and size answers how many. Every operation in this chapter — ordering slices, comparing a CT and a SEG, resampling — is a manipulation of these four.

Same dimensions imply nothing

Hold on to this fact; this chapter uses it twice:

Two images occupy the same grid only when their size and their origin + spacing + direction all agree within prespecified tolerances. That is exact-grid equality. Spatial comparability is the weaker property: two different grids in the same patient frame can still be compared, through their coordinates or by resampling. Identical array shape implies nothing about either.

It is used once to decide whether a segmentation already shares the CT grid, so that no resampling is needed (the shared-space section), and once to decide whether two volumes are co-registered. The reflex to build: when you are about to do anything voxel-wise between two objects, first compare their geometry, not their shape.

A subtraction makes the stakes concrete. Suppose you compute diff = array_A - array_B on two (100, 512, 512) CTs and get plausible-looking numbers. Before trusting diff as a real voxel-wise comparison, two questions need answers. Are they the same patient? A difference image of two different patients is meaningless. And do their grids match — the same size, origin, spacing, and direction, within prespecified tolerances? If the grids differ, voxel [k, j, i] in A is a different physical location than voxel [k, j, i] in B, so the same-index subtraction pairs voxels that occupy different physical points: an unregistered, index-mismatched comparison, and not itself a geometric warp, since nothing in the arrays was moved or resampled. Same patient, same series, re-exported at a different spacing is exactly this case — the two remain comparable, through their coordinates or by resampling, but not by pairing identical indices. The array shape supplied none of the information needed to tell.

So the array still has no location. What it has, once the DICOM metadata is read, is a recipe for acquiring one. The next section builds that recipe — the affine, with a worked numeric example.

Sources

• DICOM PS3.3, Image Plane Module (C.7.6.2) and Image Position and Image Orientation (C.7.6.2.1.1) — ImagePositionPatient, ImageOrientationPatient, and the pixel-index ↔ patient-coordinate mapping. (current edition)
• SimpleITK, Fundamental Concepts — origin, spacing, size, direction, and the index-to-physical-point mapping. (docs)
• 3D Slicer, Coordinate systems — the LPS and RAS conventions and their relationship. (docs)