Repository: pyvista/pyvista

Issue:
Direction matrix not preserved when saving and loading ImageData objects
### Describe the bug, what's wrong, and what you expected.

When saving a PyVista `ImageData` object with a custom direction matrix to a `.vtk` file and loading it back, the `direction_matrix` is not preserved.
Hence, manipulations of the `ImageData` object such as rotations and transformations are lost after loading and saving.
Any ideas how we could resolve this?

### Steps to reproduce the bug.

Minimum working example
```python
import pyvista as pv

temp_file = "dataset.vtk"

dataset = pv.ImageData(dimensions=(2,3,4), spacing=(0.5,1.0,1.5), origin=(0., 0.5, 1.0))
dataset.SetDirectionMatrix((-1, 0, 0, 0, 1, 0, 0, 0, 1))

dataset.save(temp_file)

loaded_dataset = pv.read(temp_file)

print("Original direction matrix: \n")
print(dataset.direction_matrix)

print("Direction matrix after loading from file: \n")
print(loaded_dataset.direction_matrix)
```
Output
```
Original direction matrix: 
[[-1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]

Direction matrix after loading from file: 
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
```

### System Information

```shell
--------------------------------------------------------------------------------
  Date: Fri Apr 25 12:55:15 2025 CEST
                OS : Darwin (macOS 14.7.5)
            CPU(s) : 8
           Machine : x86_64
      Architecture : 64bit
               RAM : 16.0 GiB
       Environment : Python
       File system : apfs
        GPU Vendor : Apple
      GPU Renderer : Apple M1 Pro
       GPU Version : 4.1 Metal - 88.1
     Render Window : vtkCocoaRenderWindow
  MathText Support : True
  Python 3.11.5 (main, Aug 26 2023, 06:41:45) [Clang 16.0.3 ]
           pyvista : 0.45.0
               vtk : 9.4.2
             numpy : 2.1.1
        matplotlib : 3.10.1
            scooby : 0.10.0
             pooch : 1.8.2
            pillow : 11.1.0
           imageio : 2.37.0
             scipy : 1.15.2
              tqdm : 4.67.1
--------------------------------------------------------------------------------
```

### Screenshots

_No response_

Discussion:
user27182: After reading https://docs.vtk.org/en/latest/design_documents/VTKFileFormats.html more closely, I think the convention is that the `.vtk` is the extension is used exclusively for the legacy format, whereas `.vt*` is used exclusively for XML formats.

This was not obvious to me, and I thought these two extensions were interchangeable (and up to personal preference if you want to use `.vtk` or `.vtp` for PolyData).

So these two formats were perhaps previously not ambiguous by convention. But now if we use `.vtk` to mean either legacy or XML, reading becomes ambiguous. This ambiguity is handled fine by this PR (by using an XML parser to check first), but perhaps other programs might try reading an XML-formatted `.vtk` file using the legacy reader and throw an error (e.g. reader.GetOutput() returns None).

Perhaps, it may instead be better to revert many changes from this PR, and simply issue a warning when saving `ImageData` as `.vtk` indicating that the direction matrix will be lost?
user27182: > perhaps other programs might try reading an XML-formatted `.vtk` file using the legacy reader and throw an error (e.g. reader.GetOutput() returns None).

This is indeed the case. ParaView can't open an XML-formatted `.vtk` file.

> Perhaps, it may instead be better to revert many changes from this PR, and simply issue a warning when saving ImageData as .vtk indicating that the direction matrix will be lost?

I will do this instead.


pyvista-bot: <!-- NETLIFY DEPLOY COMMENT GENERATED BY ACTIONS_NETLIFY - APP ID SHA256: d3574f413fe51079937885070831e027ca5b10f7facb9e609494e0a481a31f79 -->
🚀 Deployed on https://680c3070c0639b28ef033bcc--pyvista-dev.netlify.app

Patch:
diff --git a/pyvista/core/dataobject.py b/pyvista/core/dataobject.py
index 30d593d3..154408a4 100644
--- a/pyvista/core/dataobject.py
+++ b/pyvista/core/dataobject.py
@@ -175,6 +175,14 @@ class DataObject(_vtk.DisableVtkSnakeCase, _vtk.vtkPyVistaOverride):
                         'Use `move_nested_field_data_to_root` to store the field data with the root MultiBlock before saving.'
                     )
 
+        def _warn_imagedata_direction_matrix(mesh: pyvista.ImageData) -> None:
+            if not np.allclose(mesh.direction_matrix, np.eye(3)):
+                warnings.warn(
+                    'The direction matrix for ImageData will not be saved using the legacy `.vtk` format.\n'
+                    'See https://gitlab.kitware.com/vtk/vtk/-/issues/19663 \n'
+                    'Use the `.vti` extension instead (XML format).'
+                )
+
         def _write_vtk(mesh_: DataObject) -> None:
             writer = mesh_._WRITERS[file_ext]()
             set_vtkwriter_mode(vtk_writer=writer, use_binary=binary)
@@ -213,6 +221,8 @@ class DataObject(_vtk.DisableVtkSnakeCase, _vtk.vtkPyVistaOverride):
         # warn if data will be lost
         if isinstance(self, pyvista.MultiBlock):
             _warn_multiblock_nested_field_data(self)
+        if isinstance(self, pyvista.ImageData) and file_ext == '.vtk':
+            _warn_imagedata_direction_matrix(self)
 
         writer_exts = self._WRITERS.keys()
         if file_ext in writer_exts:


Test patch:
diff --git a/tests/core/test_grid.py b/tests/core/test_grid.py
index 9352ec45..919629da 100644
--- a/tests/core/test_grid.py
+++ b/tests/core/test_grid.py
@@ -1083,21 +1083,33 @@ def test_save_rectilinear(extension, binary, tmpdir):
 
 @pytest.mark.parametrize('binary', [True, False])
 @pytest.mark.parametrize('extension', ['.vtk', '.vti'])
-def test_save_uniform(extension, binary, tmpdir):
-    filename = str(tmpdir.mkdir('tmpdir').join(f'tmp.{extension}'))
-    ogrid = examples.load_uniform()
-    ogrid.save(filename, binary)
-    grid = pv.ImageData(filename)
-    assert grid.n_cells == ogrid.n_cells
-    assert grid.origin == ogrid.origin
-    assert grid.spacing == ogrid.spacing
-    assert grid.dimensions == ogrid.dimensions
-    grid = pv.read(filename)
-    assert isinstance(grid, pv.ImageData)
-    assert grid.n_cells == ogrid.n_cells
-    assert grid.origin == ogrid.origin
-    assert grid.spacing == ogrid.spacing
-    assert grid.dimensions == ogrid.dimensions
+@pytest.mark.parametrize('reader', [pv.ImageData, pv.read])
+@pytest.mark.parametrize('direction_matrix', [np.eye(3), np.diag((-1, 1, -1))])
+def test_save_uniform(extension, binary, tmpdir, uniform, reader, direction_matrix):
+    filename = str(tmpdir.mkdir('tmpdir').join(f'tmp{extension}'))
+    is_identity_matrix = np.allclose(direction_matrix, np.eye(3))
+    uniform.direction_matrix = direction_matrix
+
+    if extension == '.vtk' and not is_identity_matrix:
+        match = re.escape(
+            'The direction matrix for ImageData will not be saved using the legacy `.vtk` format.\n'
+            'See https://gitlab.kitware.com/vtk/vtk/-/issues/19663 \n'
+            'Use the `.vti` extension instead (XML format).'
+        )
+        with pytest.warns(UserWarning, match=match):
+            uniform.save(filename, binary)
+    else:
+        uniform.save(filename, binary)
+
+    grid = reader(filename)
+
+    if extension == '.vtk' and not is_identity_matrix:
+        # Direction matrix is lost
+        assert not np.allclose(grid.direction_matrix, uniform.direction_matrix)
+        # Add it back manually for equality check
+        grid.direction_matrix = uniform.direction_matrix
+
+    assert grid == uniform
 
 
 def test_grid_points():