Module ukat.mapping.t2

Functions

def three_param_eq(t, t2, m0, b)

Calculate the expected signal from the equation signal = M0 * exp(-t / T2) + b

Parameters

t : list
The times the signal will be calculated at
t2 : float
The T2 of the signal
m0 : float
The M0 of the signal
b : float
The baseline noise floor of the signal

Returns

signal : np.ndarray
The expected signal
def two_param_eq(t, t2, m0)

Calculate the expected signal from the equation signal = M0 * exp(-t / T2)

Parameters

t : list
The times the signal will be calculated at
t2 : float
The T2 of the signal
m0 : float
The M0 of the signal

Returns

signal : np.ndarray
The expected signal

Classes

class T2 (pixel_array, echo_list, affine, mask=None, noise_threshold=0, method='2p_exp', multithread='auto')

Attributes

t2_map : np.ndarray
The estimated T2 values in ms
t2_err : np.ndarray
The certainty in the fit of t2 in ms
m0_map : np.ndarray
The estimated M0 values
m0_err : np.ndarray
The certainty in the fit of m0
r2 : np.ndarray
The R-Squared value of the fit, values close to 1 indicate a good fit, lower values indicate a poorer fit
shape : tuple
The shape of the T2 map
n_te : int
The number of TE used to calculate the map
n_vox : int
The number of voxels in the map i.e. the product of all dimensions apart from TE

Initialise a T2 class instance.

Parameters

pixel_array : np.ndarray
An array containing the signal from each voxel at each echo time with the last dimension being time i.e. the array needed to generate a 3D T2 map would have dimensions [x, y, z, TE].
echo_list : list()
An array of the echo times used for the last dimension of the raw data. In milliseconds.
affine : np.ndarray
A matrix giving the relationship between voxel coordinates and world coordinates.
mask : np.ndarray, optional
A boolean mask of the voxels to fit. Should be the shape of the desired T2 map rather than the raw data i.e. omit the time dimension.
noise_threshold : float, optional
Default 0 Voxels with magnitude less than this threshold will not be used when fitting. This can be useful if the noise floor of the data is known.
method : {'2p_exp', '3p_exp'}, optional
Default 2p_exp The model the data is fit to. 2p_exp uses a two parameter exponential model (S = S0 * exp(-t / T2)) whereas 3p_exp uses a three parameter exponential model (S = S0 * exp(-t / T2) + b) to fit for noise/very long T2 components of the signal.
multithread : bool or 'auto', optional
Default 'auto'. If True, fitting will be distributed over all cores available on the node. If False, fitting will be carried out on a single thread. 'auto' attempts to apply multithreading where appropriate based on the number of voxels being fit.
Expand source code
class T2:
    """
    Attributes
    ----------
    t2_map : np.ndarray
        The estimated T2 values in ms
    t2_err : np.ndarray
        The certainty in the fit of `t2` in ms
    m0_map : np.ndarray
        The estimated M0 values
    m0_err : np.ndarray
        The certainty in the fit of `m0`
    r2 : np.ndarray
        The R-Squared value of the fit, values close to 1 indicate a good
        fit, lower values indicate a poorer fit
    shape : tuple
        The shape of the T2 map
    n_te : int
        The number of TE used to calculate the map
    n_vox : int
        The number of voxels in the map i.e. the product of all dimensions
        apart from TE
    """

    def __init__(self, pixel_array, echo_list, affine, mask=None,
                 noise_threshold=0, method='2p_exp', multithread='auto'):
        """Initialise a T2 class instance.

        Parameters
        ----------
        pixel_array : np.ndarray
            An array containing the signal from each voxel at each echo
            time with the last dimension being time i.e. the array needed to
            generate a 3D T2 map would have dimensions [x, y, z, TE].
        echo_list : list()
            An array of the echo times used for the last dimension of the
            raw data. In milliseconds.
        affine : np.ndarray
            A matrix giving the relationship between voxel coordinates and
            world coordinates.
        mask : np.ndarray, optional
            A boolean mask of the voxels to fit. Should be the shape of the
            desired T2 map rather than the raw data i.e. omit the time
            dimension.
        noise_threshold : float, optional
            Default 0
            Voxels with magnitude less than this threshold will not be used
            when fitting. This can be useful if the noise floor of the data
            is known.
        method : {'2p_exp', '3p_exp'}, optional
            Default `2p_exp`
            The model the data is fit to. 2p_exp uses a two parameter
            exponential model (S = S0 * exp(-t / T2)) whereas 3p_exp uses a
            three parameter exponential model (S = S0 * exp(-t / T2) + b) to
            fit for noise/very long T2 components of the signal.
        multithread : bool or 'auto', optional
            Default 'auto'.
            If True, fitting will be distributed over all cores available on
            the node. If False, fitting will be carried out on a single thread.
            'auto' attempts to apply multithreading where appropriate based
            on the number of voxels being fit.
        """
        # Some sanity checks
        assert (pixel_array.shape[-1]
                == len(echo_list)), 'Number of echoes does not match the ' \
                                    'number of time frames on the last axis ' \
                                    'of pixel_array'
        assert multithread is True \
               or multithread is False \
               or multithread == 'auto', f'multithreaded must be True,' \
                                         f'False or auto. You entered ' \
                                         f'{multithread}'

        if method != '2p_exp' and method != '3p_exp':
            raise ValueError(f'method can be 2p_exp or 3p_exp only. You '
                             f'specified {method}')

        # Normalise the data so its roughly in the same range across vendors
        self.scale = np.nanmax(pixel_array)
        self.pixel_array = pixel_array / self.scale

        self.shape = pixel_array.shape[:-1]
        self.n_te = pixel_array.shape[-1]
        self.n_vox = np.prod(self.shape)
        self.affine = affine
        # Generate a mask if there isn't one specified
        if mask is None:
            self.mask = np.ones(self.shape, dtype=bool)
        else:
            self.mask = mask.astype(bool)

        # Don't process any nan values
        self.mask[np.isnan(np.sum(pixel_array, axis=-1))] = False
        self.noise_threshold = noise_threshold / self.scale
        self.method = method
        self.echo_list = echo_list
        # Auto multithreading conditions
        if multithread == 'auto':
            if self.n_vox > 20:
                multithread = True
            else:
                multithread = False
        self.multithread = multithread

        # Fit data
        self.fitting_model = T2Model(self.pixel_array, self.echo_list,
                                self.method, self.mask, self.multithread)

        if self.noise_threshold > 0:
            self.fitting_model.threshold_noise(self.noise_threshold)
        popt, error, r2 = fitting.fit_image(self.fitting_model)
        self.t2_map = popt[0]
        self.m0_map = popt[1]
        self.t2_err = error[0]
        self.m0_err = error[1]
        self.r2 = r2

        if self.method == '3p_exp':
            self.b_map = popt[2]
            self.b_err = error[2]

        # Filter values that are very close to models upper bounds of T2 or
        # M0 out.
        threshold = 0.999  # 99.9% of the upper bound
        bounds_mask = ((self.t2_map > self.fitting_model.bounds[1][0] *
                        threshold) |
                       (self.m0_map > self.fitting_model.bounds[1][1] *
                        threshold))
        self.t2_map[bounds_mask] = 0
        self.m0_map[bounds_mask] = 0
        self.t2_err[bounds_mask] = 0
        self.m0_err[bounds_mask] = 0
        self.r2[bounds_mask] = 0
        if self.method == '3p_exp':
            self.b_map[bounds_mask] = 0
            self.b_err[bounds_mask] = 0

            self.b_map *= self.scale
            self.b_err *= self.scale

        # Scale the data back to the original scale
        self.m0_map *= self.scale
        self.m0_err *= self.scale
        self.noise_threshold *= self.scale

    def to_nifti(self, output_directory=os.getcwd(), base_file_name='Output',
                 maps='all'):
        """Exports some of the T2 class attributes to NIFTI.

        Parameters
        ----------
        output_directory : string, optional
            Path to the folder where the NIFTI files will be saved.
        base_file_name : string, optional
            Filename of the resulting NIFTI. This code appends the extension.
            Eg., base_file_name = 'Output' will result in 'Output.nii.gz'.
        maps : list or 'all', optional
            List of maps to save to NIFTI. This should either the string "all"
            or a list of maps from ["t2", "t2_err", "m0", "m0_err",
            "r2", "mask"].
        """
        os.makedirs(output_directory, exist_ok=True)
        base_path = os.path.join(output_directory, base_file_name)
        if maps == 'all' or maps == ['all']:
            maps = ['t2', 't2_err', 'm0', 'm0_err', 'r2', 'mask']
            if self.method == '3p_exp':
                maps.append('b')
                maps.append('b_err')
        if isinstance(maps, list):
            for result in maps:
                if result == 't2' or result == 't2_map':
                    t2_nifti = nib.Nifti1Image(self.t2_map, affine=self.affine)
                    nib.save(t2_nifti, base_path + '_t2_map.nii.gz')
                elif result == 't2_err':
                    t2_err_nifti = nib.Nifti1Image(self.t2_err,
                                                   affine=self.affine)
                    nib.save(t2_err_nifti, base_path + '_t2_err.nii.gz')
                elif result == 'm0' or result == 'm0_map':
                    m0_nifti = nib.Nifti1Image(self.m0_map, affine=self.affine)
                    nib.save(m0_nifti, base_path + '_m0_map.nii.gz')
                elif result == 'm0_err':
                    m0_err_nifti = nib.Nifti1Image(self.m0_err,
                                                   affine=self.affine)
                    nib.save(m0_err_nifti, base_path + '_m0_err.nii.gz')
                elif result == 'r2' or result == 'r2_map':
                    r2_nifti = nib.Nifti1Image(self.r2,
                                               affine=self.affine)
                    nib.save(r2_nifti, base_path + '_r2_map.nii.gz')
                elif result == 'mask':
                    mask_nifti = nib.Nifti1Image(self.mask.astype(np.uint16),
                                                 affine=self.affine)
                    nib.save(mask_nifti, base_path + '_mask.nii.gz')
                elif result == 'b' or result == 'b_map':
                    b_nifti = nib.Nifti1Image(self.b_map,
                                              affine=self.affine)
                    nib.save(b_nifti, base_path + '_b_map.nii.gz')
                elif result == 'b_err':
                    b_err_nifti = nib.Nifti1Image(self.b_err,
                                                  affine=self.affine)
                    nib.save(b_err_nifti, base_path + '_b_err.nii.gz')
        else:
            raise ValueError('No NIFTI file saved. The variable "maps" '
                             'should be "all" or a list of maps from '
                             '"["t2", "t2_err", "m0", "m0_err", "r2", '
                             '"mask"]".')
        return

    def get_fit_signal(self):
        """
        Get the fit signal from the model used to fit the data i.e. the
        simulated signal at each echo time given the estimated T2, M0
        (and baseline noise floor (b) if applicable).

        Returns
        -------
        fit_signal : np.ndarray
            An array containing the fit signal generated by the model
        """
        fit_signal = np.zeros((self.n_vox, self.n_te))
        if self.method == '2p_exp':
            params = np.array([self.t2_map.reshape(-1),
                               self.m0_map.reshape(-1)])
        elif self.method == '3p_exp':
            params = np.array([self.t2_map.reshape(-1),
                               self.m0_map.reshape(-1),
                               self.b_map.reshape(-1)])

        for n in range(self.n_vox):
            fit_signal[n] = self.fitting_model.t2_eq(self.echo_list,
                                                     *params[:, n])

        fit_signal = fit_signal.reshape((*self.shape, self.n_te))
        return fit_signal

Methods

def get_fit_signal(self)

Get the fit signal from the model used to fit the data i.e. the simulated signal at each echo time given the estimated T2, M0 (and baseline noise floor (b) if applicable).

Returns

fit_signal : np.ndarray
An array containing the fit signal generated by the model
def to_nifti(self, output_directory='/home/runner/work/ukat/ukat', base_file_name='Output', maps='all')

Exports some of the T2 class attributes to NIFTI.

Parameters

output_directory : string, optional
Path to the folder where the NIFTI files will be saved.
base_file_name : string, optional
Filename of the resulting NIFTI. This code appends the extension. Eg., base_file_name = 'Output' will result in 'Output.nii.gz'.
maps : list or 'all', optional
List of maps to save to NIFTI. This should either the string "all" or a list of maps from ["t2", "t2_err", "m0", "m0_err", "r2", "mask"].
class T2Model (pixel_array, te, method='2p_exp', mask=None, multithread=True)

A class containing the T2 fitting model

Parameters

pixel_array : np.ndarray
An array containing the signal from each voxel at each echo time with the last dimension being time i.e. the array needed to generate a 3D T2 map would have dimensions [x, y, z, TE].
te : np.ndarray
An array of the echo times used for the last dimension of the pixel_array. In milliseconds.
method : {'2p_exp', '3p_exp'}, optional
Default '2p_exp' The model the data is fit to. 2p_exp uses a two parameter exponential model (S = S0 * exp(-t / T2)) whereas 3p_exp uses a three parameter exponential model (S = S0 * exp(-t / T2) + b) to fit for noise/very long T2 components of the signal.
mask : np.ndarray, optional
A boolean mask of the voxels to fit. Should be the shape of the desired T2 map rather than the raw data i.e. omit the time dimension.
multithread : bool, optional
Default True If True, the fitting will be performed in parallel using all available cores
Expand source code
class T2Model(fitting.Model):
    def __init__(self, pixel_array, te, method='2p_exp', mask=None,
                 multithread=True):
        """
        A class containing the T2 fitting model

        Parameters
        ----------
        pixel_array : np.ndarray
            An array containing the signal from each voxel at each echo
            time with the last dimension being time i.e. the array needed to
            generate a 3D T2 map would have dimensions [x, y, z, TE].
        te : np.ndarray
            An array of the echo times used for the last dimension of the
            pixel_array. In milliseconds.
        method : {'2p_exp', '3p_exp'}, optional
            Default '2p_exp'
            The model the data is fit to. 2p_exp uses a two parameter
            exponential model (S = S0 * exp(-t / T2)) whereas 3p_exp uses a
            three parameter exponential model (S = S0 * exp(-t / T2) + b) to
            fit for noise/very long T2 components of the signal.
        mask : np.ndarray, optional
            A boolean mask of the voxels to fit. Should be the shape of the
            desired T2 map rather than the raw data i.e. omit the time
            dimension.
        multithread : bool, optional
            Default True
            If True, the fitting will be performed in parallel using all
            available cores
        """
        self.method = method

        if self.method == '2p_exp':
            self.t2_eq = two_param_eq
            super().__init__(pixel_array, te, self.t2_eq, mask, multithread)
            self.bounds = ([0, 0], [1000, 100])
            self.initial_guess = [20, 1]
        elif self.method == '3p_exp':
            self.t2_eq = three_param_eq
            super().__init__(pixel_array, te, self.t2_eq, mask,
                             multithread)
            self.bounds = ([0, 0, 0], [1000, 100, 1])
            self.initial_guess = [20, 1, 5e-4]

        self.generate_lists()

    def threshold_noise(self, threshold=0):
        """
        Remove voxel values below a certain threshold from the fitting
        process, useful if long echo times have been collected and thus
        thermal noise is being measured below a certain threshold rather
        than the T2 decay.

        Parameters
        ----------
        threshold : float, optional
            Default 0
            The threshold below which to remove values
        """
        for ind, (sig, te, p0) in enumerate(zip(self.signal_list,
                                                self.x_list,
                                                self.p0_list)):
            self.signal_list[ind] = np.array(
                [x for (x, b) in zip(sig, np.array(sig) > threshold) if b])
            self.x_list[ind] = np.array(
                [x for (x, b) in zip(te, np.array(sig) > threshold) if b])
            self.p0_list[ind] = np.array(
                [x for (x, b) in zip(p0, np.array(sig) > threshold) if b])

Ancestors

Methods

def threshold_noise(self, threshold=0)

Remove voxel values below a certain threshold from the fitting process, useful if long echo times have been collected and thus thermal noise is being measured below a certain threshold rather than the T2 decay.

Parameters

threshold : float, optional
Default 0 The threshold below which to remove values

Inherited members