fix issue with non-diagonal direction matrix

src/itk_torch_affine_matrix_bridge.py

@@ -94,7 +94,7 @@ def compute_reference_space_affine_matrix(image, ref_image):
     ref_direction_matrix, ref_inv_direction_matrix = [m[:ndim, :ndim].numpy() for m in compute_direction_matrix(ref_image)]
 
     # Matrix calculation
-    matrix = ref_direction_matrix @ ref_spacing_matrix @ inv_direction_matrix @ inv_spacing_matrix 
+    matrix = ref_direction_matrix @ ref_spacing_matrix @ inv_spacing_matrix @ inv_direction_matrix 
 
     # Offset calculation
     pixel_offset = -1
@@ -123,8 +123,7 @@ def itk_to_monai_affine(image, matrix, translation, center_of_rotation=None, ref
                             between the center of the image and the center of rotation.
         reference_image: The coordinate space that matrix and translation were defined
                          in respect to. If not supplied, the coordinate space of image
-                         is used. Note: it does not work properly when the direction
-                         matrix is non-diagonal.
+                         is used.
         
     Returns:
         A 4x4 MONAI affine matrix.
@@ -149,17 +148,17 @@ def itk_to_monai_affine(image, matrix, translation, center_of_rotation=None, ref
         offset_matrix, inverse_offset_matrix = compute_offset_matrix(image, center_of_rotation)
         affine_matrix = inverse_offset_matrix @ affine_matrix @ offset_matrix
 
+    # Adjust direction
+    direction_matrix, inverse_direction_matrix = compute_direction_matrix(image)
+    affine_matrix = inverse_direction_matrix @ affine_matrix @ direction_matrix
+
     # Adjust based on spacing. It is required because MONAI does not update the 
     # pixel array according to the spacing after a transformation. For example,
     # a rotation of 90deg for an image with different spacing along the two axis
     # will just rotate the image array by 90deg without also scaling accordingly.
     spacing_matrix, inverse_spacing_matrix = compute_spacing_matrix(image)
     affine_matrix = inverse_spacing_matrix @ affine_matrix @ spacing_matrix 
 
-    # Adjust direction
-    direction_matrix, inverse_direction_matrix = compute_direction_matrix(image)
-    affine_matrix = inverse_direction_matrix @ affine_matrix @ direction_matrix
-
     return affine_matrix @ reference_affine_matrix
 
 def monai_to_itk_affine(image, affine_matrix, center_of_rotation=None):
@@ -178,14 +177,14 @@ def monai_to_itk_affine(image, affine_matrix, center_of_rotation=None):
     Returns:
         The ITK matrix and the translation vector.
     """
-    # Adjust direction
-    direction_matrix, inverse_direction_matrix = compute_direction_matrix(image)
-    affine_matrix = direction_matrix @ affine_matrix @ inverse_direction_matrix 
-
     # Adjust spacing
     spacing_matrix, inverse_spacing_matrix = compute_spacing_matrix(image)
     affine_matrix = spacing_matrix @ affine_matrix @ inverse_spacing_matrix 
 
+    # Adjust direction
+    direction_matrix, inverse_direction_matrix = compute_direction_matrix(image)
+    affine_matrix = direction_matrix @ affine_matrix @ inverse_direction_matrix 
+
     # Adjust offset when center of rotation is different from center of the image
     if center_of_rotation:
         offset_matrix, inverse_offset_matrix = compute_offset_matrix(image, center_of_rotation)

src/test_cases.py

@@ -163,6 +163,12 @@ def test_cyclic_conversion(filepath):
     origin = [8.10416794, 5.4831944, 0.49211025][:ndim]
     spacing = np.array([0.7, 3.2, 1.3])[:ndim]
 
+    direction = np.array([[1.02895588, 0.22791448, 0.02429561],
+                          [0.21927512, 1.28632268, -0.14932226],
+                          [0.47455613, 0.38534345, 0.98505633]],
+                          dtype=np.float64)
+    image.SetDirection(direction[:ndim, :ndim])
+
     image.SetSpacing(spacing)
     image.SetOrigin(origin)
 
@@ -186,16 +192,16 @@ def test_use_reference_space(ref_filepath, filepath):
     image.SetSpacing([1.2, 2.0, 1.7][:ndim])
     ref_image.SetSpacing([1.9, 1.5, 1.3][:ndim])
 
-    direction = np.eye(3, dtype=np.float64)
-    direction[0, 0] = 0.68
-    direction[1, 1] = 1.05
-    direction[2, 2] = 1.83
+    direction = np.array([[1.02895588, 0.22791448, 0.02429561],
+                          [0.21927512, 1.28632268, -0.14932226],
+                          [0.47455613, 0.38534345, 0.98505633]],
+                          dtype=np.float64)
     image.SetDirection(direction[:ndim, :ndim])
 
-    ref_direction = np.eye(3, dtype=np.float64)
-    ref_direction[0, 0] = 1.25
-    ref_direction[1, 1] = 0.99 
-    ref_direction[2, 2] = 1.50
+    ref_direction = np.array([[1.26032417, -0.19243174, 0.54877414],
+                              [0.31958275, 0.9543068, 0.2720827],
+                              [-0.24106769, -0.22344502, 0.9143302]],
+                              dtype=np.float64)
     ref_image.SetDirection(ref_direction[:ndim, :ndim])
 
     image.SetOrigin([57.3, 102.0, -20.9][:ndim])
@@ -224,8 +230,4 @@ def test_use_reference_space(ref_filepath, filepath):
     print("[Min, Max] ITK: [{}, {}]".format(output_array_itk.min(), output_array_itk.max()))
     print("[Min, Max] diff: [{}, {}]".format(diff_output.min(), diff_output.max()))
 
-    itk.imwrite(itk.GetImageFromArray(diff_output), "./output/diff.tif")
-    itk.imwrite(itk.GetImageFromArray(output_array_monai), "./output/output_monai.tif")
-    itk.imwrite(itk.GetImageFromArray(output_array_itk), "./output/output_itk.tif")
-