Python image['spacing'] axis-order conflict between ITK (z,y,x) and SimpleITK (x,y,z) blocks generic interop (PR #6021)

[hjmjohnson]

ITK and SimpleITK expose the same string keys (image['spacing'], ['origin'], ['direction']) with opposite axis-order conventions — ITK returns numpy (z,y,x) order, SimpleITK returns physical (x,y,z) order — so any generic reader that consumes these keys silently corrupts geometry. This blocks PR #6021 (SimpleITK↔ITK conversion) and cannot be fixed by picking one order, because both conventions have legitimate user expectations. Proposed fix: deprecate the ambiguous bare keys and add an explicit, order-annotated geometry protocol shared by both toolkits.

Problem detail and replication

For a 3-D image with spacing (1, 2, 3) in x,y,z (itk 5.4.6, SimpleITK 2.5.5):

Interface	['spacing'] returns	Order
itk_image['spacing']	[3. 2. 1.]	numpy (z,y,x) — flipped
sitk_image['spacing']	(1.0, 2.0, 3.0)	ITK (x,y,z)
itk.dict_from_image(img)['spacing']	(1.0, 2.0, 3.0)	ITK (x,y,z)
GetSpacing() (both toolkits)	(1.0, 2.0, 3.0)	ITK (x,y,z)

• ITK is also internally inconsistent: dict(image)['spacing'] ≠ itk.dict_from_image(image)['spacing'].
• The flip is intentional (d8086fa, 2020) for consistency with np.array(image) [z,y,x] indexing, implemented in Wrapping/Generators/Python/PyBase/pyBase.i (keys/__getitem__/__setitem__).
• The 'direction' flip (np.flip(matrix, axis=None) = D[::-1,::-1] = P·D·P) is internally correct: with origin/spacing/index all listed z,y,x, origin_zyx + D_zyx @ (spacing_zyx * index_zyx) reproduces flip(TransformIndexToPhysicalPoint(index)) to machine precision (verified, 200 random indices, non-axis-aligned 3-D rotation), and __setitem__ round-trips exactly. The corruption arises only cross-convention: passing itk_image['direction'] into an x,y,z consumer (e.g. SetDirection() or SimpleITK) yields a wrong matrix with no shape error.
• Replication script: https://gist.github.com/hjmjohnson/799940ed88d2b9616239ca417f0e9986

Both expectations are individually reasonable: numpy-ecosystem users expect spacing[0] to match axis 0 of np.array(image); ITK/SimpleITK users expect physical x,y,z order matching GetSpacing(). The defect is that a bare string key carries no order semantics, so each community bound its own convention to the same name.

Proposed solution: explicit spatial-image protocol

1. Define a small documented protocol both toolkits implement, with order-explicit names:
   • spacing_xyz, origin_xyz, direction_xyz — physical (ITK) order
   • spacing_zyx, origin_zyx, direction_zyx — numpy axis order; direction_zyx is well-defined as P·D·P (= np.flip(D, axis=None)), the correct change-of-basis when physical coordinates are also listed z,y,x (matches current image['direction'] behavior, verified against TransformIndexToPhysicalPoint).
   • Pixel access via shaped __array__ / GetArrayViewFromImage.
2. Deprecate the bare 'spacing'/'origin'/'direction' keys in ITK's __getitem__/__setitem__/keys() with a DeprecationWarning pointing to the explicit names. Current behavior unchanged during the deprecation window.
3. Coordinate with SimpleITK (Add .GetITKImage() method to sitk::Image class in Python SimpleITK/SimpleITK#2531) so SimpleITK 3.x adopts the same protocol before its conventions cement.
4. Rebuild PR ENH: Add SimpleITK Image to ITK Image conversion support #6021's converter on the protocol (short-term unblock: GetSpacing()-family + GetArrayViewFromImage, which also fixes the np.array(sitk_image) flat-array failure on released SimpleITK).

Out of scope / unchanged: dict_from_image/image_from_dict (itk-wasm interchange format, ITK order, frozen), array_from_image pixel bridges (numpy order, correct), xarray bridge (named dims, already unambiguous).

Related: PR #6021 (blocked by this), SimpleITK/SimpleITK#2531, d8086fa.

[hjmjohnson]

Verified (itk 5.4.6): the numpy (z,y,x) convention is a full isomorphic coordinate frame, not just a display flip — every ITK geometry operation (TransformIndexToPhysicalPoint, TransformPhysicalPointToContinuousIndex, direction composition/inversion) has an exact one-line numpy counterpart, to machine precision (≤3e-14), provided the protocol carries five quantities: origin, spacing, direction, index, size. Today's dict keys expose only three; every numpy-facing layer silently assumes index == 0, which mis-places voxels by 26 mm in the test case below.

The complete numpy-domain algebra, with φ(v) = np.flip(v), φ(M) = np.flip(M, axis=None) = P·M·P:

# forward:  φ(TransformContinuousIndexToPhysicalPoint(i))
p = origin + direction @ (spacing * i)
# inverse:  φ(TransformPhysicalPointToContinuousIndex(p))
i = inv(direction) @ (p - origin) / spacing
# np.array(image) coordinates:  a = i - index   (a=[0,...,0] is the region start)

Verification results (26 checks, all pass)

Test geometry: direction = R_z(30°)·R_x(20°), spacing (1,2,3), origin (5,10,15), size (7,9,11); 500 random points per check; run for start index (0,0,0) and (4,−3,10).

Property	Result
φ is a consistent involution on all five quantities; np.array(image).shape == φ(GetSize())	exact
Homomorphism: φ(A@B) == φ(A)@φ(B), φ(inv(A)) == inv(φ(A)) (so direction_zyx composes/inverts natively)	exact
Forward computed in numpy domain vs TransformContinuousIndexToPhysicalPoint oracle	≤7.1e-15
Forward computed in ITK domain, mapped through φ, vs numpy domain	≤1.5e-14
Inverse computed in numpy domain vs TransformPhysicalPointToContinuousIndex oracle	≤2.2e-14
Inverse computed in ITK domain, mapped through φ, vs numpy domain	≤2.9e-14
Cross-domain round trip (ITK forward → φ → numpy inverse → φ)	≤7.2e-15
Physical location of array[0,0,0] via origin + D·S·index, both domains, vs TransformIndexToPhysicalPoint(startIndex)	exact
Setter round-trip: numpy-domain writes through __setitem__ reproduce ITK members	exact

Key negative finding: with start index (4,−3,10), reconstructing geometry from the current three dict keys alone places array[0,0,0] 26.139 mm from its true physical location. The xarray bridge (xarray.py coords) and generic dict-key converters share this blind spot; dict_from_image is the only layer carrying bufferedRegion.index.

Also corrects an earlier claim in this issue's history: np.flip(direction, axis=None) is not the transpose in general 3-D (that was a 2-D coincidence); it is P·D·P, the correct conjugation.

Implications for the protocol design

1. The protocol must include index_xyz/index_zyx and size_xyz/size_zyx alongside the three geometry keys — without index, physical/index mappings are wrong for any streamed/cropped/non-zero-origin-region image.
2. direction_zyx is well-defined and should be included: conjugation by P preserves the full group structure, so numpy-side consumers can compose and invert directions without leaving their domain.
3. The verification script (forward, inverse, both-domain oracles, non-zero-index cases) should land as the regression test in Wrapping/Generators/Python/Tests/ with the protocol PR.