From e0bd8237c72dd6ebadbd4e7667e4de84381c51de Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Sat, 21 Nov 2020 18:58:47 -0500 Subject: [PATCH 01/10] Initial gradient api --- shimmingtoolbox/coils/coordinates.py | 45 ++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/shimmingtoolbox/coils/coordinates.py b/shimmingtoolbox/coils/coordinates.py index cb8fcd443..5f8281fa8 100644 --- a/shimmingtoolbox/coils/coordinates.py +++ b/shimmingtoolbox/coils/coordinates.py @@ -4,6 +4,7 @@ import numpy as np from nibabel.affines import apply_affine +import math def generate_meshgrid(dim, affine): @@ -31,3 +32,47 @@ def generate_meshgrid(dim, affine): for i in range(3): coord_phys[i][ix, iy, iz] = coord_phys_list[i] return coord_phys + + +def phys_gradient(data, affine): + """Calculate the gradient of ``data`` along physical coordinates defined by ``affine`` + + Args: + data (numpy.ndarray): 3d array containing data to apply gradient + affine (numpy.ndarray): 4x4 array containing affine transformation + """ + + x_vox = 0 + y_vox = 1 + z_vox = 2 + + x_vox_spacing = math.sqrt((affine[x_vox, 0] ** 2) + (affine[x_vox, 1] ** 2) + (affine[x_vox, 2] ** 2)) + y_vox_spacing = math.sqrt((affine[y_vox, 0] ** 2) + (affine[y_vox, 1] ** 2) + (affine[y_vox, 2] ** 2)) + z_vox_spacing = math.sqrt((affine[z_vox, 0] ** 2) + (affine[z_vox, 1] ** 2) + (affine[z_vox, 2] ** 2)) + + if data.shape[x_vox] != 1: + x_vox_gradient = np.gradient(data, x_vox_spacing, axis=x_vox) + else: + x_vox_gradient = np.zeros_like(data) + + if data.shape[y_vox] != 1: + y_vox_gradient = np.gradient(data, y_vox_spacing, axis=y_vox) + else: + y_vox_gradient = np.zeros_like(data) + + if data.shape[z_vox] != 1: + z_vox_gradient = np.gradient(data, z_vox_spacing, axis=z_vox) + else: + z_vox_gradient = np.zeros_like(data) + + x_gradient = ((x_vox_gradient * (affine[x_vox, 0] / x_vox_spacing)) + + (y_vox_gradient * (affine[x_vox, 1] / y_vox_spacing)) + + (z_vox_gradient * (affine[x_vox, 2] / z_vox_spacing))) + y_gradient = ((x_vox_gradient * (affine[y_vox, 0] / x_vox_spacing)) + + (y_vox_gradient * (affine[y_vox, 1] / y_vox_spacing)) + + (z_vox_gradient * (affine[y_vox, 2] / z_vox_spacing))) + z_gradient = ((x_vox_gradient * (affine[z_vox, 2] / x_vox_spacing)) + + (y_vox_gradient * (affine[z_vox, 1] / y_vox_spacing)) + + (z_vox_gradient * (affine[z_vox, 2] / z_vox_spacing))) + + return x_gradient, y_gradient, z_gradient From a5a0bbace8cca7abc04e01e1e16ba10d1417e7e2 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Sat, 21 Nov 2020 18:59:04 -0500 Subject: [PATCH 02/10] Initial tests --- test/test_coordinates.py | 53 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/test/test_coordinates.py b/test/test_coordinates.py index d8d665077..d003aad11 100644 --- a/test/test_coordinates.py +++ b/test/test_coordinates.py @@ -2,8 +2,14 @@ # -*- coding: utf-8 -* import numpy as np +import math +import os +import nibabel as nib from shimmingtoolbox.coils.coordinates import generate_meshgrid +from shimmingtoolbox.coils.coordinates import phys_gradient +from shimmingtoolbox import __dir_shimmingtoolbox__ +from shimmingtoolbox import __dir_testing__ def test_generate_meshgrid(): @@ -33,3 +39,50 @@ def test_generate_meshgrid(): [-126.9379747, -126.9379747]]])] assert(np.all(np.isclose(coord, expected))) + + +def test_phys_gradient_synt(): + """Define a previously calculated matrix (matrix was defined at 45 deg for a gx=-6, gy=2)""" + img_array = np.expand_dims(np.array([[6 * math.sqrt(2), 4 * math.sqrt(2), 2 * math.sqrt(2)], + [2 * math.sqrt(2), 0, -2 * math.sqrt(2)], + [-2 * math.sqrt(2), -4 * math.sqrt(2), -6 * math.sqrt(2)]]), -1) + + scale = np.array([[1, 0, 0], + [0, 1, 0], + [0, 0, 1]]) + deg_angle = -45 + rot = np.array([[math.cos(deg_angle * math.pi / 180), -math.sin(deg_angle * math.pi / 180), 0], + [math.sin(deg_angle * math.pi / 180), math.cos(deg_angle * math.pi / 180), 0], + [0, 0, 1]]) + + m_affine = np.dot(scale, rot) + static_affine = [0, 0, 0, 1] + + affine = np.zeros([4, 4]) + affine[:3, :3] = m_affine + affine[3, :] = static_affine + + gx, gy, gz = phys_gradient(img_array, affine) # gx = -6, gy = 2, gz = 0 + + assert np.all(np.isclose(gx, -6)) and np.all(np.isclose(gy, 2)) and np.all(np.isclose(gz, 0)) + + +def test_phys_gradient_reel(): + fname_fieldmap = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', + 'sub-example_fieldmap.nii.gz') + + nii_fieldmap = nib.load(fname_fieldmap) + + affine = nii_fieldmap.affine + fmap = nii_fieldmap.get_fdata() + + g_x, g_y, g_z = phys_gradient(fmap[..., 0], affine) + + # Test against scaled but non rotated fieldmap, should get the same results as phys_gradient + x_coord, y_coord, z_coord = generate_meshgrid(fmap[..., 0].shape, affine) + gz_truth = np.gradient(fmap[..., 0], z_coord[0, :, 0], axis=1) + gy_truth = np.gradient(fmap[..., 0], y_coord[:, 0, 0], axis=0) + + assert np.all(np.isclose(g_z, gz_truth)) and np.all(np.isclose(g_y, gy_truth)) + + From 29cbbf0a65fac6c9c917437d056c9f65a8bb5fbf Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Mon, 23 Nov 2020 21:54:24 -0500 Subject: [PATCH 03/10] Fix mistakes and bugs --- shimmingtoolbox/coils/coordinates.py | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/shimmingtoolbox/coils/coordinates.py b/shimmingtoolbox/coils/coordinates.py index 5f8281fa8..b3e146615 100644 --- a/shimmingtoolbox/coils/coordinates.py +++ b/shimmingtoolbox/coils/coordinates.py @@ -46,9 +46,9 @@ def phys_gradient(data, affine): y_vox = 1 z_vox = 2 - x_vox_spacing = math.sqrt((affine[x_vox, 0] ** 2) + (affine[x_vox, 1] ** 2) + (affine[x_vox, 2] ** 2)) - y_vox_spacing = math.sqrt((affine[y_vox, 0] ** 2) + (affine[y_vox, 1] ** 2) + (affine[y_vox, 2] ** 2)) - z_vox_spacing = math.sqrt((affine[z_vox, 0] ** 2) + (affine[z_vox, 1] ** 2) + (affine[z_vox, 2] ** 2)) + x_vox_spacing = math.sqrt((affine[0, x_vox] ** 2) + (affine[1, x_vox] ** 2) + (affine[2, x_vox] ** 2)) + y_vox_spacing = math.sqrt((affine[0, y_vox] ** 2) + (affine[1, y_vox] ** 2) + (affine[2, y_vox] ** 2)) + z_vox_spacing = math.sqrt((affine[0, z_vox] ** 2) + (affine[1, z_vox] ** 2) + (affine[2, z_vox] ** 2)) if data.shape[x_vox] != 1: x_vox_gradient = np.gradient(data, x_vox_spacing, axis=x_vox) @@ -65,14 +65,8 @@ def phys_gradient(data, affine): else: z_vox_gradient = np.zeros_like(data) - x_gradient = ((x_vox_gradient * (affine[x_vox, 0] / x_vox_spacing)) + - (y_vox_gradient * (affine[x_vox, 1] / y_vox_spacing)) + - (z_vox_gradient * (affine[x_vox, 2] / z_vox_spacing))) - y_gradient = ((x_vox_gradient * (affine[y_vox, 0] / x_vox_spacing)) + - (y_vox_gradient * (affine[y_vox, 1] / y_vox_spacing)) + - (z_vox_gradient * (affine[y_vox, 2] / z_vox_spacing))) - z_gradient = ((x_vox_gradient * (affine[z_vox, 2] / x_vox_spacing)) + - (y_vox_gradient * (affine[z_vox, 1] / y_vox_spacing)) + - (z_vox_gradient * (affine[z_vox, 2] / z_vox_spacing))) + x_gradient = ((x_vox_gradient * (affine[x_vox, 0] / x_vox_spacing)) + (y_vox_gradient * (affine[x_vox, 1] / y_vox_spacing)) + (z_vox_gradient * (affine[x_vox, 2] / z_vox_spacing))) + y_gradient = ((x_vox_gradient * (affine[y_vox, 0] / x_vox_spacing)) + (y_vox_gradient * (affine[y_vox, 1] / y_vox_spacing)) + (z_vox_gradient * (affine[y_vox, 2] / z_vox_spacing))) + z_gradient = ((x_vox_gradient * (affine[z_vox, 0] / x_vox_spacing)) + (y_vox_gradient * (affine[z_vox, 1] / y_vox_spacing)) + (z_vox_gradient * (affine[z_vox, 2] / z_vox_spacing))) return x_gradient, y_gradient, z_gradient From 74c8fcbe9e6d2e7f9f479c5cc11e99234cd918e7 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Mon, 23 Nov 2020 21:58:23 -0500 Subject: [PATCH 04/10] Update tests --- test/test_coordinates.py | 49 ++++++++++++++++++++++++++++++++++------ 1 file changed, 42 insertions(+), 7 deletions(-) diff --git a/test/test_coordinates.py b/test/test_coordinates.py index d003aad11..f37ae2333 100644 --- a/test/test_coordinates.py +++ b/test/test_coordinates.py @@ -8,7 +8,6 @@ from shimmingtoolbox.coils.coordinates import generate_meshgrid from shimmingtoolbox.coils.coordinates import phys_gradient -from shimmingtoolbox import __dir_shimmingtoolbox__ from shimmingtoolbox import __dir_testing__ @@ -46,28 +45,63 @@ def test_phys_gradient_synt(): img_array = np.expand_dims(np.array([[6 * math.sqrt(2), 4 * math.sqrt(2), 2 * math.sqrt(2)], [2 * math.sqrt(2), 0, -2 * math.sqrt(2)], [-2 * math.sqrt(2), -4 * math.sqrt(2), -6 * math.sqrt(2)]]), -1) - + # Define a scale matrix scale = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + + # Define a rotation matrix deg_angle = -45 rot = np.array([[math.cos(deg_angle * math.pi / 180), -math.sin(deg_angle * math.pi / 180), 0], [math.sin(deg_angle * math.pi / 180), math.cos(deg_angle * math.pi / 180), 0], [0, 0, 1]]) - m_affine = np.dot(scale, rot) + # Calculate affine matrix + m_affine = np.dot(rot, scale) static_affine = [0, 0, 0, 1] + affine = np.zeros([4, 4]) + affine[:3, :3] = m_affine + affine[3, :] = static_affine + + g_x, g_y, g_z = phys_gradient(img_array, affine) # gx = -6, gy = 2, gz = 0 + assert np.all(np.isclose(g_x, -6)) and np.all(np.isclose(g_y, 2)) and np.all(np.isclose(g_z, 0)) + + +def test_phys_gradient_synt_scaled(): + """Define a previously calculated matrix (matrix was defined at 45 deg for a gx=-3, gy=1)""" + img_array = np.expand_dims(np.array([[6 * math.sqrt(2), 4 * math.sqrt(2), 2 * math.sqrt(2)], + [2 * math.sqrt(2), 0, -2 * math.sqrt(2)], + [-2 * math.sqrt(2), -4 * math.sqrt(2), -6 * math.sqrt(2)]]), -1) + # Define a scale matrix + scale = np.array([[2, 0, 0], + [0, 2, 0], + [0, 0, 2]]) + + # Define a rotation matrix + deg_angle = -45 + rot = np.array([[math.cos(deg_angle * math.pi / 180), -math.sin(deg_angle * math.pi / 180), 0], + [math.sin(deg_angle * math.pi / 180), math.cos(deg_angle * math.pi / 180), 0], + [0, 0, 1]]) + + # Calculate affine matrix + m_affine = np.dot(rot, scale) + static_affine = [0, 0, 0, 1] affine = np.zeros([4, 4]) affine[:3, :3] = m_affine affine[3, :] = static_affine - gx, gy, gz = phys_gradient(img_array, affine) # gx = -6, gy = 2, gz = 0 + g_x, g_y, g_z = phys_gradient(img_array, affine) # gx = -6, gy = 2, gz = 0 - assert np.all(np.isclose(gx, -6)) and np.all(np.isclose(gy, 2)) and np.all(np.isclose(gz, 0)) + assert np.all(np.isclose(g_x, -3)) and np.all(np.isclose(g_y, 1)) and np.all(np.isclose(g_z, 0)) def test_phys_gradient_reel(): + """ + Test the function on real data at 0 degrees of rotation so a ground truth can be calculated with a simple + gradient calculation since they are parallel. The reel data adds a degree of complexity since it is a sagittal image + """ + fname_fieldmap = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', 'sub-example_fieldmap.nii.gz') @@ -80,9 +114,10 @@ def test_phys_gradient_reel(): # Test against scaled but non rotated fieldmap, should get the same results as phys_gradient x_coord, y_coord, z_coord = generate_meshgrid(fmap[..., 0].shape, affine) - gz_truth = np.gradient(fmap[..., 0], z_coord[0, :, 0], axis=1) + gx_truth = np.zeros_like(fmap[..., 0]) gy_truth = np.gradient(fmap[..., 0], y_coord[:, 0, 0], axis=0) + gz_truth = np.gradient(fmap[..., 0], z_coord[0, :, 0], axis=1) - assert np.all(np.isclose(g_z, gz_truth)) and np.all(np.isclose(g_y, gy_truth)) + assert np.all(np.isclose(g_z, gz_truth)) and np.all(np.isclose(g_y, gy_truth)) and np.all(np.isclose(g_x, gx_truth)) From b5b7c92f2148e774bfcb1a212f18130adbe9fc96 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Mon, 23 Nov 2020 22:06:45 -0500 Subject: [PATCH 05/10] Add comments --- shimmingtoolbox/coils/coordinates.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/shimmingtoolbox/coils/coordinates.py b/shimmingtoolbox/coils/coordinates.py index b3e146615..0320cbe53 100644 --- a/shimmingtoolbox/coils/coordinates.py +++ b/shimmingtoolbox/coils/coordinates.py @@ -40,16 +40,23 @@ def phys_gradient(data, affine): Args: data (numpy.ndarray): 3d array containing data to apply gradient affine (numpy.ndarray): 4x4 array containing affine transformation + + Returns + numpy.ndarray: 3D matrix containing the gradient along the x direction in the physical coordinate system + numpy.ndarray: 3D matrix containing the gradient along the y direction in the physical coordinate system + numpy.ndarray: 3D matrix containing the gradient along the z direction in the physical coordinate system """ x_vox = 0 y_vox = 1 z_vox = 2 + # Calculate the spacing along the different voxel axis x_vox_spacing = math.sqrt((affine[0, x_vox] ** 2) + (affine[1, x_vox] ** 2) + (affine[2, x_vox] ** 2)) y_vox_spacing = math.sqrt((affine[0, y_vox] ** 2) + (affine[1, y_vox] ** 2) + (affine[2, y_vox] ** 2)) z_vox_spacing = math.sqrt((affine[0, z_vox] ** 2) + (affine[1, z_vox] ** 2) + (affine[2, z_vox] ** 2)) + # Compute the gradient along the different voxel axis if data.shape[x_vox] != 1: x_vox_gradient = np.gradient(data, x_vox_spacing, axis=x_vox) else: @@ -65,8 +72,15 @@ def phys_gradient(data, affine): else: z_vox_gradient = np.zeros_like(data) - x_gradient = ((x_vox_gradient * (affine[x_vox, 0] / x_vox_spacing)) + (y_vox_gradient * (affine[x_vox, 1] / y_vox_spacing)) + (z_vox_gradient * (affine[x_vox, 2] / z_vox_spacing))) - y_gradient = ((x_vox_gradient * (affine[y_vox, 0] / x_vox_spacing)) + (y_vox_gradient * (affine[y_vox, 1] / y_vox_spacing)) + (z_vox_gradient * (affine[y_vox, 2] / z_vox_spacing))) - z_gradient = ((x_vox_gradient * (affine[z_vox, 0] / x_vox_spacing)) + (y_vox_gradient * (affine[z_vox, 1] / y_vox_spacing)) + (z_vox_gradient * (affine[z_vox, 2] / z_vox_spacing))) + # Compute the gradient along the physical axis + x_gradient = ((x_vox_gradient * (affine[0, x_vox] / x_vox_spacing)) + + (y_vox_gradient * (affine[0, y_vox] / y_vox_spacing)) + + (z_vox_gradient * (affine[0, z_vox] / z_vox_spacing))) + y_gradient = ((x_vox_gradient * (affine[1, x_vox] / x_vox_spacing)) + + (y_vox_gradient * (affine[1, y_vox] / y_vox_spacing)) + + (z_vox_gradient * (affine[1, z_vox] / z_vox_spacing))) + z_gradient = ((x_vox_gradient * (affine[2, x_vox] / x_vox_spacing)) + + (y_vox_gradient * (affine[2, y_vox] / y_vox_spacing)) + + (z_vox_gradient * (affine[2, z_vox] / z_vox_spacing))) return x_gradient, y_gradient, z_gradient From 34b873cec5fc32899ccf9f4f714fd7ad085ec0c4 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Mon, 23 Nov 2020 22:09:35 -0500 Subject: [PATCH 06/10] Add comments to tests --- test/test_coordinates.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/test/test_coordinates.py b/test/test_coordinates.py index f37ae2333..e6bee19c3 100644 --- a/test/test_coordinates.py +++ b/test/test_coordinates.py @@ -41,11 +41,11 @@ def test_generate_meshgrid(): def test_phys_gradient_synt(): - """Define a previously calculated matrix (matrix was defined at 45 deg for a gx=-6, gy=2)""" + """Define a previously calculated matrix (matrix was defined at 45 deg for gx=-6, gy=2)""" img_array = np.expand_dims(np.array([[6 * math.sqrt(2), 4 * math.sqrt(2), 2 * math.sqrt(2)], [2 * math.sqrt(2), 0, -2 * math.sqrt(2)], [-2 * math.sqrt(2), -4 * math.sqrt(2), -6 * math.sqrt(2)]]), -1) - # Define a scale matrix + # Define a scaling matrix scale = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) @@ -69,11 +69,11 @@ def test_phys_gradient_synt(): def test_phys_gradient_synt_scaled(): - """Define a previously calculated matrix (matrix was defined at 45 deg for a gx=-3, gy=1)""" + """Define a previously calculated matrix (matrix was defined at 45 deg for gx=-3, gy=1)""" img_array = np.expand_dims(np.array([[6 * math.sqrt(2), 4 * math.sqrt(2), 2 * math.sqrt(2)], [2 * math.sqrt(2), 0, -2 * math.sqrt(2)], [-2 * math.sqrt(2), -4 * math.sqrt(2), -6 * math.sqrt(2)]]), -1) - # Define a scale matrix + # Define a scaling matrix scale = np.array([[2, 0, 0], [0, 2, 0], [0, 0, 2]]) @@ -112,12 +112,10 @@ def test_phys_gradient_reel(): g_x, g_y, g_z = phys_gradient(fmap[..., 0], affine) - # Test against scaled but non rotated fieldmap, should get the same results as phys_gradient + # Test against scaled, non rotated sagittal fieldmap, this should get the same results as phys_gradient x_coord, y_coord, z_coord = generate_meshgrid(fmap[..., 0].shape, affine) gx_truth = np.zeros_like(fmap[..., 0]) gy_truth = np.gradient(fmap[..., 0], y_coord[:, 0, 0], axis=0) gz_truth = np.gradient(fmap[..., 0], z_coord[0, :, 0], axis=1) - assert np.all(np.isclose(g_z, gz_truth)) and np.all(np.isclose(g_y, gy_truth)) and np.all(np.isclose(g_x, gx_truth)) - - + assert np.all(np.isclose(g_x, gx_truth)) and np.all(np.isclose(g_y, gy_truth)) and np.all(np.isclose(g_z, gz_truth)) From 5ff7ba172d37c3028c3a88e941ce5a19bb7bbb3c Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 25 Nov 2020 12:21:50 -0500 Subject: [PATCH 07/10] Add initial phys_to_vox_gradient --- shimmingtoolbox/coils/coordinates.py | 109 ++++++++++++++++++++++++--- 1 file changed, 100 insertions(+), 9 deletions(-) diff --git a/shimmingtoolbox/coils/coordinates.py b/shimmingtoolbox/coils/coordinates.py index 0320cbe53..883db7ed4 100644 --- a/shimmingtoolbox/coils/coordinates.py +++ b/shimmingtoolbox/coils/coordinates.py @@ -73,14 +73,105 @@ def phys_gradient(data, affine): z_vox_gradient = np.zeros_like(data) # Compute the gradient along the physical axis - x_gradient = ((x_vox_gradient * (affine[0, x_vox] / x_vox_spacing)) + - (y_vox_gradient * (affine[0, y_vox] / y_vox_spacing)) + - (z_vox_gradient * (affine[0, z_vox] / z_vox_spacing))) - y_gradient = ((x_vox_gradient * (affine[1, x_vox] / x_vox_spacing)) + - (y_vox_gradient * (affine[1, y_vox] / y_vox_spacing)) + - (z_vox_gradient * (affine[1, z_vox] / z_vox_spacing))) - z_gradient = ((x_vox_gradient * (affine[2, x_vox] / x_vox_spacing)) + - (y_vox_gradient * (affine[2, y_vox] / y_vox_spacing)) + - (z_vox_gradient * (affine[2, z_vox] / z_vox_spacing))) + x_gradient = (x_vox_gradient * affine[0, x_vox] / x_vox_spacing) + \ + (y_vox_gradient * affine[0, y_vox] / y_vox_spacing) + \ + (z_vox_gradient * affine[0, z_vox] / z_vox_spacing) + y_gradient = (x_vox_gradient * affine[1, x_vox] / x_vox_spacing) + \ + (y_vox_gradient * affine[1, y_vox] / y_vox_spacing) + \ + (z_vox_gradient * affine[1, z_vox] / z_vox_spacing) + z_gradient = (x_vox_gradient * affine[2, x_vox] / x_vox_spacing) + \ + (y_vox_gradient * affine[2, y_vox] / y_vox_spacing) + \ + (z_vox_gradient * affine[2, z_vox] / z_vox_spacing) return x_gradient, y_gradient, z_gradient + + +def phys_to_vox_gradient(gx, gy, gz, affine): + """ + Calculate the gradient along the voxel coordinates defined by ``affine`` with gradients in the physical + coordinate system + + Args: + gx (numpy.ndarray): 3D matrix containing the gradient along the x direction in the physical coordinate system + gy (numpy.ndarray): 3D matrix containing the gradient along the y direction in the physical coordinate system + gz (numpy.ndarray): 3D matrix containing the gradient along the z direction in the physical coordinate system + affine (numpy.ndarray): 4x4 array containing affine transformation + + Returns: + numpy.ndarray: 3D matrix containing the gradient along the x direction in the voxel coordinate system + numpy.ndarray: 3D matrix containing the gradient along the y direction in the voxel coordinate system + numpy.ndarray: 3D matrix containing the gradient along the z direction in the voxel coordinate system + """ + + x_vox = 0 + y_vox = 1 + z_vox = 2 + + # Calculate the spacing along the different voxel axis + x_vox_spacing = math.sqrt((affine[0, x_vox] ** 2) + (affine[1, x_vox] ** 2) + (affine[2, x_vox] ** 2)) + y_vox_spacing = math.sqrt((affine[0, y_vox] ** 2) + (affine[1, y_vox] ** 2) + (affine[2, y_vox] ** 2)) + z_vox_spacing = math.sqrt((affine[0, z_vox] ** 2) + (affine[1, z_vox] ** 2) + (affine[2, z_vox] ** 2)) + + x_coord, y_coord, z_coord = generate_meshgrid(gx.shape, affine) + + x_vox_is_neg = False + y_vox_is_neg = False + z_vox_is_neg = False + # if x_coord.shape[0] > 1: + # if (x_coord[1, 0, 0] - x_coord[0, 0, 0]) < 0: + # x_vox_is_neg = True + # if (y_coord[1, 0, 0] - y_coord[0, 0, 0]) < 0: + # x_vox_is_neg = True + # if (z_coord[1, 0, 0] - z_coord[0, 0, 0]) < 0: + # x_vox_is_neg = True + # if x_coord.shape[1] > 1: + # if (x_coord[0, 1, 0] - x_coord[0, 0, 0]) < 0: + # y_vox_is_neg = True + # if (y_coord[0, 1, 0] - y_coord[0, 0, 0]) < 0: + # y_vox_is_neg = True + # if (z_coord[0, 1, 0] - z_coord[0, 0, 0]) < 0: + # y_vox_is_neg = True + # if x_coord.shape[2] > 1: + # if (x_coord[0, 0, 1] - x_coord[0, 0, 0]) < 0: + # z_vox_is_neg = True + # if (y_coord[0, 0, 1] - y_coord[0, 0, 0]) < 0: + # z_vox_is_neg = True + # if (z_coord[0, 0, 1] - z_coord[0, 0, 0]) < 0: + # z_vox_is_neg = True + + # if x_coord.shape[0] > 1: + # x_point_0 = math.sqrt((x_coord[0, 0, 0] ** 2) + (y_coord[0, 0, 0] ** 2) + (z_coord[0, 0, 0] ** 2)) + # x_point_1 = math.sqrt((x_coord[1, 0, 0] ** 2) + (y_coord[1, 0, 0] ** 2) + (z_coord[1, 0, 0] ** 2)) + # if (x_point_1 - x_point_0) < 0: + # x_vox_is_neg = True + # if x_coord.shape[1] > 1: + # y_point_0 = math.sqrt((x_coord[0, 0, 0] ** 2) + (y_coord[0, 0, 0] ** 2) + (z_coord[0, 0, 0] ** 2)) + # y_point_1 = math.sqrt((x_coord[0, 1, 0] ** 2) + (y_coord[0, 1, 0] ** 2) + (z_coord[0, 1, 0] ** 2)) + # if (y_point_1 - y_point_0) < 0: + # y_vox_is_neg = True + # if x_coord.shape[2] > 1: + # z_point_0 = math.sqrt((x_coord[0, 0, 0] ** 2) + (y_coord[0, 0, 0] ** 2) + (z_coord[0, 0, 0] ** 2)) + # z_point_1 = math.sqrt((x_coord[0, 0, 1] ** 2) + (y_coord[0, 0, 1] ** 2) + (z_coord[0, 0, 1] ** 2)) + # if (z_point_1 - z_point_0) < 0: + # z_vox_is_neg = True + + if x_vox_is_neg: + x_vox_spacing *= -1 + if y_vox_is_neg: + y_vox_spacing *= -1 + if z_vox_is_neg: + z_vox_spacing *= -1 + + inv_affine = np.linalg.inv(affine[:3, :3]) + + gx_vox = (gx * inv_affine[0, x_vox] * x_vox_spacing) + \ + (gy * inv_affine[0, y_vox] * x_vox_spacing) + \ + (gz * inv_affine[0, z_vox] * x_vox_spacing) + gy_vox = (gx * inv_affine[1, x_vox] * y_vox_spacing) + \ + (gy * inv_affine[1, y_vox] * y_vox_spacing) + \ + (gz * inv_affine[1, z_vox] * y_vox_spacing) + gz_vox = (gx * inv_affine[2, x_vox] * z_vox_spacing) + \ + (gy * inv_affine[2, y_vox] * z_vox_spacing) + \ + (gz * inv_affine[2, z_vox] * z_vox_spacing) + + return gx_vox, gy_vox, gz_vox From eacbd74003acec27ed38eb7fb5de1ca9e57cb50d Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 25 Nov 2020 12:22:00 -0500 Subject: [PATCH 08/10] Add testing --- test/test_coordinates.py | 70 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/test/test_coordinates.py b/test/test_coordinates.py index e6bee19c3..9668e7e72 100644 --- a/test/test_coordinates.py +++ b/test/test_coordinates.py @@ -8,6 +8,7 @@ from shimmingtoolbox.coils.coordinates import generate_meshgrid from shimmingtoolbox.coils.coordinates import phys_gradient +from shimmingtoolbox.coils.coordinates import phys_to_vox_gradient from shimmingtoolbox import __dir_testing__ @@ -119,3 +120,72 @@ def test_phys_gradient_reel(): gz_truth = np.gradient(fmap[..., 0], z_coord[0, :, 0], axis=1) assert np.all(np.isclose(g_x, gx_truth)) and np.all(np.isclose(g_y, gy_truth)) and np.all(np.isclose(g_z, gz_truth)) + + +def test_phys_to_vox_gradient_synt(): + """Define a previously calculated matrix""" + img_array = np.expand_dims(np.array([[6 * math.sqrt(2), 4 * math.sqrt(2), 2 * math.sqrt(2)], + [2 * math.sqrt(2), 0, -2 * math.sqrt(2)], + [-2 * math.sqrt(2), -4 * math.sqrt(2), -6 * math.sqrt(2)]]), -1) + # Define a scaling matrix + x_vox_spacing = 1 + y_vox_spacing = -2 + scale = np.array([[x_vox_spacing, 0, 0], + [0, y_vox_spacing, 0], + [0, 0, 1]]) + + # Define a rotation matrix + deg_angle = 10 + rot = np.array([[math.cos(deg_angle * math.pi / 180), -math.sin(deg_angle * math.pi / 180), 0], + [math.sin(deg_angle * math.pi / 180), math.cos(deg_angle * math.pi / 180), 0], + [0, 0, 1]]) + + # Calculate affine matrix + m_affine = np.dot(rot, scale) + static_affine = [0, 0, 0, 1] + affine = np.zeros([4, 4]) + affine[:3, :3] = m_affine + affine[3, :] = static_affine + + gx_phys, gy_phys, gz_phys = phys_gradient(img_array, affine) + + gx_vox, gy_vox, gz_vox = phys_to_vox_gradient(gx_phys, gy_phys, gz_phys, affine) + + # Calculate ground truth with the original matrix + # TODO: account for scaling being negative + gx_truth = np.gradient(img_array, x_vox_spacing, axis=0) + gy_truth = np.gradient(img_array, y_vox_spacing, axis=1) + gz_truth = np.zeros_like(img_array) + + assert np.all(np.isclose(gx_vox, gx_truth)) and \ + np.all(np.isclose(gy_vox, gy_truth)) and \ + np.all(np.isclose(gz_vox, gz_truth)) + + +def test_phys_to_vox_gradient_reel(): + """ + Test the function on real data at 0 degrees of rotation so a ground truth can be calculated with a simple + gradient calculation since they are parallel. The reel data adds a degree of complexity since it is a sagittal image + """ + + fname_fieldmap = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', + 'sub-example_fieldmap.nii.gz') + + nii_fieldmap = nib.load(fname_fieldmap) + + affine = nii_fieldmap.affine + fmap = nii_fieldmap.get_fdata() + + gx_phys, gy_phys, gz_phys = phys_gradient(fmap[..., 0], affine) + + gx_vox, gy_vox, gz_vox = phys_to_vox_gradient(gx_phys, gy_phys, gz_phys, affine) + + # Test against scaled, non rotated sagittal fieldmap, this should get the same results as phys_gradient + x_coord, y_coord, z_coord = generate_meshgrid(fmap[..., 0].shape, affine) + gx_truth = np.zeros_like(fmap[..., 0]) + gy_truth = np.gradient(fmap[..., 0], y_coord[:, 0, 0], axis=0) + gz_truth = np.gradient(fmap[..., 0], z_coord[0, :, 0], axis=1) + + assert np.all(np.isclose(gx_truth, gz_vox)) and \ + np.all(np.isclose(gy_truth, gx_vox)) and \ + np.all(np.isclose(gz_truth, gy_vox)) From 53d0f6c95c76d9e0456a0c96438ce853839d5f2c Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Thu, 26 Nov 2020 12:44:38 -0500 Subject: [PATCH 09/10] Start solving for rotation matrix --- shimmingtoolbox/coils/coordinates.py | 45 ++++++++++++++++------------ 1 file changed, 26 insertions(+), 19 deletions(-) diff --git a/shimmingtoolbox/coils/coordinates.py b/shimmingtoolbox/coils/coordinates.py index 883db7ed4..5437d4cae 100644 --- a/shimmingtoolbox/coils/coordinates.py +++ b/shimmingtoolbox/coils/coordinates.py @@ -114,30 +114,34 @@ def phys_to_vox_gradient(gx, gy, gz, affine): x_coord, y_coord, z_coord = generate_meshgrid(gx.shape, affine) - x_vox_is_neg = False - y_vox_is_neg = False - z_vox_is_neg = False + # Solve for rotation matrix + Rx = np.zeros([3, 3]) + + # x_vox_is_neg_x = x_vox_is_neg_y = x_vox_is_neg_z = x_vox_spacing + # y_vox_is_neg_x = y_vox_is_neg_y = y_vox_is_neg_z = y_vox_spacing + # z_vox_is_neg_x = z_vox_is_neg_y = z_vox_is_neg_z = z_vox_spacing + # if x_coord.shape[0] > 1: # if (x_coord[1, 0, 0] - x_coord[0, 0, 0]) < 0: - # x_vox_is_neg = True + # x_vox_is_neg_x *= -1 # if (y_coord[1, 0, 0] - y_coord[0, 0, 0]) < 0: - # x_vox_is_neg = True + # x_vox_is_neg_y *= -1 # if (z_coord[1, 0, 0] - z_coord[0, 0, 0]) < 0: - # x_vox_is_neg = True + # x_vox_is_neg_z *= -1 # if x_coord.shape[1] > 1: # if (x_coord[0, 1, 0] - x_coord[0, 0, 0]) < 0: - # y_vox_is_neg = True + # y_vox_is_neg_x *= -1 # if (y_coord[0, 1, 0] - y_coord[0, 0, 0]) < 0: - # y_vox_is_neg = True + # y_vox_is_neg_y *= -1 # if (z_coord[0, 1, 0] - z_coord[0, 0, 0]) < 0: - # y_vox_is_neg = True + # y_vox_is_neg_z *= -1 # if x_coord.shape[2] > 1: # if (x_coord[0, 0, 1] - x_coord[0, 0, 0]) < 0: - # z_vox_is_neg = True + # z_vox_is_neg_x *= -1 # if (y_coord[0, 0, 1] - y_coord[0, 0, 0]) < 0: - # z_vox_is_neg = True + # z_vox_is_neg_y *= -1 # if (z_coord[0, 0, 1] - z_coord[0, 0, 0]) < 0: - # z_vox_is_neg = True + # z_vox_is_neg_z *= -1 # if x_coord.shape[0] > 1: # x_point_0 = math.sqrt((x_coord[0, 0, 0] ** 2) + (y_coord[0, 0, 0] ** 2) + (z_coord[0, 0, 0] ** 2)) @@ -155,13 +159,6 @@ def phys_to_vox_gradient(gx, gy, gz, affine): # if (z_point_1 - z_point_0) < 0: # z_vox_is_neg = True - if x_vox_is_neg: - x_vox_spacing *= -1 - if y_vox_is_neg: - y_vox_spacing *= -1 - if z_vox_is_neg: - z_vox_spacing *= -1 - inv_affine = np.linalg.inv(affine[:3, :3]) gx_vox = (gx * inv_affine[0, x_vox] * x_vox_spacing) + \ @@ -174,4 +171,14 @@ def phys_to_vox_gradient(gx, gy, gz, affine): (gy * inv_affine[2, y_vox] * z_vox_spacing) + \ (gz * inv_affine[2, z_vox] * z_vox_spacing) + # gx_vox = (gx * inv_affine[0, x_vox] * x_vox_spacing) + \ + # (gy * inv_affine[0, y_vox] * x_vox_spacing) + \ + # (gz * inv_affine[0, z_vox] * x_vox_spacing) + # gy_vox = (gx * inv_affine[1, x_vox] * y_vox_spacing) + \ + # (gy * inv_affine[1, y_vox] * y_vox_spacing) + \ + # (gz * inv_affine[1, z_vox] * y_vox_spacing) + # gz_vox = (gx * inv_affine[2, x_vox] * z_vox_spacing) + \ + # (gy * inv_affine[2, y_vox] * z_vox_spacing) + \ + # (gz * inv_affine[2, z_vox] * z_vox_spacing) + return gx_vox, gy_vox, gz_vox From c977ab49ca491920a4152790b456650e1a6ab872 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Fri, 27 Nov 2020 21:01:49 -0500 Subject: [PATCH 10/10] Remove unecessary comments and fix phys_to_vox --- shimmingtoolbox/coils/coordinates.py | 57 ---------------------------- test/test_coordinates.py | 15 ++++---- 2 files changed, 7 insertions(+), 65 deletions(-) diff --git a/shimmingtoolbox/coils/coordinates.py b/shimmingtoolbox/coils/coordinates.py index 5437d4cae..d8bcefdc8 100644 --- a/shimmingtoolbox/coils/coordinates.py +++ b/shimmingtoolbox/coils/coordinates.py @@ -112,53 +112,6 @@ def phys_to_vox_gradient(gx, gy, gz, affine): y_vox_spacing = math.sqrt((affine[0, y_vox] ** 2) + (affine[1, y_vox] ** 2) + (affine[2, y_vox] ** 2)) z_vox_spacing = math.sqrt((affine[0, z_vox] ** 2) + (affine[1, z_vox] ** 2) + (affine[2, z_vox] ** 2)) - x_coord, y_coord, z_coord = generate_meshgrid(gx.shape, affine) - - # Solve for rotation matrix - Rx = np.zeros([3, 3]) - - # x_vox_is_neg_x = x_vox_is_neg_y = x_vox_is_neg_z = x_vox_spacing - # y_vox_is_neg_x = y_vox_is_neg_y = y_vox_is_neg_z = y_vox_spacing - # z_vox_is_neg_x = z_vox_is_neg_y = z_vox_is_neg_z = z_vox_spacing - - # if x_coord.shape[0] > 1: - # if (x_coord[1, 0, 0] - x_coord[0, 0, 0]) < 0: - # x_vox_is_neg_x *= -1 - # if (y_coord[1, 0, 0] - y_coord[0, 0, 0]) < 0: - # x_vox_is_neg_y *= -1 - # if (z_coord[1, 0, 0] - z_coord[0, 0, 0]) < 0: - # x_vox_is_neg_z *= -1 - # if x_coord.shape[1] > 1: - # if (x_coord[0, 1, 0] - x_coord[0, 0, 0]) < 0: - # y_vox_is_neg_x *= -1 - # if (y_coord[0, 1, 0] - y_coord[0, 0, 0]) < 0: - # y_vox_is_neg_y *= -1 - # if (z_coord[0, 1, 0] - z_coord[0, 0, 0]) < 0: - # y_vox_is_neg_z *= -1 - # if x_coord.shape[2] > 1: - # if (x_coord[0, 0, 1] - x_coord[0, 0, 0]) < 0: - # z_vox_is_neg_x *= -1 - # if (y_coord[0, 0, 1] - y_coord[0, 0, 0]) < 0: - # z_vox_is_neg_y *= -1 - # if (z_coord[0, 0, 1] - z_coord[0, 0, 0]) < 0: - # z_vox_is_neg_z *= -1 - - # if x_coord.shape[0] > 1: - # x_point_0 = math.sqrt((x_coord[0, 0, 0] ** 2) + (y_coord[0, 0, 0] ** 2) + (z_coord[0, 0, 0] ** 2)) - # x_point_1 = math.sqrt((x_coord[1, 0, 0] ** 2) + (y_coord[1, 0, 0] ** 2) + (z_coord[1, 0, 0] ** 2)) - # if (x_point_1 - x_point_0) < 0: - # x_vox_is_neg = True - # if x_coord.shape[1] > 1: - # y_point_0 = math.sqrt((x_coord[0, 0, 0] ** 2) + (y_coord[0, 0, 0] ** 2) + (z_coord[0, 0, 0] ** 2)) - # y_point_1 = math.sqrt((x_coord[0, 1, 0] ** 2) + (y_coord[0, 1, 0] ** 2) + (z_coord[0, 1, 0] ** 2)) - # if (y_point_1 - y_point_0) < 0: - # y_vox_is_neg = True - # if x_coord.shape[2] > 1: - # z_point_0 = math.sqrt((x_coord[0, 0, 0] ** 2) + (y_coord[0, 0, 0] ** 2) + (z_coord[0, 0, 0] ** 2)) - # z_point_1 = math.sqrt((x_coord[0, 0, 1] ** 2) + (y_coord[0, 0, 1] ** 2) + (z_coord[0, 0, 1] ** 2)) - # if (z_point_1 - z_point_0) < 0: - # z_vox_is_neg = True - inv_affine = np.linalg.inv(affine[:3, :3]) gx_vox = (gx * inv_affine[0, x_vox] * x_vox_spacing) + \ @@ -171,14 +124,4 @@ def phys_to_vox_gradient(gx, gy, gz, affine): (gy * inv_affine[2, y_vox] * z_vox_spacing) + \ (gz * inv_affine[2, z_vox] * z_vox_spacing) - # gx_vox = (gx * inv_affine[0, x_vox] * x_vox_spacing) + \ - # (gy * inv_affine[0, y_vox] * x_vox_spacing) + \ - # (gz * inv_affine[0, z_vox] * x_vox_spacing) - # gy_vox = (gx * inv_affine[1, x_vox] * y_vox_spacing) + \ - # (gy * inv_affine[1, y_vox] * y_vox_spacing) + \ - # (gz * inv_affine[1, z_vox] * y_vox_spacing) - # gz_vox = (gx * inv_affine[2, x_vox] * z_vox_spacing) + \ - # (gy * inv_affine[2, y_vox] * z_vox_spacing) + \ - # (gz * inv_affine[2, z_vox] * z_vox_spacing) - return gx_vox, gy_vox, gz_vox diff --git a/test/test_coordinates.py b/test/test_coordinates.py index 9668e7e72..50123aaf4 100644 --- a/test/test_coordinates.py +++ b/test/test_coordinates.py @@ -135,7 +135,7 @@ def test_phys_to_vox_gradient_synt(): [0, 0, 1]]) # Define a rotation matrix - deg_angle = 10 + deg_angle = -10 rot = np.array([[math.cos(deg_angle * math.pi / 180), -math.sin(deg_angle * math.pi / 180), 0], [math.sin(deg_angle * math.pi / 180), math.cos(deg_angle * math.pi / 180), 0], [0, 0, 1]]) @@ -147,14 +147,13 @@ def test_phys_to_vox_gradient_synt(): affine[:3, :3] = m_affine affine[3, :] = static_affine - gx_phys, gy_phys, gz_phys = phys_gradient(img_array, affine) + gx_phys, gy_phys, gz_phys = phys_gradient(img_array, affine) # gx = -5.32, gy = 2.37, gz = 0 - gx_vox, gy_vox, gz_vox = phys_to_vox_gradient(gx_phys, gy_phys, gz_phys, affine) + gx_vox, gy_vox, gz_vox = phys_to_vox_gradient(gx_phys, gy_phys, gz_phys, affine) # gx_vox = -5.66, gy_vox = -1.41 # Calculate ground truth with the original matrix - # TODO: account for scaling being negative - gx_truth = np.gradient(img_array, x_vox_spacing, axis=0) - gy_truth = np.gradient(img_array, y_vox_spacing, axis=1) + gx_truth = np.gradient(img_array, abs(x_vox_spacing), axis=0) + gy_truth = np.gradient(img_array, abs(y_vox_spacing), axis=1) gz_truth = np.zeros_like(img_array) assert np.all(np.isclose(gx_vox, gx_truth)) and \ @@ -183,8 +182,8 @@ def test_phys_to_vox_gradient_reel(): # Test against scaled, non rotated sagittal fieldmap, this should get the same results as phys_gradient x_coord, y_coord, z_coord = generate_meshgrid(fmap[..., 0].shape, affine) gx_truth = np.zeros_like(fmap[..., 0]) - gy_truth = np.gradient(fmap[..., 0], y_coord[:, 0, 0], axis=0) - gz_truth = np.gradient(fmap[..., 0], z_coord[0, :, 0], axis=1) + gy_truth = np.gradient(fmap[..., 0], abs(y_coord[1, 0, 0] - y_coord[0, 0, 0]), axis=0) + gz_truth = np.gradient(fmap[..., 0], abs(z_coord[0, 1, 0] - z_coord[0, 0, 0]), axis=1) assert np.all(np.isclose(gx_truth, gz_vox)) and \ np.all(np.isclose(gy_truth, gx_vox)) and \