Geometry, segmentation and resampling — reference
Reference · Dense lookup material — return here later. The teaching versions live in the lessons.
Coordinate systems
| System | What it is | Order / unit |
|---|---|---|
| Array indices | voxel position in the NumPy array | [k, j, i] = slice, row, column (integers) |
| Physical patient | position in the patient | (x, y, z) LPS mm (continuous) |
| Display | screen pixels | view-dependent |
LPS: +x → Left, +y → Posterior, +z → Superior. 3D Slicer uses RAS (sign flip on x, y). NumPy has no coordinate system at all.
The index ↔ patient equation (Image Plane Module, PS3.3 C.7.6.2.1.1)
DICOM pixel indices: i = column, j = row.
ImageOrientationPatient(0020,0037)=[rx,ry,rz, cx,cy,cz]: first triplet = row direction (increasingi); second = column direction (increasingj).PixelSpacing(0028,0030)=[row spacing, col spacing].
In-plane:
P(i, j) = P0 + i·(colSpacing)·[rx,ry,rz] + j·(rowSpacing)·[cx,cy,cz]P0 = ImagePositionPatient (0020,0032) of the top-left pixel. Slice normal
n = cross([rx,ry,rz],[cx,cy,cz]); slice step = slice spacing; slice index k →
P0_k = ImagePositionPatient of slice k.
Index orders (the trap)
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)img.GetSize() is (x,y,z); arr.shape is (z,y,x). Always print
GetSpacing()/GetOrigin()/GetDirection(); convert at the boundary with
GetArrayFromImage / GetImageFromArray.
Slice ordering
t_k = dot(ImagePositionPatient_k, slice_normal); sort slices by t_k ascending.
Spacing truth = t_{k+1} − t_k; SpacingBetweenSlices is only a claim.
InstanceNumber and filenames are not geometry.
Frame of Reference decision
same FoR ≠ same voxel grid (verify origin/spacing/direction)different FoR ≠ impossible to relate (explicit registration can map them)different FoR ⇒ do NOT overlay without an explicit spatial relationshipPairing classification: EXACT_GRID / RESAMPLING_REQUIRED / INCOMPATIBLE,
decided by orientations_equivalent() / spacings_equivalent() /
positions_equivalent().
Segmentation representations
| Representation | Stores | Carries | Loses on conversion |
|---|---|---|---|
| Binary mask | voxel labels | image geometry (if kept with image) | source references; fractional edges |
| DICOM SEG | per-frame labelled planes | source SOP refs; per-frame geometry; segment identity | fractional→binary; off-grid snapping |
| RTSTRUCT | coplanar contours | contour points + reference grid | rasterisation is lossy/tool-dependent |
SEG Segmentation Type (0062,0001): BINARY / FRACTIONAL / LABELMAP.
Resampling
- Image: linear or B-spline (radiomics prefers B-spline, then clip).
- Mask: nearest-neighbour only. Never linear-interpolate a mask.
- B-spline overshoots at sharp boundaries → clip HU afterward (TRACE-CT:
[-1000, +400]). - Report the target spacing/interpolator with every result.
resampled_img = sitk.Resample(img, ref, sitk.Transform(), sitk.sitkBSpline, 0.0, sitk.sitkInt16)resampled_mask = sitk.Resample(mask, ref, sitk.Transform(), sitk.sitkNearestNeighbor, 0, sitk.sitkUInt8)What usually goes wrong
- Index vs physical mix-up — using array indices where physical mm are needed.
- LPS vs RAS confusion (DICOM/SimpleITK LPS; Slicer RAS; NumPy neither).
- Ordering by
InstanceNumberinstead of slice-normal projection. - Trusting FoR equality without verifying origin/spacing/direction.
- Linear-interpolating a mask — destroys label integrity.
- Comparing array shapes instead of physical geometry.
- B-spline overshoot at sharp boundaries left unclipped.
- Resampling to a target spacing but not reporting it.
TRACE-CT connection (read-only references)
- Slice normal —
seg_decode_resample.pycompute_slice_normal(). - Affine —
seg_decode_resample.pybuild_affine();rider_radiomics_extraction.pyaffine_to_sitk_geometry(). - Geometric ordering —
project_position()(IPP projected onto the slice normal). - Geometry equivalence —
orientations_equivalent(),spacings_equivalent(),positions_equivalent();EXACT_GRID/RESAMPLEDinclassify_transfer_pair(). - FoR — carried on every
CTVolume; checked in preflight, never trusted alone. - SEG decode + mask transfer —
decode_seg_file()(per-frame functional groups, source SOP UIDs),build_source_mask_on_ct_grid(). - Inverse NN mask resample —
resample_mask_inverse_nn()(label-preserving, bitwise-verified byverify_determinism()). - Image resample + clip —
resample_image_isotropic()(B-spline to 1 mm) +clip_intensities()inrider_radiomics_extraction.py. - Visual QC —
scripts/visual_qc/*andtools/slicer/open_qc_bundle.py.
Go deeper
- SimpleITK fundamentals docs (origin/spacing/direction, physical vs index space, resampling) — the single best reference for the affine and resampling.
- DICOM PS3.3 Image Plane Module + Multi-frame Functional Group macros — normative IOP/IPP and SEG per-frame geometry.
- dcmqi / highdicom docs — SEG encoding/decoding beyond pydicom.
- IBSI image-processing workflow (Zwanenburg et al. 2020) — which resampling/discretisation choices a radiomics pipeline must pin (Chapter 3).
What to retain (chapter summary)
- Three coordinate systems: array indices, physical patient (LPS), display. The affine connects the first two.
- Geometry = origin + spacing + direction. Slice normal =
cross(row, col); order slices by projection onto it — never byInstanceNumber. FrameOfReferenceUIDis a fast rejection / slow acceptance test; always verify geometry.- Decode SEGs by per-frame functional groups and source SOP UIDs; round-trips lose fractional edges and off-grid frames.
- Images interpolate (linear/B-spline); masks resample nearest-neighbour. Same dimensions ≠ same physical space.
- Resampling changes features. Pin and report every preprocessing choice.