Bug fix related to updating FSLeyes version and Shimming Toolbox vers…

…ion (#483)

* Allow 1D list

* Allow masks with from different geometries

* Add investigational device 7T

* Fix typo

* Add fmap and anat from B.U.F.F in dcm2bids config file

* Allow masks from other geometries

* Bring 2pi offset as close to the mean as possible for phase difference fieldmaps

* Better visuals

* Make the figure bigger

* Refactor bids config file

* Fix config file to use re

examples/demo_realtime_shimming.sh

@@ -1,8 +1,8 @@
 #!/usr/bin/env bash
 #
-# This function will generate static and dynamic (due to respiration) Gx, Gy, Gz components based on a fieldmap time
+# This function will generate static and dynamic (due to respiration) Gx, Gy, Gz components based on a field map time
 # series (magnitude and phase images) and respiratory trace information obtained from Siemens bellows. An additional
-# multi-gradient echo (MGRE) magnitude image is used to generate an ROI and resample the static and dynamic Gx, Gy, Gz
+# multi-gradient echo (MGRE) magnitude image is used to generate a ROI and resample the static and dynamic Gx, Gy, Gz
 # component maps to match the MGRE image. Lastly the average Gx, Gy, Gz values within the ROI are computed for each
 # slice.
 

examples/general_demo.py

@@ -54,13 +54,13 @@ def general_demo(path_output=os.path.join(os.path.curdir, 'output_dir')):
     dicom_to_nifti(path_dicom_unsorted, path_nifti, subject_id='sub-01')
 
     # Open phase data
-    fname_phases = glob.glob(os.path.join(path_nifti, 'sub-01', 'fmap', '*phase*.nii.gz'))
+    fname_phases = glob.glob(os.path.join(path_nifti, 'sub-01', 'anat', '*phase*.nii.gz'))
 
     nii_phase_e1, _, _ = read_nii(fname_phases[0])
     nii_phase_e2, _, _ = read_nii(fname_phases[1])
 
     # Open mag data
-    fname_mags = glob.glob(os.path.join(path_nifti, 'sub-01', 'fmap', '*magnitude*.nii.gz'))
+    fname_mags = glob.glob(os.path.join(path_nifti, 'sub-01', 'anat', '*magnitude*.nii.gz'))
 
     nii_mag_e1 = nib.load(fname_mags[0])
     nii_mag_e2 = nib.load(fname_mags[1])

shimmingtoolbox/cli/prepare_fieldmap.py

@@ -33,7 +33,7 @@
                    "standard data, it would be preferable to set this option to False and input your phase data from "
                    "-pi to pi to avoid unwanted rescaling")
 @click.option('--mask', 'fname_mask', type=click.Path(exists=True),
-              help="Input path for a mask. Mask must be the same shape as the array of each PHASE input.")
+              help="Input path for a mask.")
 @click.option('--threshold', 'threshold', type=float, show_default=True, default=0.05,
               help="Threshold for masking if no mask is provided. Allowed range: [0, 1] where all scaled values lower "
                    "than the threshold are set to 0.")
@@ -74,8 +74,7 @@ def prepare_fieldmap_uncli(phase, fname_mag, unwrapper='prelude',
         autoscale (bool): Tells whether to auto rescale phase inputs according to manufacturer standards. If you have
                           non siemens data not automatically converted from dcm2niix, you should set this to False and
                           input phase data from -pi to pi.
-        fname_mask (str): Input path for a mask. Mask must be the same shape as the array of each PHASE input.
-                          Used for PRELUDE
+        fname_mask (str): Input path for a mask. Used for PRELUDE
         threshold (float): Threshold for masking. Used for: PRELUDE
         fname_save_mask (str): Filename of the mask calculated by the unwrapper
         gaussian_filter (bool): Gaussian filter for B0 map
@@ -106,8 +105,7 @@ def prepare_fieldmap_cli_inputs(phase, fname_mag, unwrapper, autoscale, fname_ma
         autoscale (bool): Tells whether to auto rescale phase inputs according to manufacturer standards. If you have
                           non siemens data not automatically converted from dcm2niix, you should set this to False and
                           input phase data from -pi to pi.
-        fname_mask (str): Input path for a mask. Mask must be the same shape as the array of each PHASE input.
-                          Used for PRELUDE
+        fname_mask (str): Input path for a mask. Used for PRELUDE
         threshold (float): Threshold for masking. Used for: PRELUDE
         gaussian_filter (bool): Gaussian filter for B0 map
         sigma (float): Standard deviation of gaussian filter. Used for: gaussian_filter
@@ -152,12 +150,12 @@ def prepare_fieldmap_cli_inputs(phase, fname_mag, unwrapper, autoscale, fname_ma
 
     # Import mask
     if fname_mask is not None:
-        mask = nib.load(fname_mask).get_fdata()
+        nii_mask = nib.load(fname_mask)
     else:
-        mask = None
+        nii_mask = None
 
     fieldmap_hz, save_mask = prepare_fieldmap(list_nii_phase, echo_times, mag=mag, unwrapper=unwrapper,
-                                              mask=mask, threshold=threshold, gaussian_filter=gaussian_filter,
+                                              nii_mask=nii_mask, threshold=threshold, gaussian_filter=gaussian_filter,
                                               sigma=sigma, fname_save_mask=fname_save_mask)
 
     # Create nii fieldmap

shimmingtoolbox/coils/coil.py

@@ -196,9 +196,28 @@ def get_scanner_constraints(manufacturers_model_name, order=2):
                                                        [-3551.29, 3551.29],
                                                        [-3487.302, 3487.302]])
 
+    elif manufacturers_model_name == "Investigational_Device_7T":
+        constraints = {
+            "name": "Investigational_Device_7T",
+            "coef_channel_minmax": [],
+            "coef_sum_max": None
+        }
+        if order >= 0:
+            pass
+            # todo: f0 min and max is wrong
+        constraints["coef_channel_minmax"].append([None, None])
+        if order >= 1:
+            for _ in range(3):
+                constraints["coef_channel_minmax"].append([-5000, 5000])
+        if order >= 2:
+            constraints["coef_channel_minmax"].extend([[-1839.63, 1839.63],
+                                                       [-791.84, 791.84],
+                                                       [-791.84, 791.84],
+                                                       [-615.87, 615.87],
+                                                       [-615.87, 615.87]])
     else:
         logger.warning(f"Scanner: {manufacturers_model_name} constraints not yet implemented, constraints might not be "
-                       f"respected.")
+                       "respected.")
         constraints = {
             "name": "Unknown",
             "coef_sum_max": None

shimmingtoolbox/prepare_fieldmap.py

@@ -9,16 +9,17 @@
 from skimage.filters import gaussian
 from sklearn.linear_model import LinearRegression
 
-from shimmingtoolbox.unwrap.unwrap_phase import unwrap_phase
+from shimmingtoolbox.coils.coordinates import resample_from_to
 from shimmingtoolbox.masking.threshold import threshold as mask_threshold
+from shimmingtoolbox.unwrap.unwrap_phase import unwrap_phase
 from shimmingtoolbox.utils import fill
 
 logger = logging.getLogger(__name__)
 # Threshold to apply when creating a mask to remove 2*pi offsets
 VALIDITY_THRESHOLD = 0.2
 
 
-def prepare_fieldmap(list_nii_phase, echo_times, mag, unwrapper='prelude', mask=None, threshold=0.05,
+def prepare_fieldmap(list_nii_phase, echo_times, mag, unwrapper='prelude', nii_mask=None, threshold=0.05,
                      gaussian_filter=False, sigma=1, fname_save_mask=None):
     """ Creates fieldmap (in Hz) from phase images. This function accommodates multiple echoes (2 or more) and phase
     difference. This function also accommodates 4D phase inputs, where the 4th dimension represents the time, in case
@@ -31,7 +32,7 @@ def prepare_fieldmap(list_nii_phase, echo_times, mag, unwrapper='prelude', mask=
                            echoes. It input is a phasediff (1 phase), input 2 echotimes.
         unwrapper (str): Unwrapper to use for phase unwrapping. Supported: ``prelude``, ``skimage``.
         mag (numpy.ndarray): Array containing magnitude data relevant for ``phase`` input. Shape must match phase[echo].
-        mask (numpy.ndarray): Mask for masking output fieldmap. Must match shape of phase[echo].
+        nii_mask (nib.Nifti1Image): Mask for masking output fieldmap.
         threshold: Threshold for masking if no mask is provided. Allowed range: [0, 1] where all scaled values lower
                    than the threshold are set to 0.
         gaussian_filter (bool): Option of using a Gaussian filter to smooth the fieldmaps (boolean)
@@ -73,12 +74,36 @@ def prepare_fieldmap(list_nii_phase, echo_times, mag, unwrapper='prelude', mask=
         raise ValueError(f"Threshold should range from 0 to 1. Input value was: {threshold}")
 
     # Make sure mask has the right shape
-    if mask is None:
+    if nii_mask is None:
         # Define the mask using the threshold
         mask = mask_threshold(mag, threshold, scaled_thr=True)
     else:
-        if mask.shape != phase[0].shape:
-            raise ValueError("Shape of mask and phase must match.")
+        # Check that the mask is the right shape
+        if not np.all(nii_mask.shape == list_nii_phase[0].shape) or not np.all(
+                nii_mask.affine == list_nii_phase[0].affine):
+            logger.debug("Resampling mask on the target anat")
+            if list_nii_phase[0].ndim == 4:
+                nii_tmp_target = nib.Nifti1Image(list_nii_phase[0].get_fdata()[..., 0], list_nii_phase[0].affine,
+                                                 header=list_nii_phase[0].header)
+                if nii_mask.ndim == 3:
+                    tmp_mask = np.repeat(nii_mask.get_fdata()[..., np.newaxis], list_nii_phase[0].shape[-1], axis=-1)
+                    nii_tmp_mask = nib.Nifti1Image(tmp_mask, nii_mask.affine, header=nii_mask.header)
+                elif nii_mask.ndim == 4:
+                    nii_tmp_mask = nii_mask
+                else:
+                    raise ValueError("Mask must be 3D or 4D")
+            else:
+                # If it's not in 4d, assume it's a 3d mask
+                nii_tmp_target = list_nii_phase[0]
+                nii_tmp_mask = nii_mask
+
+            nii_mask_soft = resample_from_to(nii_tmp_mask, nii_tmp_target, order=1, mode='grid-constant')
+            tmp_mask = nii_mask_soft.get_fdata()
+            # Change soft mask into binary mask
+            mask = mask_threshold(tmp_mask, thr=0.001, scaled_thr=True)
+        else:
+            mask = nii_mask.get_fdata()
+
 
         logger.debug("A mask was provided, ignoring threshold value")
 
@@ -90,9 +115,9 @@ def prepare_fieldmap(list_nii_phase, echo_times, mag, unwrapper='prelude', mask=
         # Run the unwrapper
         phasediff_unwrapped = unwrap_phase(nii_phasediff, unwrapper=unwrapper, mag=mag, mask=mask,
                                            fname_save_mask=fname_save_mask)
-        # If it's 4d (i.e. there are timepoints)
-        if len(phasediff_unwrapped.shape) == 4:
-            phasediff_unwrapped = correct_2pi_offset(phasediff_unwrapped, mag, mask, VALIDITY_THRESHOLD)
+
+        # If 4d, correct 2pi offset between timepoints, if it's 3d, bring offset closest to the mean
+        phasediff_unwrapped = correct_2pi_offset(phasediff_unwrapped, mag, mask, VALIDITY_THRESHOLD)
 
         # Divide by echo time
         fieldmap_rad = phasediff_unwrapped / echo_time_diff  # [rad / s]
@@ -191,9 +216,10 @@ def correct_2pi_offset(unwrapped, mag, mask, validity_threshold):
         'correct' offset is assumed to be at time 0.
 
     Args:
-        unwrapped (numpy.ndarray): 4d array of the spatially unwrapped phase
-        mag (numpy.ndarray): 4d array containing the magnitude values of the phase
-        mask (numpy.ndarray): 4d mask of the unwrapped phase array
+        unwrapped (numpy.ndarray): Array of the spatially unwrapped phase. If there is a time dimension, the offset is
+                                   corrected in time, if unwrapped is 3D, the offset closest to 0 is chosen.
+        mag (numpy.ndarray): Array containing the magnitude values of the phase. Same shape as unwrapped.
+        mask (numpy.ndarray): Mask of the unwrapped phase array. Same shape as unwrapped.
         validity_threshold (float): Threshold to create a mask on each timepoints and assume as reliable phase data
 
     Returns:
@@ -208,29 +234,40 @@ def correct_2pi_offset(unwrapped, mag, mask, validity_threshold):
     # Logical and with the mask used for calculating the fieldmap
     validity_masks = np.logical_and(mask, validity_masks)
 
-    for i_time in range(1, unwrapped_cp.shape[3]):
-        # Take the region where both masks intersect
-        validity_mask = np.logical_and(validity_masks[..., i_time - 1], validity_masks[..., i_time])
+    if unwrapped.ndim == 4:
+        for i_time in range(1, unwrapped_cp.shape[3]):
+            # Take the region where both masks intersect
+            validity_mask = np.logical_and(validity_masks[..., i_time - 1], validity_masks[..., i_time])
 
-        # Calculate the means in the same validity mask
-        ma_0 = np.ma.array(unwrapped_cp[..., i_time - 1], mask=validity_mask == False)
-        mean_0 = np.ma.mean(ma_0)
-        ma_1 = np.ma.array(unwrapped_cp[..., i_time], mask=validity_mask == False)
-        mean_1 = np.ma.mean(ma_1)
+            # Calculate the means in the same validity mask
+            ma_0 = np.ma.array(unwrapped_cp[..., i_time - 1], mask=validity_mask == False)
+            mean_0 = np.ma.mean(ma_0)
+            ma_1 = np.ma.array(unwrapped_cp[..., i_time], mask=validity_mask == False)
+            mean_1 = np.ma.mean(ma_1)
 
-        # Calculate the number of offset by rounding to the nearest integer.
-        n_offsets_float = (mean_1 - mean_0) / (2 * np.pi)
-        n_offsets = round(n_offsets_float)
+            # Calculate the number of offset by rounding to the nearest integer.
+            n_offsets_float = (mean_1 - mean_0) / (2 * np.pi)
+            n_offsets = np.round(n_offsets_float)
 
-        if 0.3 < (n_offsets_float % 1) < 0.7:
-            logger.warning("The number of 2*pi offsets when calculating the fieldmap is close to "
-                           "ambiguous, verify the output fieldmap.")
+            if 0.3 < (n_offsets_float % 1) < 0.7:
+                logger.warning("The number of 2*pi offsets when calculating the fieldmap is close to "
+                               "ambiguous, verify the output fieldmap.")
 
+            if n_offsets != 0:
+                logger.info(f"Correcting for n*2pi phase offset, 'n' was: {n_offsets_float}")
+
+            logger.debug(f"Offset was: {n_offsets_float}")
+            # Remove n_offsets to unwrapped[..., i_time] only in the masked region
+            unwrapped_cp[..., i_time] -= mask[..., i_time] * n_offsets * (2 * np.pi)
+    else:
+        # unwrapped.ndim == 3:
+        ma_0 = np.ma.array(unwrapped_cp, mask=validity_masks == False)
+        mean_0 = np.ma.mean(ma_0)
+        n_offsets_float = mean_0 / (2 * np.pi)
+        n_offsets = np.round(n_offsets_float)
         if n_offsets != 0:
             logger.info(f"Correcting for n*2pi phase offset, 'n' was: {n_offsets_float}")
 
-        logger.debug(f"Offset was: {n_offsets_float}")
-        # Remove n_offsets to unwrapped[..., i_time] only in the masked region
-        unwrapped_cp[..., i_time] -= mask[..., i_time] * n_offsets * (2 * np.pi)
+        unwrapped_cp -= mask * n_offsets * (2 * np.pi)
 
     return unwrapped_cp

shimmingtoolbox/shim/sequencer.py

@@ -289,7 +289,7 @@ def get_mask(self, nii_mask_anat):
             nii_mask_anat_soft = resample_from_to(nii_mask_anat, self.nii_anat, order=1, mode='grid-constant')
             tmp_mask = nii_mask_anat_soft.get_fdata()
             # Change soft mask into binary mask
-            tmp_mask = threshold(tmp_mask, thr=0.001)
+            tmp_mask = threshold(tmp_mask, thr=0.001, scaled_thr=True)
             nii_mask_anat = nib.Nifti1Image(tmp_mask, nii_mask_anat_soft.affine, header=nii_mask_anat_soft.header)
             if logger.level <= getattr(logging, 'DEBUG') and self.path_output is not None:
                 nib.save(nii_mask_anat, os.path.join(self.path_output, "mask_static_resampled_on_anat.nii.gz"))
@@ -610,15 +610,15 @@ def plot_full_mask(self, unshimmed, shimmed_masked, mask):
         min_value = min(mt_unshimmed_masked.min(), mt_shimmed_masked.min())
         max_value = max(mt_unshimmed_masked.max(), mt_shimmed_masked.max())
 
-        fig = Figure(figsize=(8, 5))
+        fig = Figure(figsize=(15, 9))
         fig.suptitle(f"Fieldmaps\nFieldmap Coordinate System")
 
         ax = fig.add_subplot(1, 2, 1)
         ax.imshow(mt_unshimmed, cmap='gray')
         mt_unshimmed_masked[mt_unshimmed_masked == 0] = np.nan
         im = ax.imshow(mt_unshimmed_masked, vmin=min_value, vmax=max_value, cmap='viridis')
-        ax.set_title(f"Before shimming\nstd: {metric_unshimmed_std:.3}, mean: {metric_unshimmed_mean:.3}\n"
-                     f"mae: {metric_unshimmed_mae:.3}, rmse: {metric_unshimmed_rmse:.3}")
+        ax.set_title(f"Before shimming\nstd: {metric_unshimmed_std:.1f}, mean: {metric_unshimmed_mean:.1f}\n"
+                     f"mae: {metric_unshimmed_mae:.1f}, rmse: {metric_unshimmed_rmse:.1f}")
         ax.get_xaxis().set_visible(False)
         ax.get_yaxis().set_visible(False)
         divider = make_axes_locatable(ax)
@@ -629,8 +629,8 @@ def plot_full_mask(self, unshimmed, shimmed_masked, mask):
         ax.imshow(mt_unshimmed, cmap='gray')
         mt_shimmed_masked[mt_shimmed_masked == 0] = np.nan
         im = ax.imshow(mt_shimmed_masked, vmin=min_value, vmax=max_value, cmap='viridis')
-        ax.set_title(f"After shimming\nstd: {metric_shimmed_std:.3}, mean: {metric_shimmed_mean:.3}\n"
-                     f"mae: {metric_shimmed_mae:.3}, rmse: {metric_shimmed_rmse:.3}")
+        ax.set_title(f"After shimming\nstd: {metric_shimmed_std:.1f}, mean: {metric_shimmed_mean:.1f}\n"
+                     f"mae: {metric_shimmed_mae:.1f}, rmse: {metric_shimmed_rmse:.1f}")
         ax.get_xaxis().set_visible(False)
         ax.get_yaxis().set_visible(False)
         divider = make_axes_locatable(ax)

shimmingtoolbox/shim/shim_utils.py

@@ -160,7 +160,7 @@ def phys_to_shim_cs(coefs, manufacturer):
     """Convert a list of coefficients from RAS to the Shim Coordinate System
 
     Args:
-        coefs (np.ndarray): 1d list of coefficients in the physical RAS coordinate system of the manufacturer. The first
+        coefs (np.ndarray): Coefficients in the physical RAS coordinate system of the manufacturer. The first
                             dimension represents the different channels. (indexes 0, 1, 2 --> x, y, z...). If there are
                             more coefficients, they are of higher order and must correspond to the implementation of the
                             manufacturer. i.e. Siemens: *X, Y, Z, Z2, ZX, ZY, X2-Y2, XY*

test/test_dicom_to_nifti.py

@@ -34,7 +34,7 @@ def test_dicom_to_nifti():
             for modality in ['phase', 'magnitude']:
                 for ext in ['nii.gz', 'json']:
                     assert os.path.exists(os.path.join(path_nifti, f"sub-{subject_id}",
-                                                       'fmap', f"sub-{subject_id}_{modality}{i}.{ext}"))
+                                                       'anat', f"sub-{subject_id}_{modality}{i}.{ext}"))
 
 
 @pytest.mark.dcm2niix
@@ -67,7 +67,7 @@ def test_dicom_to_nifti_realtime_zshim(test_dcm2niix_installation):
             for ext in ['nii.gz', 'json']:
                 assert os.path.exists(
                     os.path.join(path_nifti, f"sub-{subject_id}", sequence_type,
-                                 f"sub-{subject_id}_unshimmed_e{i+1}.{ext}"))
+                                 f"sub-{subject_id}_magnitude{i+1}.{ext}"))
 
         sequence_type = 'func'
         for ext in ['nii.gz', 'json']: