Add ability to create fieldmaps (#169)

* Initial set up for CLI and test

* Add st_prepare_fieldmap as an entry_point

* Added comments and started to implement

* Added read_nii and echo_time_diff

* Now checking the input phase is a phase diff

For the "single-echo" case, which corresponds to phase diff input

* Added TODOs, added call to prelude, added unwrapper case

* Fixed wrong input for mag

* Deal with 4d phase data

* Divide by echo_time, output fieldmap, add support for 2d fieldmaps

* start implementing dual echo

* Seperate single and dual echo test

* Update dual echo method

* Add magnitude and remove unecessary debug lines

* Move unwrapper to seperate API

* Add unwrap_phase

* Add tests for unwrap_phase and update unwrap_phase to correct found errors

* Cleanup

* Transfer to an API

* Rename dli test to follow convention

* Add tests for prepare_fieldmap API and fix errors found with tests

* Add new tests for error handling

* Add comments

* Updated help

* Added TODO

* Minor comments

* Added integrity test for prepare fieldmap

* More intuitive if, and add NotImplementedError test

* Add and fix comments

* Add not implemented error

* Switch phasediff variable to phase_img

* Refactor unwrap_phase

* Fix merge with master, make mag optional in unwrap phase and prepare fieldmap and fix comments

Co-authored-by: Julien Cohen-Adad <jcohen@polymtl.ca>

setup.py

@@ -21,6 +21,7 @@
         'console_scripts': [
             "st_download_data=shimmingtoolbox.cli.download_data:download_data",
             "st_dicom_to_nifti=shimmingtoolbox.cli.dicom_to_nifti:dicom_to_nifti_cli",
+            "st_prepare_fieldmap= shimmingtoolbox.cli.prepare_fieldmap:prepare_fieldmap_cli"
         ]
     },
     packages=find_packages(exclude=["docs"]),

shimmingtoolbox/cli/prepare_fieldmap.py

@@ -0,0 +1,76 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*
+
+import click
+import os
+import math
+import numpy as np
+import nibabel as nib
+
+from nibabel import load as load_nib
+
+from shimmingtoolbox.load_nifti import read_nii
+from shimmingtoolbox.prepare_fieldmap import prepare_fieldmap
+
+CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help'])
+
+
+@click.command(
+    context_settings=CONTEXT_SETTINGS,
+)
+@click.argument('phase', nargs=-1, type=click.Path(exists=True), required=True)
+@click.option('-mag', 'fname_mag', type=click.Path(exists=True), required=False, help="Input path of mag nifti file")
+@click.option('-unwrapper', type=click.Choice(['prelude']), default='prelude', help="Algorithm for unwrapping")
+@click.option('-output', 'fname_output', type=click.Path(), default=os.curdir, help="Output filename for the fieldmap")
+@click.option('-mask', 'fname_mask', type=click.Path(exists=True), help="Input path for a mask. Used for PRELUDE")
+@click.option('-threshold', 'threshold', type=float, help="Threshold for masking. Used for: PRELUDE")
+def prepare_fieldmap_cli(phase, fname_mag, unwrapper, fname_output, fname_mask, threshold):
+    """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.
+
+    phase: Input path of phase nifti file(s), in ascending order: echo1, echo2, etc.
+    """
+
+    # Import phase
+    list_phase = []
+    echo_times = []
+    for i_echo in range(len(phase)):
+        nii_phase, json_phase, phase_img = read_nii(phase[i_echo], auto_scale=True)
+        # Add pi since read_nii returns phase between 0 and 2pi whereas prepare_fieldmap accepts between -pi to pi
+        phase_img -= math.pi
+
+        list_phase.append(phase_img)
+        # Special case for echo_times if input is a phasediff
+        if len(phase) == 1:
+            # Check that the input phase is indeed a phasediff, by checking the existence of two echo times in the
+            # metadata
+            if not ('EchoTime1' in json_phase) or not ('EchoTime2' in json_phase):
+                raise RuntimeError(
+                    "The JSON file of the input phase should include the fields EchoTime1 and EchoTime2 if"
+                    "it is a phase difference.")
+            echo_times = [json_phase['EchoTime1'], json_phase['EchoTime2']]  # [s]
+        else:
+            echo_times.append(json_phase['EchoTime'])
+
+    # Get affine from nii
+    affine = nii_phase.affine
+
+    # If fname_mag is not an input define mag as None
+    if fname_mag is not None:
+        mag = load_nib(fname_mag).get_fdata()
+    else:
+        mag = None
+
+    # Import mask
+    if fname_mask is not None:
+        mask = nib.load(fname_mask).get_fdata()
+    else:
+        mask = None
+
+    fieldmap_hz = prepare_fieldmap(list_phase, echo_times, affine, mag=mag, unwrapper=unwrapper, mask=mask,
+                                   threshold=threshold)
+
+    # Save NIFTI
+    nii_fieldmap = nib.Nifti1Image(fieldmap_hz, affine)
+    nib.save(nii_fieldmap, fname_output)

shimmingtoolbox/prepare_fieldmap.py

@@ -0,0 +1,84 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*
+
+import math
+import numpy as np
+
+from shimmingtoolbox.unwrap.unwrap_phase import unwrap_phase
+
+
+def prepare_fieldmap(phase, echo_times, affine, unwrapper='prelude', mag=None, mask=None, threshold=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.
+
+    Args:
+        phase (list): List of phase values in a numpy.ndarray. The numpy array can be [x, y], [x, y, z] or [x, y, z, t].
+                      The values must range from [-pi to pi].
+        echo_times (list): List of echo times in seconds for each echo. The number of echotimes must match the number of
+                           echoes. It input is a phasediff (1 phase), input 2 echotimes.
+        affine (numpy.ndarray): 4x4 affine matrix.
+        unwrapper (str): Unwrapper to use for phase unwrapping. Supported: prelude.
+        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].
+        threshold: Prelude parameter used for masking.
+
+    Returns
+        numpy.ndarray: Unwrapped fieldmap in Hz.
+    """
+    # Check inputs
+    for i_echo in range(len(phase)):
+        # Check that the output phase is in radian (Note: the test below is not 100% bullet proof)
+        if (phase[i_echo].max() > math.pi) or (phase[i_echo].min() < -math.pi):
+            raise RuntimeError("read_nii must range from -pi to pi.")
+
+    # Check that the input echotimes are the appropriate size by looking at phase
+    is_phasediff = (len(phase) == 1 and len(echo_times) == 2)
+    if not is_phasediff:
+        if len(phase) != len(echo_times) or (len(phase) == 1 and len(echo_times) == 1):
+            raise RuntimeError("Phasediff must have 2 echotime points. Otherwise the number of echoes must match the"
+                               " number of echo times.")
+
+    # Make sure mag is the right shape
+    if mag is not None:
+        if mag.shape != phase[0].shape:
+            raise RuntimeError("mag and phase must have the same dimensions.")
+
+    # Make sure mask has the right shape
+    if mask is not None:
+        if mask.shape != phase[0].shape:
+            raise RuntimeError("Shape of mask and phase must match.")
+
+    # Get the time between echoes and calculate phase difference depending on number of echoes
+    if len(phase) == 1:
+        # phase should be a phasediff
+        phasediff = phase[0]
+        echo_time_diff = echo_times[1] - echo_times[0]  # [s]
+
+    elif len(phase) == 2:
+        echo_0 = phase[0]
+        echo_1 = phase[1]
+
+        # Calculate phasediff using complex difference
+        comp_0 = np.ones_like(echo_0) * np.exp(-1j * echo_0)
+        comp_1 = np.ones_like(echo_1) * np.exp(1j * echo_1)
+        phasediff = np.angle(comp_0 * comp_1)
+
+        # Calculate the echo time difference
+        echo_time_diff = echo_times[1] - echo_times[0]  # [s]
+
+    else:
+        # TODO: More echoes
+        # TODO: Add method once multiple methods are implemented
+        raise NotImplementedError(f"This number of phase input is not supported: {len(phase)}.")
+
+    # Run the unwrapper
+    phasediff_unwrapped = unwrap_phase(phasediff, affine, unwrapper=unwrapper, mag=mag, mask=mask, threshold=threshold)
+
+    # TODO: correct for potential wraps between time points
+
+    # Divide by echo time
+    fieldmap_rad = phasediff_unwrapped / echo_time_diff  # [rad / s]
+    fieldmap_hz = fieldmap_rad / (2 * math.pi)  # [Hz]
+
+    return fieldmap_hz

shimmingtoolbox/unwrap/unwrap_phase.py

@@ -0,0 +1,72 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+""" Wrapper to different unwrapping algorithms. """
+
+import numpy as np
+from shimmingtoolbox.unwrap.prelude import prelude
+
+
+def unwrap_phase(phase, affine, unwrapper='prelude', mag=None, mask=None, threshold=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.
+
+    Args:
+        phase (numpy.ndarray): 2D, 3D or 4D radian values [-pi to pi] to perform phase unwrapping.
+                               Supported shapes: [x, y], [x, y, z] or [x, y, z, t].
+        affine (numpy.ndarray): 2D array (4x4) containing the transformation coefficients. Can be acquired by :
+            nii = nib.load("nii_path")
+            affine = nii.affine
+        unwrapper (str, optional): Unwrapper algorithm name. Possible values: ``prelude``.
+        mag (numpy.ndarray): 2D, 3D or 4D magnitude data corresponding to phase data. Shape must be the same as
+                     ``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.
+
+    Returns:
+        numpy.ndarray: Unwrapped phase image.
+    """
+
+    if unwrapper == 'prelude':
+        mag4d = None
+        mask4d = None
+        if phase.ndim == 2:
+            phase4d = np.expand_dims(phase, -1)
+            if mag is not None:
+                mag4d = np.expand_dims(mag, -1)
+            if mask is not None:
+                mask4d = np.expand_dims(mask, -1)
+
+            mag = mag4d
+            mask = mask4d
+
+            phase3d_unwrapped = prelude(phase4d, affine, mag=mag, mask=mask, threshold=threshold)
+
+            phase_unwrapped = phase3d_unwrapped[..., 0]
+
+        elif phase.ndim == 3:
+            phase_unwrapped = prelude(phase, affine, mag=mag, mask=mask, threshold=threshold)
+
+        elif phase.ndim == 4:
+            phase_unwrapped = np.zeros_like(phase)
+            for i_t in range(phase.shape[3]):
+                mask3d = None
+                mag3d = None
+
+                phase3d = phase[..., i_t]
+                if mag is not None:
+                    mag3d = mag[..., i_t]
+                if mask is not None:
+                    mask3d = mask[..., i_t]
+
+                mask_input = mask3d
+                mag_input = mag3d
+
+                phase_unwrapped[..., i_t] = prelude(phase3d, affine, mag=mag_input, mask=mask_input, threshold=threshold)
+
+        else:
+            raise RuntimeError("Shape of input phase is not supported.")
+
+    else:
+        raise NotImplementedError(f'The unwrap function {unwrapper} is not implemented.')
+
+    return phase_unwrapped

test/cli/test_cli_prepare_fieldmap.py

@@ -0,0 +1,47 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*
+
+from click.testing import CliRunner
+import os
+import pathlib
+import tempfile
+import pytest
+
+from shimmingtoolbox.cli.prepare_fieldmap import prepare_fieldmap_cli
+from shimmingtoolbox import __dir_testing__
+
+
+@pytest.mark.prelude
+def test_cli_prepare_fieldmap_1_echo():
+    with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp:
+        runner = CliRunner()
+
+        fname_phasediff = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap',
+                                       'sub-example_phasediff.nii.gz')
+        fname_mag = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap',
+                                 'sub-example_magnitude1.nii.gz')
+        fname_output = os.path.join(tmp, 'fieldmap.nii.gz')
+
+        result = runner.invoke(prepare_fieldmap_cli, [fname_phasediff, '-mag', fname_mag, '-output', fname_output],
+                               catch_exceptions=False)
+
+        assert result.exit_code == 0
+        assert os.path.isfile(fname_output)
+
+
+@pytest.mark.prelude
+def test_cli_prepare_fieldmap_2_echos():
+    with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp:
+        runner = CliRunner()
+
+        fname_phase1 = os.path.join(__dir_testing__, 'sub-fieldmap', 'fmap', 'sub-fieldmap_phase1.nii.gz')
+        fname_phase2 = os.path.join(__dir_testing__, 'sub-fieldmap', 'fmap', 'sub-fieldmap_phase2.nii.gz')
+        fname_output = os.path.join(tmp, 'fieldmap.nii.gz')
+
+        fname_mag = os.path.join(__dir_testing__, 'sub-fieldmap', 'fmap', 'sub-fieldmap_magnitude1.nii.gz')
+
+        result = runner.invoke(prepare_fieldmap_cli, [fname_phase1, fname_phase2, '-mag', fname_mag,
+                                                      '-output', fname_output], catch_exceptions=False)
+
+        assert result.exit_code == 0
+        assert os.path.isfile(fname_output)