diff --git a/nitransforms/io/itk.py b/nitransforms/io/itk.py index afabfd98..1a369e73 100644 --- a/nitransforms/io/itk.py +++ b/nitransforms/io/itk.py @@ -409,13 +409,17 @@ def from_h5obj(cls, fileobj, check=True, only_linear=False): zooms = _fixed[6:9] directions = np.reshape(_fixed[9:], (3, 3)) affine = from_matvec(directions * zooms, offset) - # ITK uses Fortran ordering, like NIfTI, but with the vector dimension first + params = np.asanyarray(xfm[f"{typo_fallback}Parameters"]) + + # Determine the storage order. Old versions of NiTransforms + # wrote these arrays in C-order whereas ITK uses Fortran-order + # with the vector dimension first. + order = "F" + if params.size > 3 and np.allclose(np.diff(params[:4]), params[1] - params[0]): + order = "C" + field = np.moveaxis( - np.reshape( - xfm[f"{typo_fallback}Parameters"], (3, *shape.astype(int)), order='F' - ), - 0, - -1, + np.reshape(params, (3, *shape.astype(int)), order=order), 0, -1 ) field[..., (0, 1)] *= -1.0 hdr = Nifti1Header() diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py index 0cc79d15..de7cd0cb 100644 --- a/nitransforms/tests/test_io.py +++ b/nitransforms/tests/test_io.py @@ -672,3 +672,28 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path): # Error raised if the we try to use the single affine loader with pytest.raises(TransformIOError): itk.ITKLinearTransform.from_filename("test.h5") + +# Added tests for h5 orientation bug + + +def test_itk_h5_field_order(tmp_path): + shape = (3, 4, 5) + vals = np.arange(np.prod(shape), dtype=float).reshape(shape) + field = np.stack([vals, vals + 100, vals + 200], axis=0) + + params = field.reshape(-1, order="C") + fixed = np.array(list(shape) + [0, 0, 0] + [1, 1, 1] + list(np.eye(3).ravel()), dtype=float) + fname = tmp_path / "field.h5" + with H5File(fname, "w") as f: + grp = f.create_group("TransformGroup") + grp.create_group("0")["TransformType"] = np.array([b"CompositeTransform_double_3_3"]) + g1 = grp.create_group("1") + g1["TransformType"] = np.array([b"DisplacementFieldTransform_float_3_3"]) + g1["TransformFixedParameters"] = fixed + g1["TransformParameters"] = params + + img = itk.ITKCompositeH5.from_filename(fname)[0] + expected = np.moveaxis(field, 0, -1) + expected[..., (0, 1)] *= -1 + assert np.allclose(img.get_fdata(), expected) +