DICOM coordinates and gamma, illustrated#
A dose value and its position belong together. A file may store voxels from left to right, from right to left, or with rows and columns exchanged. Reordering those values for calculation must preserve which dose belongs at which physical position. It does not move dose within the patient.
This notebook explains that distinction for a physicist who need not read Python. Read the text and figures in order; expand the hidden code to inspect a calculation. All dose fields are synthetic and generated locally. This is the main explanation of the geometry, historical impact and plotting. The separate validation record contains test evidence, fixture provenance, reproduction commands and remaining limitations.
The main examples cover the DICOM definition (section 1), independent checks (section 2), historical errors and edge cases (section 3), gamma and plotting (sections 4–5), and checks on earlier results (section 6). The appendix explains the former interpolation failures; performance is covered separately. Historical claims below concern PyMedPhys 0.41.0; they do not establish the behaviour of every earlier release.
Start with a coordinate–dose pair#
The example below shows three voxels in two array orders. Follow the 10 Gy voxel: its index changes from 0 to 2, but it remains at x = 100 mm. The arrows follow coordinate–dose pairs, not physical motion.
Blue denotes the current conversion, orange the previous conversion, and dashed black lines an independent DICOM calculation. Examples use two Numba threads to suit documentation build hosts.
as supplied: [100. 99. 98.] [10. 20. 30.]
prepared: [ 98. 99. 100.] [30. 20. 10.] (interpolator: pymedphys)
Patient directions and file terminology#
For the human (BIPED) DICOM patient coordinate system, +x is towards the patient’s left, +y towards the posterior, and +z towards the head. These anatomical directions do not change when the patient is prone or enters the scanner feet first. A common coordinate system is distinct from anatomical registration: this conversion does not align scans or account for patient movement.
Abbreviation |
Patient position |
|---|---|
HFS / HFP |
Head first supine / head first prone |
FFS / FFP |
Feet first supine / feet first prone |
HFDL / HFDR |
Head first decubitus left / right (lying on the left / right side) |
FFDL / FFDR |
Feet first decubitus left / right |
The row direction runs along a row, towards increasing column index. The column direction runs towards increasing row index. In contrast, the first value in PixelSpacing is the distance between rows. The equation below makes these pairings explicit.
We write raw data as pixel_array[k, row, column], where k selects a stored slice. After conversion we write dose[i_z, i_y, i_x], whose coordinates are (z[i_z], y[i_y], x[i_x]). These indices describe different array layouts.
The hidden setup builds DICOM datasets and evaluates the standard independently. Distinct stored values identify each voxel, so the checks follow its dose through the reordering.
The historical comparison uses the 0.41.0 coordinate function bodies, with a legacy_ prefix and docstrings omitted. The dose wrapper is simplified: it scales pixel_array directly and omits the original transfer-syntax helper, because every synthetic dataset already declares its transfer syntax. This reproduces coordinate behaviour, not the whole historical dose-loading path.
1. From stored pixels to patient coordinates#
DICOM defines the patient position of the stored voxel in slice \(k\), row \(r_i\), and column \(c_i\) of pixel_array[k, row, column] as
where:
\(\mathbf{S}\) is
ImagePositionPatient, the centre of the first stored voxel;\(\hat{\mathbf{r}}\) and \(\hat{\mathbf{c}}\) are the first and second triplets of
ImageOrientationPatient, the directions of increasing column index \(c_i\) and increasing row index \(r_i\);\(\Delta_r\) is
PixelSpacing[0], the spacing between rows, and \(\Delta_c\) isPixelSpacing[1], the spacing between columns;\(g_k\) is
GridFrameOffsetVector[k], the offset of slice \(k\) along \(\hat{\mathbf{n}}\). A vector whose first element is not zero holds absolute z coordinates instead, which the standard permits only forImageOrientationPatient\(= (1, 0, 0, 0, 1, 0)\), with the first element equal to \(S_z\); subtracting \(S_z\) converts it to relative offsets. A positive offset follows the normal; decreasing offsets reverse traversal.
Equivalently, as a homogeneous matrix,
Changing the orientation changes the direction cosines, never the sign of \(\mathbf{S}\). The figure below draws a stored grid of four rows and five columns in each of the eight supported orientations, in the transverse plane as it is usually displayed: patient left to the right, posterior downwards.
Sources: DICOM voxel position, pixel spacing, and RT Dose slice offsets.
For every supported orientation, \(\hat{\mathbf{r}}\), \(\hat{\mathbf{c}}\), and \(\hat{\mathbf{n}}\) are signed unit vectors along the patient axes, so the matrix is a signed permutation with scaling. Each patient coordinate therefore depends on exactly one stored index. The table below is computed from ImageOrientationPatient alone; with increasing slice offsets it matches the table in the validation note.
orientation columns run rows run positive offset
HFS +x +y +z
HFP -x -y +z
FFS -x +y -z
FFP +x -y -z
HFDL -y +x +z
HFDR +y -x +z
FFDL +y +x -z
FFDR -y -x -z
1.1 How PyMedPhys applies the definition#
Because the matrix is a signed permutation, PyMedPhys evaluates one axis per patient coordinate instead of a coordinate for every voxel: \(x\), \(y\), and \(z\) each come from the single stored index they depend on, in \(O(\text{rows} + \text{columns} + \text{slices})\) operations. pymedphys.dicom.zyx_and_dose_from_dataset then
transposes the stored (slice, row, column) dimensions into patient \((z, y, x)\) order, which swaps rows and columns for the decubitus orientations; and
reverses each axis that decreases, together with the matching dimension of the dose.
The result satisfies dose[i_z, i_y, i_x] \(\leftrightarrow\) (z[i_z], y[i_y], x[i_x]) with every axis ascending. Nothing is interpolated: each stored value is moved, unchanged, to the array position that matches its coordinates. Below, each stored voxel of a small grid is labelled with its storage index, so the rearrangement can be read directly.
Voxel 0 is the first stored voxel, at \(\mathbf{S} = (40, -30, 12)\) mm. For head first prone the columns run towards \(-x\) and the rows towards \(-y\), so both in-plane axes are reversed and voxel 0 moves to the corner with the largest \(x\) and \(y\). For head first decubitus left the rows run along \(+x\) and the columns along \(-y\), so the dimensions are also swapped and the returned array has shape (3, 5, 4) rather than (3, 4, 5).
2. The current conversion places every voxel correctly#
2.1 Every voxel against the DICOM definition#
The first check converts a grid with an off-centre origin, unequal row and column spacing, unequal numbers of rows and columns, and unevenly spaced slices, stored in each orientation with increasing and with decreasing slice offsets. It reports the largest distance between any returned voxel and the position that voxel_positions calculates for it.
orientation slice offsets current previous
HFS increasing 0 mm 0.0 mm
HFS decreasing 0 mm 0.0 mm
HFP increasing 0 mm 446.3 mm
HFP decreasing 0 mm 446.3 mm
FFS increasing 0 mm 256.5 mm
FFS decreasing 0 mm 261.2 mm
FFP increasing 0 mm 394.2 mm
FFP decreasing 0 mm 397.3 mm
HFDL increasing 0 mm axes did not fit the dose
HFDL decreasing 0 mm axes did not fit the dose
HFDR increasing 0 mm axes did not fit the dose
HFDR decreasing 0 mm axes did not fit the dose
FFDL increasing 0 mm axes did not fit the dose
FFDL decreasing 0 mm axes did not fit the dose
FFDR increasing 0 mm axes did not fit the dose
FFDR decreasing 0 mm axes did not fit the dose
The current conversion places every voxel exactly: the differences are zero, not merely small, because both calculations add the same signed multiples of the spacing to the same origin. The previous conversion was exact only for head first supine. For the other non-decubitus orientations it displaced the grid by hundreds of millimetres, and for the decubitus orientations its axes did not even have the lengths of the dose dimensions they were paired with.
2.2 Randomised grids#
The same comparison over 800 random grids, 100 per orientation, with 2 to 8 voxels along each dimension, spacings from 0.2 mm to 5 mm, origins within ±1500 mm, and uneven slice offsets that increase or decrease. Some head first supine grids use absolute slice offsets.
Largest current placement error in 800 grids: 0 mm
The vertical scale is linear from 0 to 1 mm and logarithmic above 1 mm, so zero and large errors can share a panel. These deliberately broad synthetic origins stress the calculation; this is not a survey of clinical errors.
Head first supine grids with relative offsets were placed exactly before as well; the displaced head first supine grids are those with absolute slice offsets (section 3.4). Decubitus grids with equal numbers of rows and columns were accepted by the previous conversion but transposed in the transverse plane, so they appear as finite displacements.
2.3 One dose distribution, eight storage orders#
A stronger test starts from a known physical dose distribution and stores it in each orientation. The storage directions below are tabulated by hand from section 1, and the encoding uses neither voxel_positions nor PyMedPhys. The distribution is shaped like the letter F, which has no mirror or rotational symmetry, so any flip or transposition would be visible. It also rises asymmetrically along z, so a slice reversal changes the pattern. Its edges come from a continuous Gaussian-smoothed field evaluated at the sample coordinates; a requested 3 mm shift remains 3 mm even on a 2.5 mm grid.
The hand-tabulated storage directions agree with ImageOrientationPatient.
16 of 16 encodings recover exactly the same axes and dose
FFS: stored slice k=0 is at z=27 mm; k=4 is at z=15 mm. Both vertical axes increase upwards here.
The stored images differ by rotations, reflections, and transpositions; the returned arrays are identical to each other and to the original, bit for bit, in every orientation and for both slice orders.
2.4 Where the previous conversion placed the same dose#
The figure projects the central slice’s isodose at 50% of the whole-volume maximum onto the xy plane. The previous and correct placements can have different z coordinates: both appear in the panel titles. This compares in-plane placement; it is not an overlay on one common anatomical slice.
The current conversion (blue) coincides with the DICOM position (dashed) in every orientation. The previous conversion (orange):
placed head first supine correctly;
for the uniformly spaced head first prone, feet first supine, and feet first prone grids shown here, reflected the grid centre through the origin along each reversed axis, so the dose was translated, not mirrored, by twice the centre coordinate: \(x' = x - 2c_x\) on a reversed x axis, with \(c_x\) the centre of the grid;
for the decubitus orientations, also transposed the transverse plane, because it paired the stored rows with \(y\) and the columns with \(x\).
Why the translation? For a uniformly spaced reversed axis the true coordinates in storage order are \(S - j\Delta\). The previous code formed the image-aligned axis \(-S + j\Delta\), appropriate for its IEC fixed output, and then reversed the array for DICOM output. Reversing restores the direction of travel but not the sign of the origin, so storage index \(j\) received \(-S + (n - 1 - j)\Delta\), which differs from the truth by \((n - 1)\Delta - 2S = -2c\).
3. Edge cases#
3.1 Cropping and shifting a grid#
Because the error is \(-2c\), it depends on where the grid is and how far it extends. The smallest example is a single head first prone row of five 1 mm columns starting at \(x = 100\) mm, with and without its last two columns.
five columns:
DICOM x of stored columns 0, 1, ... [100. 99. 98. 97. 96.]
previous x of stored columns 0, 1, ... [ -96. -97. -98. -99. -100.] error -196 mm
current ascending x [ 96. 97. 98. 99. 100.], holding stored columns [4 3 2 1 0]
three columns:
DICOM x of stored columns 0, 1, ... [100. 99. 98.]
previous x of stored columns 0, 1, ... [ -98. -99. -100.] error -198 mm
current ascending x [ 98. 99. 100.], holding stored columns [2 1 0]
The previous error is \(-196\) mm for the full row and \(-198\) mm for the cropped one, so the two grids, both head first prone, disagree by 2 mm. A common constant translation preserves every reference/evaluation distance, so that particular coordinate error cancels in gamma. Merely sharing metadata does not establish this: uneven reversed slices can distort distances, and decubitus errors can alter in-plane geometry. A shared error also does not repair the former interpolator. Crops, different extents and shifts need their own checks.
The next figure repeats this on a three-dimensional head first prone grid with a smoothed 50 mm × 40 mm field: once with 10 mm removed from its \(+x\) edge, and once with the whole grid, and its dose, moved 15 mm towards \(+x\).
The current conversion reports the crop as no movement and the shift as \(+15\) mm, as it should. The previous conversion moved the cropped dose by 10 mm and moved the shifted dose by \(-15\) mm, in the wrong direction; both previous profiles also sit far from their true positions (dotted).
Gamma inherits these errors when its coordinates come from the conversion. Below, the cropped grid is the reference and the full grid the evaluation, at 3%/3 mm with a 10% lower dose cutoff. The previous coordinates are passed to the current pymedphys.gamma, which accepts their descending axes; with its default interpolator, version 0.41.0 raised an error for this evaluation grid instead (Appendix A), and with interp_algo="scipy" it returned the same values as below.
For a reference position \(\mathbf{r}\), gamma balances spatial separation and dose difference:
Here \(\mathbf{e}\) ranges over the interpolated evaluation field, \(\Delta d=3\) mm, and global \(\Delta D=0.03\max(D_{\mathrm{ref}})\). A value \(\gamma\leq1\) passes these criteria. Examples use global normalisation unless labelled local. A 10% cutoff excludes reference points below \(0.10\max(D_{\mathrm{ref}})\).
Blank map pixels are NaN: reference points below the cutoff, which gamma does not analyse. Every analysed point receives a value. Earlier releases could also leave points that the search failed to reach as NaN, so the check below still compares the number of finite results with the number eligible above the cutoff. A finite value alone is not a pass. The colour scale ends at 2, and with max_gamma=2 larger computed values are also reported as 2. The maps show one slice; the quoted pass rate below covers the whole reference volume.
previous: 4347/4347 eligible points have finite gamma
current: 4347/4347 eligible points have finite gamma
With the current conversion, gamma of a grid against its own crop is zero at every eligible reference point. With the previous one, the artificial 10 mm shift fails the high-gradient edges.
Uneven slices need a different error model#
For a feet-first grid with relative offsets \(g_k\), the true slice coordinate is \(S_z-g_k\). The former converter attached \(-S_z+g_{N-1-k}\) to slice k, giving error
This is constant only when the paired offset sums are constant, as with evenly spaced slices. For \(S_z=55\) mm and offsets \((0,2.5,6)\) mm, the errors are \((-104,-105,-104)\) mm: successive separations have been exchanged. Two fields with these same grid headers can therefore have different gamma values after the distortion. This example isolates the coordinate effect by using the current SciPy gamma path for both calculations.
Stored slice True z (mm) Previous z (mm) Error (mm)
0 55.0 -49.0 -104.0
1 52.5 -52.5 -105.0
2 49.0 -55.0 -104.0
Middle-slice gamma: correct coordinates 0.841; previous coordinates 1.176
3.2 Crops of every axis in every orientation#
The same comparison, summarised: for each orientation, 10 mm is removed from either end of each patient axis, and the table reports the largest change in where the conversion places the voxels that the two grids share.
Largest current shift over all 48 crops: 0 mm
A crop moved the previous placement by the width removed whenever it was taken along a reversed axis, and not otherwise. For decubitus grids, an in-plane crop made the numbers of rows and columns unequal, and the previous axes no longer fitted the dose at all.
3.3 Decreasing slice offsets#
A head first supine grid may list its slices in decreasing z. The previous conversion placed these slices correctly but returned a descending z axis, which the default gamma interpolator of 0.41.0 could not use as an evaluation grid (Appendix A). The current conversion returns an ascending axis with the slices reversed to match.
DICOM z of stored slices 0, 1, 2, 3: [30. 27. 24. 21.]
previous z axis: [30. 27. 24. 21.]
current z axis: [21. 24. 27. 30.]
stored slice at each current z: [3 2 1 0]
3.4 Absolute slice offsets#
When the first element of GridFrameOffsetVector is not zero, the vector holds absolute z coordinates, which is permitted only for ImageOrientationPatient \(= (1, 0, 0, 0, 1, 0)\), with the first element equal to \(S_z\). The previous conversion added them to \(S_z\) as if they were offsets, displacing every slice by \(S_z\).
relative offsets: DICOM z [-250. -247. -244. -241.]
previous z [-250. -247. -244. -241.]
current z [-250. -247. -244. -241.]
absolute offsets: DICOM z [-250. -247. -244. -241.]
previous z [-500. -497. -494. -491.]
current z [-250. -247. -244. -241.]
Absolute offsets with another orientation are rejected:
GridFrameOffsetVector does not start at zero, so it must hold absolute z coordinates. That form is only valid when Image Orientation (Patient) is [1, 0, 0, 0, 1, 0] and its first element matches the z value of Image Position (Patient) within 0.01 mm. Got orientation [-1.0, 0.0, 0.0, 0.0, -1.0, 0.0], first offset -250.0 and IPP z -250.0.
3.5 Single-slice RT Dose#
A planar RT Dose, such as an export for a detector-array measurement, has one slice. pydicom returns its pixel data as a two-dimensional array, and the standard allows GridFrameOffsetVector to be absent. The previous conversion raised an error without the vector; with a single-valued vector, it returned z as a scalar alongside a two-dimensional dose, which gamma and dicom_dose_interpolate then rejected. The current conversion returns a three-dimensional array with a length-one z axis at \(S_z\).
Gamma treats the length-one z axis of an evaluation dose as a fixed coordinate: it searches only along y and x within the plane, and adds each reference point’s distance from the plane in quadrature. For coplanar grids, the three-dimensional arrays returned by the conversion therefore give the same gamma as the equivalent two-dimensional arrays, and reference points off the plane keep their perpendicular distance. Both are shown below. Earlier releases searched such planes with three-dimensional shells, which met a zero-thickness plane only where rounding happened to land on it, so gamma depended on the coordinate origin (#2070); comparing two-dimensional arrays was then the workaround.
pixel_array shape: (21, 25)
previous conversion: AttributeError - 'Dataset' object has no attribute 'GridFrameOffsetVector'
current conversion: dose shape (1, 21, 25) with z axis [42.]
two-dimensional gamma, 3%/3 mm: 100.0% of points pass
gamma_3d = pymedphys.gamma(
(z, y, x), dose, (z, y, x), measured[None], 3, 3, lower_percent_dose_cutoff=10
)
np.testing.assert_allclose(gamma_3d[0], gamma_2d, rtol=0, atol=1e-12, equal_nan=True)
print("three-dimensional single-slice arrays: the same gamma as the 2D arrays")
# The same measured plane 2 mm further along z. Every reference point is then
# 2 mm from the evaluation plane, so its gamma is at least 2/3 at 3 mm. Gamma
# warns that the reference points lie outside the evaluation grid, as intended.
with warnings.catch_warnings():
warnings.filterwarnings("ignore", ".*outside the evaluation grid", UserWarning)
gamma_offset = pymedphys.gamma(
(z, y, x),
dose,
(z + 2.0, y, x),
measured[None],
3,
3,
lower_percent_dose_cutoff=10,
)
offset_valid = gamma_offset[~np.isnan(gamma_offset)]
assert offset_valid.min() >= 2 / 3 - 1e-12
print(f"evaluation plane 2 mm away: minimum gamma {offset_valid.min():.3f} (at least 2/3)")
three-dimensional single-slice arrays: the same gamma as the 2D arrays
evaluation plane 2 mm away: minimum gamma 0.669 (at least 2/3)
3.6 Supported orientations and rounding#
This implementation supports the eight transverse, axis-aligned orientations above. Other cardinal orientations, such as sagittal and coronal grids, can also have three independent patient axes, but are not supported by this conversion. A genuinely oblique grid needs a more general representation or resampling.
Direction cosines within \(10^{-4}\) of a supported orientation are snapped to cardinal directions when extracting separable axes. This approximation grows in positional effect with grid extent; the exact-placement demonstrations use exact cardinal cosines. Grid equality retains the original cosines when checking physical disagreement (section 3.7). A clearly oblique grid is rejected:
Dose grid orientation is not supported. Dose grid slices must be aligned along the superoinferior axis of patient.
3.7 Summing doses needs the same pixel mapping#
Summing RT Dose files, as the experimental Sum Coincident DICOM Doses app does, adds the stored pixel arrays element by element, which is valid only if equal indices mean equal positions. The previous check compared the axis values and the array shape. Below, a head first supine grid and a feet first decubitus left grid describe the same positions, centred on \(z = 0\), but the decubitus grid stores \(x\) along its rows. Their axis values agree, yet adding the arrays adds each dose to its own transpose.
previous check, coordinates equal: True
current check, coordinates equal: False
The current check combines origin, spacing and the original encoded direction cosines into the maximum 3D distance between corresponding voxel centres, across every pair of datasets. Matching coordinate sets alone does not justify adding raw arrays.
Purpose |
Rule used here |
|---|---|
Floating-point arithmetic |
A separate 10⁻⁹ mm allowance at the decision boundaries |
Small physical mismatch |
Accept silently through 0.01 mm; warn but still treat as coincident above 0.01 mm through 0.1 mm; reject larger mismatches |
Clinical significance |
Not determined by this geometry check |
Acceptance does not resample either grid or establish mathematical identity. The warning describes a geometric discrepancy, not a judgement that it is clinically negligible. Absolute-offset metadata consistency remains a separate 0.01 mm rule. Unlike the old relative tolerance, these physical limits do not grow with distance from the origin.
4. How gamma treats axis order#
pymedphys.gamma(axes_reference, dose_reference, axes_evaluation, dose_evaluation, ...) uses its two grids differently.
The reference grid supplies the points at which gamma is evaluated. It is used as given, in either order, and gamma is returned with the shape and index order of
dose_reference:gamma[i_z, i_y, i_x]belongs to the same position asdose_reference[i_z, i_y, i_x].The evaluation grid is interpolated on shells of increasing radius around each reference point. The default interpolator needs strictly ascending, evenly spaced axes, so gamma reverses each descending evaluation axis together with the matching dimension of the evaluation dose before the search starts. Unevenly spaced evaluation axes are passed to the SciPy interpolator instead, with a warning. An evaluation axis with a single value is a fixed coordinate rather than an axis to interpolate along: the search moves along the other axes only, and the distance from the plane, line or point is added (section 3.5).
4.1 Keep coordinates with their dose#
The opening diagram showed this operation: descending evaluation axes are reversed with the matching dose dimension. The reference axes retain their supplied order.
4.2 Storage order does not change gamma#
Below, the same evaluation dose is supplied in four storage orders. The four gamma results are identical, element for element. A reference supplied with a descending axis gives the same values in its own order, so it must be plotted against its own axes.
all four evaluation orders give identical gamma: True
reference with x descending gives the same values in its own order: True
4.3 The search is bounded by the two grids#
The search starts at the smallest distance between the reference and evaluation grids and stops beyond the largest, or at max_gamma times the largest distance threshold if that is smaller. Grids that do not overlap still have valid gamma values: below, every reference point lies at least 4 mm from the evaluation grid, so the search starts at 4 mm, and every gamma value is at least \(4/3\) at 3 mm.
closest possible distance 4 mm, farthest 30 mm
/tmp/ipykernel_1313/156867632.py:6: UserWarning: 11 of 11 analysed reference points (100.0%) lie more than one search step (0.3) outside the evaluation grid. Each is compared with the evaluation data up to the grid's edge, so these points may fail. To analyse only the region the evaluation grid covers, crop the reference grid to it.
gamma_1d = pymedphys.gamma(
5. Plotting gamma after the reordering#
Gamma is returned on the reference grid, so plot it against the reference axes: gamma[i_z, i_y, i_x] is at (z_ref[i_z], y_ref[i_y], x_ref[i_x]), exactly like dose_reference[i_z, i_y, i_x]. The example below uses a head first prone reference and a head first decubitus right evaluation of the same F, shifted by 3 mm and scaled by 2%.
To overlay gamma on the stored pixel data instead, rearrange it into stored order first. to_storage_order applies the inverse of the conversion’s transposition and reversals, using only ImageOrientationPatient and the direction of the slice offsets.
to_storage_order restores pixel_array exactly for 16 of 16 encodings
In the top row, gamma and the reference dose share the reference axes, and the failing regions follow the edges of the F that run along y, where the 3 mm shift along x matters. In the bottom left, the same gamma array is drawn over the stored pixel data by index: for head first prone the stored rows and columns both run backwards, so the gamma pattern is rotated by 180° relative to the dose. For decubitus orientations the stored array is also transposed and can have a different shape.
In practice:
keep each axis tuple with its own array, and index gamma with the reference axes;
choose planes by position, for example
k = np.argmin(np.abs(z_ref - z_target)), because the evaluation grid has its own axes and a shared index need not mean a shared position;subtract doses only where their sample positions coincide, or after interpolating onto a common grid;
to overlay gamma on
pixel_array, rearrange gamma into stored order, asto_storage_orderdoes; andtreat
invert_yaxisand similar calls as display choices: they change the view, not the positions of the voxels.
The red contour marks gamma = 1 only where gamma is finite. It does not close an artificial boundary through points excluded by the cutoff.
6. Checking an earlier calculation#
The helper below checks 0.41.0 geometry and interpolation risks using headers read with pydicom.dcmread(path, stop_before_pixels=True). It reports placement, array representability and suitability of the evaluation axes separately. Set interp_algo="scipy" if that was actually used; the default checks the custom interpolator introduced in 0.41.0. Other versions need a separate historical reproduction.
It samples four corners of every slice rather than allocating a full coordinate volume. Those corners bound placement error because displacement within each slice is affine in row and column. Every slice is retained, since an uneven-offset error can peak internally.
Headers can expose a geometric problem but cannot establish its effect on a particular dose comparison, anatomical registration or clinical decision. Even when placement is correct, uneven evaluation spacing could produce wrong interpolation in the former default implementation.
HFS grid and its crop
reference: constant placement error [0. 0. 0.] mm
evaluation: constant placement error [0. 0. 0.] mm
relative translation (evaluation minus reference): [0. 0. 0.] mm
0.41.0 pymedphys interpolation: no listed axis risk identified; this is not validation of a dose result
HFP grid and its crop
reference: constant placement error [-40. 0. 0.] mm
evaluation: constant placement error [-30. 0. 0.] mm
relative translation (evaluation minus reference): [10. 0. 0.] mm
0.41.0 pymedphys interpolation: descending evaluation axis: error, NaN or non-termination in the former default path
HFP reference, HFS evaluation
reference: constant placement error [-40. 0. 0.] mm
evaluation: constant placement error [0. 0. 0.] mm
relative translation (evaluation minus reference): [40. 0. 0.] mm
0.41.0 pymedphys interpolation: no listed axis risk identified; this is not validation of a dose result
HFDL reference, HFS evaluation
reference: transverse row/column pairing was transposed
reference: index-dependent placement error, up to 136.199 mm
evaluation: constant placement error [0. 0. 0.] mm
a common-translation cancellation cannot be established
0.41.0 pymedphys interpolation: no listed axis risk identified; this is not validation of a dose result
HFS with uneven offsets
reference: constant placement error [0. 0. 0.] mm
evaluation: constant placement error [0. 0. 0.] mm
relative translation (evaluation minus reference): [0. 0. 0.] mm
0.41.0 pymedphys interpolation: uneven evaluation spacing: the former default interpolation was incorrect
FFS with uneven offsets
reference: slice-dependent placement error, up to 105.019 mm
evaluation: slice-dependent placement error, up to 105.019 mm
a common-translation cancellation cannot be established
0.41.0 pymedphys interpolation: descending evaluation axis: error, NaN or non-termination in the former default path; uneven evaluation spacing: the former default interpolation was incorrect
Single slice with scalar offset
reference: historical conversion/array representation failed: single-slice conversion returned a scalar z axis
evaluation: historical conversion/array representation failed: single-slice conversion returned a scalar z axis
a common-translation cancellation cannot be established
interpolation not assessed: no usable historical evaluation grid
reference: constant placement error [0. 0. 0.] mm
evaluation: constant placement error [0. 0. 0.] mm
relative translation (evaluation minus reference): [0. 0. 0.] mm
0.41.0 scipy interpolation: no listed axis risk identified; this is not validation of a dose result
The HFS crop example has no placement discrepancy and evenly spaced evaluation axes. The uneven HFS example also has correct placement but exposes a separate interpolation risk. The uneven FFS example has slice-dependent errors, and the single-slice example reports an array-representation failure. These are distinct findings, not one binary verdict about an earlier result.
7. Remaining limitations#
Only the supported transverse cardinal orientations are converted; genuinely oblique geometry needs a different representation (section 3.6).
Gamma samples a finite set of search shells, so it approximates the continuous minimum to within the search resolution, the distance threshold divided by
interp_fraction. Every analysed reference point receives a value, including points outside the evaluation grid and off a single-slice evaluation plane (section 3.5). Earlier releases could report such points as NaN or with origin-dependent values (#2070).The private IEC patient option raises
NotImplementedError; IEC fixed retains its previous convention.Synthetic orientation invariance validates these examples, not anatomical registration, every vendor encoding, or clinical suitability of a gamma criterion.
Appendix A. Why descending evaluation axes failed before#
The interpolation kernels locate a point by testing axis[0] <= point <= axis[-1] and searching the axis in ascending order. On a descending axis no point passes that test, so every interpolated value is the fill value, which gamma sets to infinity:
descending axis: [inf inf inf]
ascending axis: [25. 20. 15.]
A search that never finds a finite value returned NaN when max_gamma was set, and otherwise never terminated. In version 0.41.0 most affected DICOM grids did not get that far. The previous conversion produced its reversed axes with numpy.flip, which returns a view with negative strides, and passed them to the Numba kernel alongside contiguous axes. Numba cannot iterate over a tuple of arrays with different memory layouts, so the call failed:
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
No implementation of function Function(<class 'enumerate'>) found for signature:
When all three axes were reversed, as for feet first decubitus right, the layouts matched and the search ran without finding a value. The current gamma converts every evaluation axis to a contiguous ascending array before the search, so neither failure can occur.
Further reading: gamma performance#
The background workstation study compares old/new PyMedPhys and SciPy interpolation paths within a two-hour total budget and produces a results bundle with explicit coverage limits. The recorded workstation evidence retains verified partial measurements from the earlier long run. The fixed-grid performance notebook explains the avoided work with a smaller recorded experiment. Both use ascending regular grids to separate performance from the coordinate corrections illustrated here.