Use --savemask option from prelude to calculate the output mask in st…

…_prepare_fieldmap (#393)

* Save mask calculated by prelude

* Remove import

shimmingtoolbox/cli/prepare_fieldmap.py

@@ -64,6 +64,11 @@ def prepare_fieldmap_cli(phase, fname_mag, unwrapper, fname_output, autoscale, f
     # Prepare the output
     create_output_dir(fname_output_v2, is_file=True)
 
+    # Save mask
+    if fname_save_mask is not None:
+        # If it is a path, add the default filename and create output directory
+        fname_save_mask = create_fname_from_path(fname_save_mask, MASK_OUTPUT_DEFAULT)
+
     # Import phase
     list_nii_phase = []
     echo_times = []
@@ -87,7 +92,7 @@ def prepare_fieldmap_cli(phase, fname_mag, unwrapper, fname_output, autoscale, f
     affine = nii_phase.affine
 
     # Magnitude image
-    mag = nib.load(fname_mag).get_fdata()
+    _, json_mag, mag = read_nii(fname_mag)
 
     # Import mask
     if fname_mask is not None:
@@ -97,13 +102,13 @@ def prepare_fieldmap_cli(phase, fname_mag, unwrapper, fname_output, autoscale, f
 
     fieldmap_hz, save_mask = prepare_fieldmap(list_nii_phase, echo_times, mag=mag, unwrapper=unwrapper,
                                               mask=mask, threshold=threshold, gaussian_filter=gaussian_filter,
-                                              sigma=sigma)
+                                              sigma=sigma, fname_save_mask=fname_save_mask)
 
     # Save fieldmap
     nii_fieldmap = nib.Nifti1Image(fieldmap_hz, affine, header=nii_phase.header)
     nib.save(nii_fieldmap, fname_output_v2)
 
-    # Save json
+    # Save fieldmap json
     json_fieldmap = json_phase
     if len(phase) > 1:
         for i_echo in range(len(echo_times)):
@@ -112,17 +117,10 @@ def prepare_fieldmap_cli(phase, fname_mag, unwrapper, fname_output, autoscale, f
     with open(fname_json, 'w') as outfile:
         json.dump(json_fieldmap, outfile, indent=2)
 
-    # Save mask
+    # save mask json
     if fname_save_mask is not None:
-        # If it is a path, add the default filename and create output directory
-        fname_save_mask = create_fname_from_path(fname_save_mask, MASK_OUTPUT_DEFAULT)
-        create_output_dir(fname_save_mask, is_file=True)
-
-        if fname_save_mask[-4:] != '.nii' and fname_save_mask[-7:] != '.nii.gz':
-            raise ValueError("Output filename must have one of the following extensions: '.nii', '.nii.gz'")
-
-        nii_fieldmap = nib.Nifti1Image(save_mask, affine, header=nii_phase.header)
-        nib.save(nii_fieldmap, fname_save_mask)
-        logger.info(f"Filename of the output mask is: {fname_save_mask}")
+        fname_mask_json = fname_save_mask.rsplit('.nii', 1)[0] + '.json'
+        with open(fname_mask_json, 'w') as outfile:
+            json.dump(json_mag, outfile, indent=2)
 
     logger.info(f"Filename of the fieldmap is: {fname_output_v2}")

shimmingtoolbox/prepare_fieldmap.py

@@ -16,7 +16,7 @@
 
 
 def prepare_fieldmap(list_nii_phase, echo_times, mag, unwrapper='prelude', mask=None, threshold=0.05,
-                     gaussian_filter=False, sigma=1):
+                     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
     multiple field maps are acquired across time for the purpose of real-time shimming experiments.
@@ -33,6 +33,8 @@ def prepare_fieldmap(list_nii_phase, echo_times, mag, unwrapper='prelude', mask=
                    than the threshold are set to 0.
         gaussian_filter (bool): Option of using a Gaussian filter to smooth the fieldmaps (boolean)
         sigma (float): Standard deviation of gaussian filter.
+        fname_save_mask (str): Filename of the mask calculated by the unwrapper
+
     Returns
         numpy.ndarray: Unwrapped fieldmap in Hz.
     """
@@ -102,7 +104,8 @@ def prepare_fieldmap(list_nii_phase, echo_times, mag, unwrapper='prelude', mask=
         raise NotImplementedError(f"This number of phase input is not supported: {len(phase)}.")
 
     # Run the unwrapper
-    phasediff_unwrapped = unwrap_phase(nii_phasediff, unwrapper=unwrapper, mag=mag, mask=mask)
+    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:

shimmingtoolbox/unwrap/prelude.py

@@ -17,7 +17,7 @@
 logger = logging.getLogger(__name__)
 
 
-def prelude(nii_wrapped_phase, mag=None, mask=None, threshold=None, is_unwrapping_in_2d=False):
+def prelude(nii_wrapped_phase, mag=None, mask=None, threshold=None, is_unwrapping_in_2d=False, fname_save_mask=None):
     """wrapper to FSL prelude
 
     This function enables phase unwrapping by calling FSL prelude on the command line. A mask can be provided to mask
@@ -31,6 +31,7 @@ def prelude(nii_wrapped_phase, mag=None, mask=None, threshold=None, is_unwrappin
                                         unwrapping
         threshold: Threshold value for automatic mask generation (Use either mask or threshold, not both)
         is_unwrapping_in_2d (bool, optional): prelude parameter to unwrap slice by slice
+        fname_save_mask (str): Filename of the mask calculated by the unwrapper
 
     Returns:
         numpy.ndarray: 3D array with the shape of `complex_array` of the unwrapped phase output from prelude
@@ -71,6 +72,13 @@ def prelude(nii_wrapped_phase, mag=None, mask=None, threshold=None, is_unwrappin
 
             nib.save(nii_mask, fname_mask)
 
+        # Save mask
+        if fname_save_mask is not None:
+            if fname_save_mask[-4:] != '.nii' and fname_save_mask[-7:] != '.nii.gz':
+                raise ValueError("Output filename must have one of the following extensions: '.nii', '.nii.gz'")
+
+            options.append(f'--savemask={fname_save_mask}')
+
         if threshold is not None:
             options += ['-t', str(threshold)]
             if mask is not None:

shimmingtoolbox/unwrap/unwrap_phase.py

@@ -11,7 +11,7 @@
 logger = logging.getLogger(__name__)
 
 
-def unwrap_phase(nii_phase_wrapped, unwrapper='prelude', mag=None, mask=None, threshold=None):
+def unwrap_phase(nii_phase_wrapped, unwrapper='prelude', mag=None, mask=None, threshold=None, fname_save_mask=None):
     """ Calls different unwrapping algorithms according to the specified `unwrapper` parameter. The function also
     allows to call the different unwrappers with more flexibility regarding input shape.
 
@@ -23,6 +23,7 @@ def unwrap_phase(nii_phase_wrapped, unwrapper='prelude', mag=None, mask=None, th
                      ``phase``.
         mask (numpy.ndarray): numpy array of booleans with shape of ``phase`` to mask during phase unwrapping.
         threshold (float): Prelude parameter, see prelude for more detail.
+        fname_save_mask (str): Filename of the mask calculated by the unwrapper
 
     Returns:
         numpy.ndarray: Unwrapped phase image.
@@ -45,19 +46,22 @@ def unwrap_phase(nii_phase_wrapped, unwrapper='prelude', mag=None, mask=None, th
             nii_2d = nib.Nifti1Image(phase2d, nii_phase_wrapped.affine, header=nii_phase_wrapped.header)
 
             logger.info(f"Unwrapping 1 volume")
-            phase3d_unwrapped = prelude(nii_2d, mag=mag, mask=mask, threshold=threshold)
+            phase3d_unwrapped = prelude(nii_2d, mag=mag, mask=mask, threshold=threshold,
+                                        fname_save_mask=fname_save_mask)
 
             phase_unwrapped = phase3d_unwrapped[..., 0]
 
         elif phase.ndim == 3:
             logger.info("Unwrapping 1 volume")
-            phase_unwrapped = prelude(nii_phase_wrapped, mag=mag, mask=mask, threshold=threshold)
+            phase_unwrapped = prelude(nii_phase_wrapped, mag=mag, mask=mask, threshold=threshold,
+                                      fname_save_mask=fname_save_mask)
 
         elif phase.ndim == 4:
 
             logger.info(f"Unwrapping {phase.shape[3]} volumes")
             phase_unwrapped = np.zeros_like(phase)
-            for i_t in range(phase.shape[3]):
+            n_t = phase.shape[3]
+            for i_t in range(n_t):
                 mask4d = None
                 mag4d = None
 
@@ -74,6 +78,11 @@ def unwrap_phase(nii_phase_wrapped, unwrapper='prelude', mag=None, mask=None, th
 
                 phase_unwrapped[..., i_t] = prelude(nii_4d, mag=mag_input, mask=mask_input, threshold=threshold)
 
+                # If it's the first volume, call it using save mask to save the mask only once
+                if i_t == 0:
+                    phase_unwrapped[..., i_t] = prelude(nii_4d, mag=mag_input, mask=mask_input, threshold=threshold,
+                                                        fname_save_mask=fname_save_mask)
+
         else:
             raise ValueError("Shape of input phase is not supported.")