Compute gradient in physical and voxel coordinate systems

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,96 @@ 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
+
+    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:
+        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)
+
+    # 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
+
+
+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))
+
+    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

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.coils.coordinates import phys_to_vox_gradient
+from shimmingtoolbox import __dir_testing__
 
 
 def test_generate_meshgrid():
@@ -33,3 +39,152 @@ 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 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 scaling 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]])
+
+    # 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 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 scaling 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
+
+    g_x, g_y, g_z = phys_gradient(img_array, affine)  # gx = -6, gy = 2, 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')
+
+    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, 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_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 = -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 = -5.66, gy_vox = -1.41
+
+    # Calculate ground truth with the original matrix
+    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 \
+           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], 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 \
+           np.all(np.isclose(gz_truth, gy_vox))