The NumPy array has no physical location

Tonight · ~20 min · read · energy: low · setup: none

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. This lesson is about the single fact that prevents most geometry bugs: a NumPy array is a box of numbers with no position in the patient.

The anchor: 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.

Three coordinate systems

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

flowchart LR
    A["(1) Array / voxel indices<br/>(col i, row j, slice k)<br/>integers, shape W,H,D"]
    B["(2) Physical / patient coords<br/>(x, y, z) in mm, LPS<br/>continuous, real numbers"]
    C["(3) Display coords<br/>screen pixels<br/>view-dependent"]
    A -- affine matrix --> B
    B -- window/level, zoom --> C

We care about the first two. The transform between them is the affine matrix (built explicitly in lesson 2).

Array indices

pixel_array[k, j, i] (image order: slice, row, column) is indexed by integers. The shape is (slices, rows, columns). Indices tell you nothing about physical size. Changing the array’s shape (cropping, resampling) changes its indices without changing any real-world position unless you also update the geometry.

Physical patient coordinates (LPS)

DICOM defines the patient coordinate system as LPS:

• +x → Left
• +y → Posterior
• +z → Superior

Every ImagePositionPatient and ImageOrientationPatient value is in **LPS millimetres`. This is where geometry becomes comparable across objects.

The two traps waiting for you

• 3D Slicer works internally in RAS (Right-Anterior-Superior), converting LPS↔RAS at its boundary. The sign flip on x and y bites people regularly.
• 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/direction. SimpleITK stays in LPS, matching DICOM, which is why we use it as the bridge between arrays and physical space.

Geometry is three things

A volume’s geometry is fully described by three quantities, each defined in patient space:

• Origin — the physical coordinate of voxel index (0,0,0).
• Spacing — the physical size of one voxel: (sx, sy, sz) in mm along the three axes.
• Direction — a 3×3 matrix whose columns are the direction cosines of the volume’s three axes in patient space. Identity if the volume is axis-aligned with the scanner; non-identity if the volume is oblique (tilted gantry).

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

“Same dimensions” implies nothing

Hold on to this fact; we will use it twice:

Two images are spatially comparable only when their origin + spacing + direction agree (within tolerance). Identical array shape implies nothing about spatial equivalence.

We use it once to decide whether a segmentation already shares the CT grid (no resample needed — lesson 5), 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.

Stop and think — then reveal

You compute diff = array_A - array_B on two (100, 512, 512) CTs and get plausible-looking numbers. What two questions must you answer before trusting diff as a real voxel-wise comparison?

1. Are they the same patient? (A difference image of two different patients is meaningless.)
2. Do they share geometry — same origin, spacing, and direction? If the spacings differ, voxel [k,j,i] in A is a different physical location than voxel [k,j,i] in B, so the subtraction aligns the wrong anatomy. Even the same patient, same series, re-exported at a different spacing, makes this subtraction a geometric warp rather than a comparison. The array shape told you none of this.

What to retain

1. A NumPy array has no physical location — shape tells you nothing about space.
2. Three coordinate systems matter: array indices, physical patient (LPS), display. The affine connects the first two.
3. Geometry = origin + spacing + direction. That trio is the entire vocabulary of medical-image space.
4. Two objects are spatially comparable only when their geometry agrees — never trust array shape as a proxy for alignment.

Next: build the thing that connects indices to patient coordinates — the affine, with a worked numeric example.