Merge 45f8b91 into 4f4869d

setup.py

@@ -55,6 +55,7 @@
         "dataclasses",
         "raven",
         "joblib",
+        "quadprog",
     ],
     extras_require={
         'docs': ["sphinx>=1.7", "sphinx_rtd_theme>=1.2.2", "sphinx-click"],

shimmingtoolbox/cli/b0shim.py

@@ -67,10 +67,10 @@ def b0shim_cli():
                    "'--slice-factor' value is '3', then with the 'sequential' mode, shimming will be performed "
                    "independently on the following groups: {0,1,2}, {3,4,5}, etc. With the mode 'interleaved', "
                    "it will be: {0,2,4}, {1,3,5}, etc.")
-@click.option('--optimizer-method', 'method', type=click.Choice(['least_squares', 'pseudo_inverse',
-                                                                 ]), required=False,
-              default='least_squares', show_default=True,
-              help="Method used by the optimizer. LS will respect the constraints, PS will not respect the constraints")
+@click.option('--optimizer-method', 'method', type=click.Choice(['least_squares', 'pseudo_inverse', 'quad_prog']),
+              required=False, default='quad_prog', show_default=True,
+              help="Method used by the optimizer. LS, and QP will respect the constraints,"
+                  "PS will not respect the constraints")
 @click.option('--regularization-factor', 'reg_factor', type=click.FLOAT, required=False, default=0.0, show_default=True,
               help="Regularization factor for the current when optimizing. A higher coefficient will penalize higher "
                    "current values while 0 provides no regularization. Not relevant for 'pseudo-inverse' "
@@ -505,9 +505,11 @@ def _save_to_text_file_static(coil, coefs, list_slices, path_output, o_format, o
                    "'--slice-factor' value is '3', then with the 'sequential' mode, shimming will be performed "
                    "independently on the following groups: {0,1,2}, {3,4,5}, etc. With the mode 'interleaved', "
                    "it will be: {0,2,4}, {1,3,5}, etc.")
-@click.option('--optimizer-method', 'method', type=click.Choice(['least_squares', 'pseudo_inverse']), required=False,
-              default='least_squares', show_default=True,
-              help="Method used by the optimizer. LS will respect the constraints, PS will not respect the constraints")
+@click.option('--optimizer-method', 'method', type=click.Choice(['least_squares', 'pseudo_inverse',
+                                                                 'quad_prog']), required=False,
+              default='quad_prog', show_default=True,
+              help="Method used by the optimizer. LS and QP will respect the constraints,"
+                   "PS will not respect the constraints")
 @click.option('--optimizer-criteria', 'opt_criteria', type=click.Choice(['mse', 'mae']), required=False,
               default='mse', show_default=True,
               help="Criteria of optimization for the optimizer 'least_squares'."

shimmingtoolbox/optimizer/basic_optimizer.py

@@ -22,9 +22,9 @@ class Optimizer(object):
 
     Attributes:
         coils (ListCoil): List of Coil objects containing the coil profiles and related constraints
-        unshimmed (numpy.ndarray): 3d array of unshimmed volume
-        unshimmed_affine (numpy.ndarray): 4x4 array containing the qform affine transformation for the unshimmed array
-        merged_coils (numpy.ndarray): 4d array containing all coil profiles resampled onto the target unshimmed array
+        unshimmed (np.ndarray): 3d array of unshimmed volume
+        unshimmed_affine (np.ndarray): 4x4 array containing the qform affine transformation for the unshimmed array
+        merged_coils (np.ndarray): 4d array containing all coil profiles resampled onto the target unshimmed array
                                       concatenated on the 4th dimension. See self.merge_coils() for more details
         merged_bounds (list): list of bounds corresponding to each merged coils: merged_bounds[3] is the (min, max)
                               bound for merged_coils[..., 3]
@@ -36,8 +36,8 @@ def __init__(self, coils: ListCoil, unshimmed, affine):
 
         Args:
             coils (ListCoil): List of Coil objects containing the coil profiles and related constraints
-            unshimmed (numpy.ndarray): 3d array of unshimmed volume
-            affine (numpy.ndarray): 4x4 array containing the affine transformation for the unshimmed array
+            unshimmed (np.ndarray): 3d array of unshimmed volume
+            affine (np.ndarray): 4x4 array containing the affine transformation for the unshimmed array
         """
         # Logging
         self.logger = logging.getLogger()
@@ -55,8 +55,8 @@ 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
+            unshimmed (np.ndarray): 3d array of unshimmed volume
+            affine: (np.ndarray): 4x4 array containing the qform affine transformation for the unshimmed array
         """
         # Check dimensions of unshimmed map
         if unshimmed.ndim != 3:
@@ -89,38 +89,57 @@ def optimize(self, mask):
         Optimize unshimmed volume by varying current to each channel
 
         Args:
-            mask (numpy.ndarray): 3d array of integers marking volume for optimization. Must be the same shape as
+            mask (np.ndarray): 3d array of integers marking volume for optimization. Must be the same shape as
                                   unshimmed
 
         Returns:
-            numpy.ndarray: Coefficients corresponding to the coil profiles that minimize the objective function.
+            np.ndarray: Coefficients corresponding to the coil profiles that minimize the objective function.
                            The shape of the array returned has shape corresponding to the total number of channels
         """
+        coil_mat, unshimmed_vec = self.get_coil_mat_and_unshimmed(mask)
+
+        output = -1 * scipy.linalg.pinv(coil_mat) @ unshimmed_vec  # N x mV' @ mV'
+
+        return output
+
+    def get_coil_mat_and_unshimmed(self, mask):
+        """
+        Returns the coil matrix, and the unshimmed vector used for the optimization
+
+        Args:
+            mask (np.ndarray): 3d array of integers marking volume for optimization. Must be the same shape as
+                              unshimmed
+        Returns:
+            (tuple) : tuple containing:
+                * np.ndarray: 2D flattened array (point, channel) of masked coils
+                              (axis 0 must align with unshimmed_vec)
+                * np.ndarray: 1D flattened array (point) of the masked unshimmed map
+
+        """
         # Check for sizing errors
         self._check_sizing(mask)
-
-        # Optimize
+        # Define coil profiles
+        n_channels = self.merged_coils.shape[3]
         mask_vec = mask.reshape((-1,))
+        # Reshape coil profile: X, Y, Z, N --> N, X, Y, Z --> N, [mask.shape]
+        # --> N, mask.size --> mask.size, N --> masked points, N
+        merged_coils_reshaped = np.reshape(np.transpose(self.merged_coils, axes=(3, 0, 1, 2)),
+                                           (n_channels, -1))
+        masked_points_indices = np.where(mask_vec != 0)
 
-        # Simple pseudo-inverse optimization
-        # Reshape coil profile: X, Y, Z, N --> [mask.shape], N
-        #   --> N, [mask.shape] --> N, mask.size --> mask.size, N --> masked points, N
-        coil_mat = np.reshape(np.transpose(self.merged_coils, axes=(3, 0, 1, 2)),
-                              (self.merged_coils.shape[3], -1)).T[mask_vec != 0, :]  # masked points x N
-        unshimmed_vec = np.reshape(self.unshimmed, (-1,))[mask_vec != 0]  # mV'
+        coil_mat = merged_coils_reshaped[:, masked_points_indices[0]].T  # masked points x N
+        unshimmed_vec = np.reshape(self.unshimmed, (-1,))[masked_points_indices[0]]  # mV'
 
-        output = -1 * scipy.linalg.pinv(coil_mat) @ unshimmed_vec  # N x mV' @ mV'
-
-        return output
+        return coil_mat, unshimmed_vec
 
     def merge_coils(self, unshimmed, affine):
         """
         Uses the list of coil profiles to return a resampled concatenated list of coil profiles matching the
         unshimmed image. Bounds are also concatenated and returned.
 
         Args:
-            unshimmed (numpy.ndarray): 3d array of unshimmed volume
-            affine (numpy.ndarray): 4x4 array containing the affine transformation for the unshimmed array
+            unshimmed (np.ndarray): 3d array of unshimmed volume
+            affine (np.ndarray): 4x4 array containing the affine transformation for the unshimmed array
         """
 
         coil_profiles_list = []
@@ -168,7 +187,7 @@ def _check_sizing(self, mask):
         Helper function to check array sizing
 
         Args:
-            mask (numpy.ndarray): 3d array of integers marking volume for optimization. Must be the same shape as
+            mask (np.ndarray): 3d array of integers marking volume for optimization. Must be the same shape as
                                   unshimmed
         """