Merge 3ae81e0 into cdd1ebb

docs/source/api_reference/api.rst

@@ -56,9 +56,9 @@ Shim
 .. Error while importing sklearn and skimage
     .. automodule:: shimmingtoolbox.shim.realtime_shim
         :members:
-
-.. automodule:: shimmingtoolbox.optimizer.sequential
-   :members:
+.. Error while importing sklearn (Linear Regression)
+    .. automodule:: shimmingtoolbox.shim.sequencer
+       :members:
 
 Optimizer
 ---------

shimmingtoolbox/coils/coordinates.py

@@ -158,7 +158,7 @@ def resample_from_to(nii_from_img, nii_to_vox_map, order=2, mode='nearest', cval
 
     from_img = nii_from_img.get_fdata()
     if from_img.ndim == 2:
-        nii_from_img_3d = nib.Nifti1Image(np.expand_dims(from_img, -1), nii_from_img.affine)
+        nii_from_img_3d = nib.Nifti1Image(np.expand_dims(from_img, -1), nii_from_img.affine, header=nii_from_img.header)
         nii_resampled = nib_resample_from_to(nii_from_img_3d, nii_to_vox_map, order=order, mode=mode, cval=cval,
                                              out_class=out_class)
 
@@ -176,15 +176,15 @@ def resample_from_to(nii_from_img, nii_to_vox_map, order=2, mode='nearest', cval
         cpus = mp.cpu_count()
         if cpus == 1:
             for it in range(nt):
-                nii_from_img_3d = nib.Nifti1Image(from_img[..., it], nii_from_img.affine)
+                nii_from_img_3d = nib.Nifti1Image(from_img[..., it], nii_from_img.affine, header=nii_from_img.header)
                 nii_resampled_3d = nib_resample_from_to(nii_from_img_3d, nii_to_vox_map, order=order, mode=mode,
                                                         cval=cval, out_class=out_class)
                 resampled_4d[..., it] = nii_resampled_3d.get_fdata()
         else:
             # Create inputs for multiprocessing
             inputs = []
             for it in range(nt):
-                nii_from_img_3d = nib.Nifti1Image(from_img[..., it], nii_from_img.affine)
+                nii_from_img_3d = nib.Nifti1Image(from_img[..., it], nii_from_img.affine, header=nii_from_img.header)
                 inputs.append((nii_from_img_3d, nii_to_vox_map, order, mode, cval, out_class))
 
             with mp.Pool(cpus - 1) as pool:
@@ -194,7 +194,7 @@ def resample_from_to(nii_from_img, nii_to_vox_map, order=2, mode='nearest', cval
             for it in range(from_img.shape[3]):
                 resampled_4d[..., it] = output[it]
 
-        nii_resampled = nib.Nifti1Image(resampled_4d, nii_to_vox_map.affine)
+        nii_resampled = nib.Nifti1Image(resampled_4d, nii_to_vox_map.affine, header=nii_to_vox_map.header)
 
     else:
         raise NotImplementedError("Dimensions of input can only be 2D, 3D or 4D")

shimmingtoolbox/masking/mask_utils.py

@@ -0,0 +1,191 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+
+import nibabel as nib
+import numpy as np
+from scipy.ndimage import binary_dilation
+from scipy.ndimage import binary_opening
+from scipy.ndimage import generate_binary_structure
+from scipy.ndimage import iterate_structure
+
+from shimmingtoolbox.coils.coordinates import resample_from_to
+
+
+def resample_mask(nii_mask_from, nii_target, from_slices):
+    """
+    Select the appropriate slices from ``nii_mask_from`` using ``from_slices`` and resample onto ``nii_target``
+
+    Args:
+        nii_mask_from (nib.Nifti1Image): Mask to resample from. False or 0 signifies not included.
+        nii_target (nib.Nifti1Image): Target image to resample onto.
+        from_slices (tuple): Tuple containing the slices to select from nii_mask_from.
+
+    Returns:
+        nib.Nifti1Image: Mask resampled with nii_target.shape and nii_target.affine.
+    """
+
+    mask_from = nii_mask_from.get_fdata()
+
+    # Initialize a sliced mask and select the slices from from_slices
+    sliced_mask = np.full_like(mask_from, fill_value=False)
+    sliced_mask[:, :, from_slices] = mask_from[:, :, from_slices]
+
+    # Create nibabel object
+    nii_mask = nib.Nifti1Image(sliced_mask.astype(int), nii_mask_from.affine, header=nii_mask_from.header)
+
+    # Resample the mask onto nii_target
+    nii_mask_target = resample_from_to(nii_mask, nii_target, order=0, mode='grid-constant', cval=0)
+
+    # dilate the mask to add more pixels in particular directions
+    mask_dilated = dilate_binary_mask(nii_mask_target.get_fdata(), 'line', 3)
+    nii_mask_dilated = nib.Nifti1Image(mask_dilated, nii_mask_target.affine, header=nii_mask_target.header)
+
+    #######
+    # Debug TODO: REMOVE
+    import os
+    nib.save(nii_mask, os.path.join(os.curdir, f"fig_mask_{from_slices[0]}.nii.gz"))
+    nib.save(nii_mask_from, os.path.join(os.curdir, "fig_mask_roi.nii.gz"))
+    nib.save(nii_mask_target, os.path.join(os.curdir, f"fig_mask_res{from_slices[0]}.nii.gz"))
+    nib.save(nii_mask_dilated, os.path.join(os.curdir, f"fig_mask_dilated{from_slices[0]}.nii.gz"))
+    #######
+
+    return nii_mask_dilated
+
+
+def dilate_binary_mask(mask, shape='cross', size=3):
+    """
+    Dilates a binary mask according to different shapes and kernel size
+
+    Args:
+        mask (numpy.ndarray): 3d array containing the binary mask.
+        shape (str): 3d kernel to perform the dilation. Allowed shapes are: 'sphere', 'cross', 'line', 'cube'.
+                     'line' uses 3 line kernels to extend in each directions by "(size - 1) / 2" only if that direction
+                     is smaller than (size - 1) / 2
+        size (int): length of a side of the 3d kernel.
+
+    Returns:
+        numpy.ndarray: Dilated mask.
+
+    Notes:
+        Kernels for
+            'cross':
+                array([[[0., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 0.]],
+                       [[0., 1., 0.],
+                        [1., 1., 1.],
+                        [0., 1., 0.]],
+                       [[0., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 0.]]])
+
+            'sphere' size 5:
+                [[[False False False False False]
+                  [False False False False False]
+                  [False False  True False False]
+                  [False False False False False]
+                  [False False False False False]]
+                 [[False False False False False]
+                  [False False  True False False]
+                  [False  True  True  True False]
+                  [False False  True False False]
+                  [False False False False False]]
+                 [[False False  True False False]
+                  [False  True  True  True False]
+                  [ True  True  True  True  True]
+                  [False  True  True  True False]
+                  [False False  True False False]]
+                 [[False False False False False]
+                  [False False  True False False]
+                  [False  True  True  True False]
+                  [False False  True False False]
+                  [False False False False False]]
+                 [[False False False False False]
+                  [False False False False False]
+                  [False False  True False False]
+                  [False False False False False]
+                  [False False False False False]]]
+
+            'cube':
+                array([[[ True,  True,  True],
+                        [ True,  True,  True],
+                        [ True,  True,  True]],
+                       [[ True,  True,  True],
+                        [ True,  True,  True],
+                        [ True,  True,  True]],
+                       [[ True,  True,  True],
+                        [ True,  True,  True],
+                        [ True,  True,  True]]])
+
+            'line':
+                array([[[0., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 0.]],
+                       [[0., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 0.]],
+                       [[0., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 0.]]])
+
+    """
+
+    if size % 2 == 0 or size < 3:
+        raise ValueError("Size must be odd and greater or equal to 3")
+    # Find the middle pixel, will always work since we check size is odd
+    mid_pixel = int((size - 1) / 2)
+
+    if shape == 'sphere':
+        # Define kernel to perform the dilation
+        struct_sphere_size1 = generate_binary_structure(3, 1)
+        struct = iterate_structure(struct_sphere_size1, mid_pixel)
+
+        # Dilate
+        mask_dilated = binary_dilation(mask, structure=struct)
+
+    elif shape == 'cross':
+        struct = np.zeros([size, size, size])
+        struct[:, mid_pixel, mid_pixel] = 1
+        struct[mid_pixel, :, mid_pixel] = 1
+        struct[mid_pixel, mid_pixel, :] = 1
+
+        # Dilate
+        mask_dilated = binary_dilation(mask, structure=struct)
+
+    elif shape == 'cube':
+        # Define kernel to perform the dilation
+        struct_cube_size1 = generate_binary_structure(3, 3)
+        struct = iterate_structure(struct_cube_size1, mid_pixel)
+
+        # Dilate
+        mask_dilated = binary_dilation(mask, structure=struct)
+
+    elif shape == 'line':
+
+        struct_dim1 = np.zeros([size, size, size])
+        struct_dim1[:, mid_pixel, mid_pixel] = 1
+        # Finds where the structure fits
+        open1 = binary_opening(mask, structure=struct_dim1)
+        # Select Everything that does not fit within the structure and erode along a dim
+        dim1 = binary_dilation(np.logical_and(np.logical_not(open1), mask), structure=struct_dim1)
+
+        struct_dim2 = np.zeros([size, size, size])
+        struct_dim2[mid_pixel, :, mid_pixel] = 1
+        # Finds where the structure fits
+        open2 = binary_opening(mask, structure=struct_dim2)
+        # Select Everything that does not fit within the structure and erode along a dim
+        dim2 = binary_dilation(np.logical_and(np.logical_not(open2), mask), structure=struct_dim2)
+
+        struct_dim3 = np.zeros([size, size, size])
+        struct_dim3[mid_pixel, mid_pixel, :] = 1
+        # Finds where the structure fits
+        open3 = binary_opening(mask, structure=struct_dim3)
+        # Select Everything that does not fit within the structure and erode along a dim
+        dim3 = binary_dilation(np.logical_and(np.logical_not(open3), mask), structure=struct_dim3)
+
+        mask_dilated = np.logical_or(np.logical_or(np.logical_or(dim1, dim2), dim3), mask)
+
+    else:
+        raise ValueError("Use of non supported algorithm for dilating the mask")
+
+    return mask_dilated.astype(int)

shimmingtoolbox/optimizer/basic_optimizer.py

@@ -17,7 +17,7 @@ class Optimizer(object):
     """
     Optimizer object that stores coil profiles and optimizes an unshimmed volume given a mask.
     Use optimize(args) to optimize a given mask.
-    For basic optimizer, uses unbounded pseudo-inverse.
+    For basic optimizer, uses *unbounded* pseudo-inverse.
 
     Attributes:
         coils (ListCoil): List of Coil objects containing the coil profiles and related constraints
@@ -43,20 +43,44 @@ def __init__(self, coils: ListCoil, unshimmed, affine):
         logging.basicConfig(filename='test_optimizer.log', filemode='w', level=logging.DEBUG)
 
         self.coils = coils
+        self.unshimmed = np.array([])
+        self.unshimmed_affine = []
+        self.merged_coils = []
+        self.merged_bounds = []
+        self.set_unshimmed(unshimmed, affine)
 
+    def set_unshimmed(self, unshimmed, affine):
+        """
+        Set the unshimmed array to a new array. Resamples coil profiles accordingly.
+
+        Args:
+            unshimmed (numpy.ndarray): 3d array of unshimmed volume
+            affine: (numpy.ndarray): 4x4 array containing the qform affine transformation for the unshimmed array
+        """
         # Check dimensions of unshimmed map
         if unshimmed.ndim != 3:
             raise ValueError(f"Unshimmed profile has {unshimmed.ndim} dimensions, expected 3 (dim1, dim2, dim3)")
-        self.unshimmed = unshimmed
 
         # Check dimensions of affine
         if affine.shape != (4, 4):
             raise ValueError("Shape of affine matrix should be 4x4")
+
+        # Define coil profiles if unshimmed or affine is different than previously
+        if (self.unshimmed.shape != unshimmed.shape) or not np.all(self.unshimmed_affine == affine):
+            self.merged_coils, self.merged_bounds = self.merge_coils(unshimmed, affine)
+
+        self.unshimmed = unshimmed
         self.unshimmed_affine = affine
 
-        # Define coil profiles
-        merged_coils, merged_bounds = self.merge_coils(unshimmed, affine)
-        self.merged_coils = merged_coils
+    def set_merged_bounds(self, merged_bounds):
+        """
+        Changes the default bounds set in the coil profile
+
+        Args:
+            merged_bounds: Concatenated coil profile bounds
+        """
+        if len(self.merged_bounds) != len(merged_bounds):
+            raise ValueError(f"Size of merged bounds: must match the number of total channel: {len(self.merged_bounds)}")
         self.merged_bounds = merged_bounds
 
     def optimize(self, mask):
@@ -99,7 +123,6 @@ def merge_coils(self, unshimmed, affine):
         """
 
         coil_profiles_list = []
-        bounds = []
 
         # Define the nibabel unshimmed array
         nii_unshimmed = nib.Nifti1Image(unshimmed, affine)
@@ -109,15 +132,48 @@ def merge_coils(self, unshimmed, affine):
 
             # Resample a coil on the unshimmed image
             resampled_coil = resample_from_to(nii_coil, nii_unshimmed).get_fdata()
+            coil_profiles_list.append(resampled_coil)
 
+        coil_profiles = np.concatenate(coil_profiles_list, axis=3)
+
+        bounds = self.merge_bounds()
+
+        return coil_profiles, bounds
+
+    def merge_bounds(self):
+        """
+        Merge the coil profile bounds into a single array.
+
+        Returns:
+            list: list of bounds corresponding to each merged coils
+        """
+
+        bounds = []
+        for coil in self.coils:
             # Concat coils and bounds
-            coil_profiles_list.append(resampled_coil)
             for a_bound in coil.coef_channel_minmax:
                 bounds.append(a_bound)
 
-        coil_profiles = np.concatenate(coil_profiles_list, axis=3)
+        return bounds
 
-        return coil_profiles, bounds
+    def initial_guess_mean_bounds(self):
+        """
+        Calculates the initial guess from the bounds, sets it to the mean of the bounds
+
+        Returns:
+            np.ndarray: 1d array (n_channels) of coefficient representing the initial guess
+
+        """
+        current_0 = []
+        for bounds in self.merged_bounds:
+            avg = np.mean(bounds)
+
+            if np.isnan(avg):
+                current_0.append(0)
+            else:
+                current_0.append(avg)
+
+        return np.array(current_0)
 
     def _check_sizing(self, mask):
         """

shimmingtoolbox/optimizer/lsq_optimizer.py

@@ -8,6 +8,10 @@
 
 
 class LsqOptimizer(Optimizer):
+    """ Optimizer object that stores coil profiles and optimizes an unshimmed volume given a mask.
+        Use optimize(args) to optimize a given mask. The algorithm uses a least squares solver to find the best shim.
+        It supports bounds for each channel as well as a bound for the absolute sum of the channels.
+    """
 
     def _residuals(self, coef, unshimmed_vec, coil_mat):
         """
@@ -68,7 +72,7 @@ def _apply_sum_constraint(inputs, indexes, coef_sum_max):
             start_index = end_index
 
         # Set up output currents
-        currents_0 = np.zeros(n_channels)
+        currents_0 = self.initial_guess_mean_bounds()
 
         # Optimize
         currents_sp = opt.minimize(self._residuals, currents_0,