Merge 34b873c into f9cd9c9

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,55 @@ 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

test/test_coordinates.py

@@ -2,8 +2,13 @@
 # -*- 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_testing__
 
 
 def test_generate_meshgrid():
@@ -33,3 +38,84 @@ 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))