From 37cc4f99b44ab25f0cf5961c07018642ced92b29 Mon Sep 17 00:00:00 2001 From: Konstantinos Ntatsis Date: Sat, 4 Feb 2023 01:59:58 +0100 Subject: [PATCH] fix issue with non-diagonal direction matrix --- src/itk_torch_affine_matrix_bridge.py | 21 ++++++++++----------- src/test_cases.py | 26 ++++++++++++++------------ 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/src/itk_torch_affine_matrix_bridge.py b/src/itk_torch_affine_matrix_bridge.py index 6af7602..8e33fc9 100644 --- a/src/itk_torch_affine_matrix_bridge.py +++ b/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,6 +148,10 @@ 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 @@ -156,10 +159,6 @@ def itk_to_monai_affine(image, matrix, translation, center_of_rotation=None, ref 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) diff --git a/src/test_cases.py b/src/test_cases.py index 512f9a4..a18fa55 100644 --- a/src/test_cases.py +++ b/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") -