diff --git a/docs/source/documentation/btensor_scripts.rst b/docs/source/documentation/btensor_scripts.rst
new file mode 100644
index 000000000..6c6c82478
--- /dev/null
+++ b/docs/source/documentation/btensor_scripts.rst
@@ -0,0 +1,20 @@
+Instructions for tensor-valued dMRI scripts (b-tensor)
+======================================================
+
+
+The scripts for multi-encoding multi-shell multi-tissue CSD (memsmt-CSD) are based on P. Karan et al., Bridging the gap between constrained spherical deconvolution and diffusional variance decomposition via tensor-valued diffusion MRI. Medical Image Analysis (2022). We recommend reading it to understand the scope of the memsmt-CSD problem.
+
+If you want to do CSD with b-tensor data, you should start by computing the fiber response functions. This script should run fast (less than 5 minutes on a full brain).
+::
+
+    scil_compute_memsmt_frf.py wm_frf.txt gm_frf.txt csf_frf.txt --in_dwis dwi_linear.nii.gz dwi_planar.nii.gz dwi_spherical.nii.gz --in_bvals dwi_linear.bval dwi_planar.bval dwi_spherical.bval --in_bvecs dwi_linear.bvec dwi_planar.bvec dwi_spherical.bvec --in_bdeltas 1 -0.5 0 --mask mask.nii.gz --mask_wm wm_mask.nii.gz --mask_gm gm_mask.nii.gz --mask_csf csf_mask.nii.gz -f
+
+Then, you should compute the fODFs and volume fractions. The following command will save a fODF file for each tissue and a volume fractions file. This script should run in about 1-2 hours for a full brain.
+::
+
+    scil_compute_memsmt_fodf.py wm_frf.txt gm_frf.txt csf_frf.txt --in_dwis dwi_linear.nii.gz dwi_planar.nii.gz dwi_spherical.nii.gz --in_bvals dwi_linear.bval dwi_planar.bval dwi_spherical.bval --in_bvecs dwi_linear.bvec dwi_planar.bvec dwi_spherical.bvec --in_bdeltas 1 -0.5 0 --mask mask.nii.gz --processes 8 -f
+
+If you want to do DIVIDE with b-tensor data, you should use the following command. It will save files for the MD, uFA, OP, MK_I, MK_A and MK_T. This script should run in about 1-2 hours for a full brain.
+::
+
+    scil_compute_divide.py --in_dwis dwi_linear.nii.gz dwi_planar.nii.gz dwi_spherical.nii.gz --in_bvals dwi_linear.bval dwi_planar.bval dwi_spherical.bval --in_bvecs dwi_linear.bvec dwi_planar.bvec dwi_spherical.bvec --in_bdeltas 1 -0.5 0 --mask mask.nii.gz --fa fa.nii.gz --processes 8 -f
\ No newline at end of file
diff --git a/scilpy/image/utils.py b/scilpy/image/utils.py
index e0f177eff..c3d35cba0 100644
--- a/scilpy/image/utils.py
+++ b/scilpy/image/utils.py
@@ -76,6 +76,25 @@ def volume_iterator(img, blocksize=1, start=0, end=0):
             yield list(range(stop, end)), img.dataobj[..., stop:end]
 
 
+def extract_affine(input_files):
+    """Extract the affine from a list of nifti files.
+
+    Parameters
+    ----------
+    input_files : list of strings (file paths)
+        Diffusion data files.
+
+    Returns
+    -------
+    affine : np.ndarray
+        Affine of the nifti volume.
+    """
+    for input_file in input_files:
+        if input_file:
+            vol = nib.load(input_file)
+            return vol.affine
+
+
 def check_slice_indices(vol_img, axis_name, slice_ids):
     """Check that the given volume can be sliced at the given slice indices
     along the requested axis.
diff --git a/scilpy/io/fetcher.py b/scilpy/io/fetcher.py
index fb4277f52..77ce86a9d 100644
--- a/scilpy/io/fetcher.py
+++ b/scilpy/io/fetcher.py
@@ -106,6 +106,9 @@ def get_testing_files_dict():
             'anatomical_filtering.zip':
             ['1Li8DdySnMnO9Gich4pilhXisjkjz1-Dy',
              '6f0eff5154ff0973a3dc26db00e383ea'],
+            'btensor_testdata.zip':
+            ['1AMsKlbOZyPnT9TAbxcFzHS1b29aJWKDg',
+             '7c68524fca01268203dc8bfee340f037'],
             'fodf_filtering.zip':
             ['1iyoX2ltLOoLer-v-49LHOzopHCFZ_Tv6',
              'e79c4291af584fdb25814aa7b403a6ce']}
diff --git a/scilpy/reconst/b_tensor_utils.py b/scilpy/reconst/b_tensor_utils.py
new file mode 100644
index 000000000..0ea5d1cb6
--- /dev/null
+++ b/scilpy/reconst/b_tensor_utils.py
@@ -0,0 +1,155 @@
+import logging
+
+from dipy.core.gradients import (gradient_table,
+                                 unique_bvals_tolerance, get_bval_indices)
+from dipy.io.gradients import read_bvals_bvecs
+import nibabel as nib
+import numpy as np
+
+from scilpy.utils.bvec_bval_tools import (normalize_bvecs, is_normalized_bvecs,
+                                          extract_dwi_shell)
+
+
+bshapes = {0: "STE", 1: "LTE", -0.5: "PTE", 0.5: "CTE"}
+
+
+def generate_btensor_input(in_dwis, in_bvals, in_bvecs,
+                           in_bdeltas, force_b0_threshold,
+                           do_pa_signals=False, tol=20):
+    """Generate b-tensor input from an ensemble of data, bvals and bvecs files.
+    This generated input is mandatory for all scripts using b-tensor encoding
+    data. Also generate the powder-averaged (PA) data if set.
+
+    Parameters
+    ----------
+    in_dwis : list of strings (file paths)
+        Diffusion data files for each b-tensor encodings.
+    in_bvals : list of strings (file paths)
+        All of the bvals files associated.
+    in_bvecs : list of strings (file paths)
+        All of the bvecs files associated.
+    in_bdeltas : list of floats
+        All of the b_deltas (describing the type of encoding) files associated.
+    force_b0_threshold : bool, optional
+        If set, will continue even if the minimum bvalue is suspiciously high.
+    do_pa_signals : bool, optional
+        If set, will compute the powder_averaged input instead of the regular
+        one. This means that the signal is averaged over all directions for
+        each bvals.
+    tol : int
+        tolerance gap for b-values clustering. Defaults to 20
+
+    Returns
+    -------
+    gtab_full : GradientTable
+        A single gradient table containing the information of all encodings.
+    data_full : np.ndarray (4d)
+        All concatenated diffusion data from the different encodings.
+    ubvals_full : array
+        All the unique bvals from the different encodings, but with a single
+        b0. If two or more encodings have the same bvalue, then they are
+        differentiate by +1.
+    ub_deltas_full : array
+        All the b_delta values associated with `ubvals_full`.
+    pa_signals : np.ndarray (4d) (if `do_pa_signals`)
+        Powder-averaged diffusion data.
+    gtab_infos : np.ndarray (if `do_pa_signals`)
+        Contains information about the gtab, such as the unique bvals, the
+        encoding types, the number of directions and the acquisition index.
+    """
+    data_full = np.empty(0)
+    bvals_full = np.empty(0)
+    bvecs_full = np.empty(0)
+    b_shapes = np.empty(0)
+    ubvals_full = np.empty(0)
+    ub_deltas_full = np.empty(0)
+    nb_bvecs_full = np.empty(0)
+    acq_index_full = np.empty(0)
+    ubvals_divide = np.empty(0)
+    acq_index = 0
+    for inputf, bvalsf, bvecsf, b_delta in zip(in_dwis, in_bvals,
+                                               in_bvecs, in_bdeltas):
+        if inputf:  # verifies if the input file exists
+            vol = nib.load(inputf)
+            bvals, bvecs = read_bvals_bvecs(bvalsf, bvecsf)
+            if np.sum([bvals > tol]) != 0:
+                bvals = np.round(bvals)
+            if not is_normalized_bvecs(bvecs):
+                logging.warning('Your b-vectors do not seem normalized...')
+                bvecs = normalize_bvecs(bvecs)
+            ubvals = unique_bvals_tolerance(bvals, tol=tol)
+            for ubval in ubvals:  # Loop over all unique bvals
+                # Extracting the data for the ubval shell
+                indices, shell_data, _, output_bvecs = \
+                    extract_dwi_shell(vol, bvals, bvecs, [ubval], tol=tol)
+                nb_bvecs = len(indices)
+                # Adding the current data to each arrays of interest
+                acq_index_full = np.concatenate([acq_index_full,
+                                                 [acq_index]]) \
+                    if acq_index_full.size else np.array([acq_index])
+                ubvals_divide = np.concatenate([ubvals_divide, [ubval]]) \
+                    if ubvals_divide.size else np.array([ubval])
+                while np.isin(ubval, ubvals_full):  # Differentiate the bvals
+                    ubval += 1
+                ubvals_full = np.concatenate([ubvals_full, [ubval]]) \
+                    if ubvals_full.size else np.array([ubval])
+                ub_deltas_full = np.concatenate([ub_deltas_full, [b_delta]]) \
+                    if ub_deltas_full.size else np.array([b_delta])
+                nb_bvecs_full = np.concatenate([nb_bvecs_full, [nb_bvecs]]) \
+                    if nb_bvecs_full.size else np.array([nb_bvecs])
+                data_full = np.concatenate([data_full, shell_data], axis=-1) \
+                    if data_full.size else shell_data
+                bvals_full = np.concatenate([bvals_full,
+                                             np.repeat([ubval], nb_bvecs)]) \
+                    if bvals_full.size else np.repeat([ubval], nb_bvecs)
+                bvecs_full = np.concatenate([bvecs_full, output_bvecs]) \
+                    if bvecs_full.size else output_bvecs
+                b_shapes = np.concatenate([b_shapes,
+                                           np.repeat([bshapes[b_delta]],
+                                                     nb_bvecs)]) \
+                    if b_shapes.size else np.repeat([bshapes[b_delta]],
+                                                    nb_bvecs)
+            acq_index += 1
+    # In the case that the PA data is wanted, there is a different return
+    if do_pa_signals:
+        pa_signals = np.zeros(((data_full.shape[:-1])+(len(ubvals_full),)))
+        for i, ubval in enumerate(ubvals_full):
+            indices = get_bval_indices(bvals_full, ubval, tol=0)
+            pa_signals[..., i] = np.nanmean(data_full[..., indices], axis=-1)
+        gtab_infos = np.ndarray((4, len(ubvals_full)))
+        gtab_infos[0] = ubvals_divide
+        gtab_infos[1] = ub_deltas_full
+        gtab_infos[2] = nb_bvecs_full
+        gtab_infos[3] = acq_index_full
+        if np.sum([ubvals_full < tol]) < acq_index - 1:
+            gtab_infos[3] *= 0
+        return(pa_signals, gtab_infos)
+    # Removing the duplicate b0s from ubvals_full
+    duplicate_b0_ind = np.union1d(np.argwhere(ubvals_full == min(ubvals_full)),
+                                  np.argwhere(ubvals_full > tol))
+    ubvals_full = ubvals_full[duplicate_b0_ind]
+    ub_deltas_full = ub_deltas_full[duplicate_b0_ind]
+    # Sorting the data by bvals
+    sorted_indices = np.argsort(bvals_full, axis=0)
+    bvals_full = np.take_along_axis(bvals_full, sorted_indices, axis=0)
+    bvals_full[bvals_full < tol] = min(ubvals_full)
+    bvecs_full = np.take_along_axis(bvecs_full,
+                                    sorted_indices.reshape(len(bvals_full), 1),
+                                    axis=0)
+    b_shapes = np.take_along_axis(b_shapes, sorted_indices, axis=0)
+    data_full = np.take_along_axis(data_full,
+                                   sorted_indices.reshape(1, 1, 1,
+                                                          len(bvals_full)),
+                                   axis=-1)
+    # Sorting the ubvals
+    sorted_indices = np.argsort(np.asarray(ubvals_full), axis=0)
+    ubvals_full = np.take_along_axis(np.asarray(ubvals_full), sorted_indices,
+                                     axis=0)
+    ub_deltas_full = np.take_along_axis(np.asarray(ub_deltas_full),
+                                        sorted_indices, axis=0)
+    # Creating the corresponding gtab
+    gtab_full = gradient_table(bvals_full, bvecs_full,
+                               b0_threshold=bvals_full.min(),
+                               btens=b_shapes)
+
+    return(gtab_full, data_full, ubvals_full, ub_deltas_full)
diff --git a/scilpy/reconst/divide_fit.py b/scilpy/reconst/divide_fit.py
new file mode 100644
index 000000000..794901886
--- /dev/null
+++ b/scilpy/reconst/divide_fit.py
@@ -0,0 +1,259 @@
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.special import erf
+
+
+def get_bounds():
+    """Define the lower (lb) and upper (ub) boundaries of the fitting
+    parameters, being the signal without diffusion weighting (S0), the mean
+    diffusivity (MD), the isotropic variance (V_I) and the anisotropic variance
+    (V_A).
+
+    Returns
+    -------
+    lb : list of floats
+        Lower boundaries of the fitting parameters.
+    ub : list of floats
+        Upper boundaries of the fitting parameters.
+    """
+    S0 = [0, 10]
+    MD = [1e-12, 4e-9]
+    V_I = [1e-24, 5e-18]
+    V_A = [1e-24, 5e-18]
+    lb = [S0[0], MD[0], V_I[0], V_A[0]]
+    ub = [S0[1], MD[1], V_I[1], V_A[1]]
+    return lb, ub
+
+
+def random_p0(signal, gtab_infos, lb, ub, weight, n_iter):
+    """Produce a guess of initial parameters for the fit, by calculating the
+    signals of a given number of random sets of parameters and keeping the one
+    closest to the input signal.
+
+    Parameters
+    ----------
+    signal : np.array
+        Diffusion data of a single voxel.
+    gtab_infos : np.ndarray
+        Contains information about the gtab, such as the unique bvals, the
+        encoding types, the number of directions and the acquisition index.
+        Obtained as output of the function
+        `reconst.b_tensor_utils.generate_btensor_input`.
+    lb : list of floats
+        Lower boundaries of the fitting parameters.
+    ub : list of floats
+        Upper boundaries of the fitting parameters.
+    weight : np.array
+        Gives a different weight to each element of `signal`.
+    n_iter : int
+        Number of random sets of parameters tested.
+
+    Returns
+    -------
+    guess : np.array
+        Array containing the guessed initial parameters.
+    """
+    guess = []
+    thr = np.inf
+
+    for i in range(n_iter):
+        params_rand = lb + (ub - lb) * np.random.rand(len(lb))
+        signal_rand = gamma_fit2data(gtab_infos, params_rand)
+        residual_rand = np.sum(((signal - signal_rand) * weight)**2)
+
+        if residual_rand < thr:
+            thr = residual_rand
+            guess = params_rand
+
+    return guess
+
+
+def gamma_data2fit(signal, gtab_infos, fit_iters=1, random_iters=50,
+                   do_weight_bvals=False, do_weight_pa=False,
+                   do_multiple_s0=False):
+    """Fit the gamma model to data
+
+    Parameters
+    ----------
+    signal : np.array
+        Diffusion data of a single voxel.
+    gtab_infos : np.ndarray
+        Contains information about the gtab, such as the unique bvals, the
+        encoding types, the number of directions and the acquisition index.
+        Obtained as output of the function
+        `reconst.b_tensor_utils.generate_btensor_input`.
+    fit_iters : int, optional
+        Number of iterations in the gamma fit. Defaults to 1.
+    random_iters : int, optional
+        Number of random sets of parameters tested to find the initial
+        parameters. Defaults to 50.
+    do_weight_bvals : bool , optional
+        If set, does a weighting on the bvalues in the gamma fit.
+    do_weight_pa : bool, optional
+        If set, does a powder averaging weighting in the gamma fit.
+    do_multiple_s0 : bool, optional
+        If set, takes into account multiple baseline signals.
+
+    Returns
+    -------
+    best_params : np.array
+        Array containing the parameters of the fit.
+    """
+    if np.sum(gtab_infos[3]) > 0 and do_multiple_s0 is True:
+        ns = len(np.unique(gtab_infos[3])) - 1
+    else:
+        ns = 0
+
+    unit_to_SI = np.array([np.max(signal), 1e-9, 1e-18, 1e-18])
+    unit_to_SI = np.concatenate((unit_to_SI, np.ones(ns)))
+
+    def weight_bvals(sthr, mdthr, wthr):
+        """Compute an array weighting the different components of the signal
+        array based on the bvalue.
+        """
+        bthr = -np.log(sthr) / mdthr
+        weight = 0.5 * (1 - erf(wthr * (gtab_infos[0] - bthr) / bthr))
+        return weight
+
+    def weight_pa():
+        """Compute an array weighting the different components of the signal
+        array based on the number of directions.
+        """
+        weight = np.sqrt(gtab_infos[2] / np.max(gtab_infos[2]))
+        return weight
+
+    def my_gamma_fit2data(gtab_infos, *args):
+        """Compute a signal from gtab infomations and fit parameters.
+        """
+        params_unit = args
+        params_SI = params_unit * unit_to_SI
+        signal = gamma_fit2data(gtab_infos, params_SI)
+        return signal * weight
+
+    lb_SI, ub_SI = get_bounds()
+    lb_SI = np.concatenate((lb_SI, 0.5 * np.ones(ns)))
+    ub_SI = np.concatenate((ub_SI, 2.0 * np.ones(ns)))
+    lb_SI[0] *= np.max(signal)
+    ub_SI[0] *= np.max(signal)
+
+    lb_unit = lb_SI / unit_to_SI
+    ub_unit = ub_SI / unit_to_SI
+
+    bounds_unit = ([lb_unit, ub_unit])
+
+    res_thr = np.inf
+
+    for i in range(fit_iters):
+        weight = np.ones(len(signal))
+        if do_weight_bvals:
+            weight *= weight_bvals(0.07, 1e-9, 2)
+        if do_weight_pa:
+            weight *= weight_pa()
+
+        p0_SI = random_p0(signal, gtab_infos, lb_SI, ub_SI, weight,
+                          random_iters)
+        p0_unit = p0_SI / unit_to_SI
+        params_unit, params_cov = curve_fit(my_gamma_fit2data, gtab_infos,
+                                            signal * weight, p0=p0_unit,
+                                            bounds=bounds_unit, method="trf",
+                                            ftol=1e-8, xtol=1e-8, gtol=1e-8)
+
+        if do_weight_bvals:
+            weight = weight_bvals(0.07, params_unit[1] * unit_to_SI[1], 2)
+            if do_weight_pa:
+                weight *= weight_pa()
+
+            params_unit, params_cov = curve_fit(my_gamma_fit2data, gtab_infos,
+                                                signal * weight,
+                                                p0=params_unit,
+                                                bounds=bounds_unit,
+                                                method="trf")
+
+        signal_fit = gamma_fit2data(gtab_infos, params_unit * unit_to_SI)
+        residual = np.sum(((signal - signal_fit) * weight) ** 2)
+        if residual < res_thr:
+            res_thr = residual
+            params_best = params_unit
+            # params_cov_best = params_cov
+
+    params_best[0] = params_best[0] * unit_to_SI[0]
+    return params_best[0:4]
+
+
+def gamma_fit2data(gtab_infos, params):
+    """Compute a signal from gtab infomations and fit parameters.
+
+    Parameters
+    ----------
+    gtab_infos : np.ndarray
+        Contains information about the gtab, such as the unique bvals, the
+        encoding types, the number of directions and the acquisition index.
+        Obtained as output of the function
+        `reconst.b_tensor_utils.generate_btensor_input`.
+    params : np.array
+        Array containing the parameters of the fit.
+
+    Returns
+    -------
+    signal : np.array
+        Array containing the signal produced by the gamma model.
+    """
+    S0 = params[0]
+    MD = params[1]
+    V_I = params[2]
+    V_A = params[3]
+    RS = params[4:]  # relative signal
+    if len(RS) != 0:
+        RS = np.concatenate(([1], RS))
+        RS_tile = np.tile(RS, len(gtab_infos[0])).reshape((len(gtab_infos[0]),
+                                                           len(RS)))
+        RS_index = np.zeros((len(gtab_infos[0]), len(RS)))
+        for i in range(len(gtab_infos[0])):
+            j = gtab_infos[3][i]
+            RS_index[i][int(j)] = 1
+        RS_matrix = RS_tile * RS_index
+        SW = S0 * np.sum(RS_matrix, axis=1)
+    else:
+        SW = S0
+
+    V_D = V_I + V_A * (gtab_infos[1] ** 2)
+    signal = SW * ((1 + gtab_infos[0] * V_D / MD) ** (-(MD ** 2) / V_D))
+
+    return np.real(signal)
+
+
+def gamma_fit2metrics(params):
+    """Compute metrics from fit parameters. This is the only function that
+    takes the full brain.
+
+    Parameters
+    ----------
+    params : np.ndarray
+        Array containing the parameters of the fit for the whole brain.
+
+    Returns
+    -------
+    microFA : np.ndarray
+        MicroFA values for the whole brain.
+    MK_I : np.ndarray
+        Isotropic mean kurtosis values for the whole brain.
+    MK_A : np.ndarray
+        Anisotropic mean kurtosis values for the whole brain.
+    MK_T : np.ndarray
+        Total mean kurtosis values for the whole brain.
+    """
+    # S0 = params[..., 0]
+    MD = params[..., 1]
+    V_I = params[..., 2]
+    V_A = params[..., 3]
+    V_T = V_I + V_A
+    V_L = 5 / 2. * V_A
+
+    MK_I = 3 * V_I / (MD ** 2)
+    MK_A = 3 * V_A / (MD ** 2)
+    MK_T = 3 * V_T / (MD ** 2)
+    microFA2 = (3/2.) * (V_L / (V_I + V_L + (MD ** 2)))
+    microFA = np.real(np.sqrt(microFA2))
+    microFA[np.isnan(microFA)] = 0
+
+    return microFA, MK_I, MK_A, MK_T
diff --git a/scilpy/reconst/frf.py b/scilpy/reconst/frf.py
index 03daddfbd..92da02243 100644
--- a/scilpy/reconst/frf.py
+++ b/scilpy/reconst/frf.py
@@ -94,8 +94,8 @@ def compute_ssst_frf(data, bvals, bvecs, mask=None, mask_wm=None,
             "estimation of the fiber response function to ensure no invalid "
             "voxel was used.")
 
-    # Iteratively trying to fit at least min_nvox voxels. Lower the FA threshold
-    # when it doesn't work. Fail if the fa threshold is smaller than
+    # Iteratively trying to fit at least min_nvox voxels. Lower the FA
+    # threshold when it doesn't work. Fail if the fa threshold is smaller than
     # the min_threshold.
     # We use an epsilon since the -= 0.05 might incur numerical imprecision.
     nvox = 0
@@ -132,16 +132,17 @@ def compute_ssst_frf(data, bvals, bvecs, mask=None, mask_wm=None,
     return full_response
 
 
-def compute_msmt_frf(data, bvals, bvecs, data_dti=None, bvals_dti=None,
-                     bvecs_dti=None, mask=None, mask_wm=None, mask_gm=None,
-                     mask_csf=None, fa_thr_wm=0.7, fa_thr_gm=0.2,
-                     fa_thr_csf=0.1, md_thr_gm=0.0007, md_thr_csf=0.003,
-                     min_nvox=300, roi_radii=10, roi_center=None,
-                     tol=20, force_b0_threshold=False):
-    """Compute a single-shell (under b=1500), single-tissue single Fiber
+def compute_msmt_frf(data, bvals, bvecs, btens=None, data_dti=None,
+                     bvals_dti=None, bvecs_dti=None, btens_dti=None,
+                     mask=None, mask_wm=None, mask_gm=None, mask_csf=None,
+                     fa_thr_wm=0.7, fa_thr_gm=0.2, fa_thr_csf=0.1,
+                     md_thr_gm=0.0007, md_thr_csf=0.003, min_nvox=300,
+                     roi_radii=10, roi_center=None, tol=20,
+                     force_b0_threshold=False):
+    """Compute a multi-shell, multi-tissue single Fiber
     Response Function from a DWI volume.
     A DTI fit is made, and voxels containing a single fiber population are
-    found using a threshold on the FA.
+    found using a threshold on the FA and MD.
 
     Parameters
     ----------
@@ -222,7 +223,7 @@ def compute_msmt_frf(data, bvals, bvecs, data_dti=None, bvals_dti=None,
 
     b0_thr = check_b0_threshold(force_b0_threshold, bvals.min(), bvals.min())
 
-    gtab = gradient_table(bvals, bvecs, b0_threshold=b0_thr)
+    gtab = gradient_table(bvals, bvecs, btens=btens, b0_threshold=b0_thr)
 
     if data_dti is None and bvals_dti is None and bvecs_dti is None:
         logging.warning(
@@ -237,14 +238,14 @@ def compute_msmt_frf(data, bvals, bvecs, data_dti=None, bvals_dti=None,
                                      csf_fa_thr=fa_thr_csf,
                                      gm_md_thr=md_thr_gm,
                                      csf_md_thr=md_thr_csf)
-    elif data_dti is not None and bvals_dti is not None and bvecs_dti is not None:
+    elif (data_dti is not None and bvals_dti is not None
+          and bvecs_dti is not None):
         if not is_normalized_bvecs(bvecs_dti):
             logging.warning('Your b-vectors do not seem normalized...')
             bvecs_dti = normalize_bvecs(bvecs_dti)
 
-        b0_thr = check_b0_threshold(
-            force_b0_threshold, bvals_dti.min(), bvals_dti.min())
-        gtab_dti = gradient_table(bvals_dti, bvecs_dti, b0_threshold=b0_thr)
+        check_b0_threshold(force_b0_threshold, bvals_dti.min())
+        gtab_dti = gradient_table(bvals_dti, bvecs_dti, btens=btens_dti)
 
         wm_frf_mask, gm_frf_mask, csf_frf_mask \
             = mask_for_response_msmt(gtab_dti, data_dti,
diff --git a/scilpy/reconst/multi_processes.py b/scilpy/reconst/multi_processes.py
index ca5a7beb1..492ad4952 100644
--- a/scilpy/reconst/multi_processes.py
+++ b/scilpy/reconst/multi_processes.py
@@ -5,13 +5,15 @@
 from scilpy.direction.peaks import peak_directions_asym
 from scipy.sparse.linalg import ArpackNoConvergence
 from dipy.direction.peaks import peak_directions
-from dipy.reconst.mcsd import MSDeconvFit
 from dipy.reconst.multi_voxel import MultiVoxelFit
 from dipy.reconst.odf import gfa
-from dipy.reconst.shm import order_from_ncoef, sh_to_sf_matrix
-from dipy.utils.optpkg import optional_package
+from dipy.reconst.shm import sh_to_sf_matrix, order_from_ncoef
+from dipy.reconst.mcsd import MSDeconvFit
 import numpy as np
 
+from scilpy.reconst.divide_fit import gamma_data2fit
+
+from dipy.utils.optpkg import optional_package
 cvx, have_cvxpy, _ = optional_package("cvxpy")
 
 
@@ -28,9 +30,6 @@ def fit_from_model_parallel(args):
             except cvx.error.SolverError:
                 coeff = np.full((len(model.n)), np.NaN)
                 sub_fit_array[i] = MSDeconvFit(model, coeff, None)
-            except ArpackNoConvergence:
-                coeff = np.full((len(model.n)), np.NaN)
-                sub_fit_array[i] = MSDeconvFit(model, coeff, None)
 
     return chunk_id, sub_fit_array
 
@@ -85,7 +84,7 @@ def fit_from_model(model, data, mask=None, nbr_processes=None):
     fit_array = np.zeros(data_shape[0:3], dtype='object')
     tmp_fit_array = np.zeros((np.count_nonzero(mask)), dtype='object')
     for i, fit in results:
-        tmp_fit_array[chunk_len[i]:chunk_len[i + 1]] = fit
+        tmp_fit_array[chunk_len[i]:chunk_len[i+1]] = fit
 
     fit_array[mask] = tmp_fit_array
     fit_array = MultiVoxelFit(model, fit_array, mask)
@@ -232,9 +231,9 @@ def peaks_from_sh(shm_coeff, sphere, mask=None, relative_peak_threshold=0.5,
     tmp_peak_values_array = np.zeros((np.count_nonzero(mask), npeaks))
     tmp_peak_indices_array = np.zeros((np.count_nonzero(mask), npeaks))
     for i, peak_dirs, peak_values, peak_indices in results:
-        tmp_peak_dirs_array[chunk_len[i]:chunk_len[i + 1], :, :] = peak_dirs
-        tmp_peak_values_array[chunk_len[i]:chunk_len[i + 1], :] = peak_values
-        tmp_peak_indices_array[chunk_len[i]:chunk_len[i + 1], :] = peak_indices
+        tmp_peak_dirs_array[chunk_len[i]:chunk_len[i+1], :, :] = peak_dirs
+        tmp_peak_values_array[chunk_len[i]:chunk_len[i+1], :] = peak_values
+        tmp_peak_indices_array[chunk_len[i]:chunk_len[i+1], :] = peak_indices
 
     peak_dirs_array[mask] = tmp_peak_dirs_array
     peak_values_array[mask] = tmp_peak_values_array
@@ -386,12 +385,12 @@ def maps_from_sh(shm_coeff, peak_dirs, peak_values, peak_indices, sphere,
         all_time_max_odf = max(all_time_global_max, max_odf)
         all_time_global_max = max(all_time_global_max, global_max)
 
-        tmp_nufo_map_array[chunk_len[i]:chunk_len[i + 1]] = nufo_map
-        tmp_afd_max_array[chunk_len[i]:chunk_len[i + 1]] = afd_max
-        tmp_afd_sum_array[chunk_len[i]:chunk_len[i + 1]] = afd_sum
-        tmp_rgb_map_array[chunk_len[i]:chunk_len[i + 1], :] = rgb_map
-        tmp_gfa_map_array[chunk_len[i]:chunk_len[i + 1]] = gfa_map
-        tmp_qa_map_array[chunk_len[i]:chunk_len[i + 1], :] = qa_map
+        tmp_nufo_map_array[chunk_len[i]:chunk_len[i+1]] = nufo_map
+        tmp_afd_max_array[chunk_len[i]:chunk_len[i+1]] = afd_max
+        tmp_afd_sum_array[chunk_len[i]:chunk_len[i+1]] = afd_sum
+        tmp_rgb_map_array[chunk_len[i]:chunk_len[i+1], :] = rgb_map
+        tmp_gfa_map_array[chunk_len[i]:chunk_len[i+1]] = gfa_map
+        tmp_qa_map_array[chunk_len[i]:chunk_len[i+1], :] = qa_map
 
     nufo_map_array[mask] = tmp_nufo_map_array
     afd_max_array[mask] = tmp_afd_max_array
@@ -409,8 +408,8 @@ def maps_from_sh(shm_coeff, peak_dirs, peak_values, peak_indices, sphere,
             or np.array_equal(np.array([1]), afd_unique):
         logging.warning('All AFD_max values are 1. The peaks seem normalized.')
 
-    return (nufo_map_array, afd_max_array, afd_sum_array,
-            rgb_map_array, gfa_map_array, qa_map_array)
+    return(nufo_map_array, afd_max_array, afd_sum_array,
+           rgb_map_array, gfa_map_array, qa_map_array)
 
 
 def convert_sh_basis_parallel(args):
@@ -489,7 +488,7 @@ def convert_sh_basis(shm_coeff, sphere, mask=None,
     shm_coeff_array = np.zeros(data_shape)
     tmp_shm_coeff_array = np.zeros((np.count_nonzero(mask), data_shape[3]))
     for i, new_shm_coeff in results:
-        tmp_shm_coeff_array[chunk_len[i]:chunk_len[i + 1], :] = new_shm_coeff
+        tmp_shm_coeff_array[chunk_len[i]:chunk_len[i+1], :] = new_shm_coeff
 
     shm_coeff_array[mask] = tmp_shm_coeff_array
 
@@ -583,3 +582,98 @@ def convert_sh_to_sf(shm_coeff, sphere, mask=None, dtype="float32",
     sf_array[mask] = tmp_sf_array
 
     return sf_array
+
+
+def fit_gamma_parallel(args):
+    data = args[0]
+    gtab_infos = args[1]
+    fit_iters = args[2]
+    random_iters = args[3]
+    do_weight_bvals = args[4]
+    do_weight_pa = args[5]
+    do_multiple_s0 = args[6]
+    chunk_id = args[7]
+
+    sub_fit_array = np.zeros((data.shape[0], 4))
+    for i in range(data.shape[0]):
+        if data[i].any():
+            sub_fit_array[i] = gamma_data2fit(data[i], gtab_infos, fit_iters,
+                                              random_iters, do_weight_bvals,
+                                              do_weight_pa, do_multiple_s0)
+
+    return chunk_id, sub_fit_array
+
+
+def fit_gamma(data, gtab_infos, mask=None, fit_iters=1, random_iters=50,
+              do_weight_bvals=False, do_weight_pa=False, do_multiple_s0=False,
+              nbr_processes=None):
+    """Fit the gamma model to data
+
+    Parameters
+    ----------
+    data : np.ndarray (4d)
+        Diffusion data, powder averaged. Obtained as output of the function
+        `reconst.b_tensor_utils.generate_powder_averaged_data`.
+    gtab_infos : np.ndarray
+        Contains information about the gtab, such as the unique bvals, the
+        encoding types, the number of directions and the acquisition index.
+        Obtained as output of the function
+        `reconst.b_tensor_utils.generate_powder_averaged_data`.
+    mask : np.ndarray, optional
+        If `mask` is provided, only the data inside the mask will be
+        used for computations.
+    fit_iters : int, optional
+        Number of iterations in the gamma fit. Defaults to 1.
+    random_iters : int, optional
+        Number of random sets of parameters tested to find the initial
+        parameters. Defaults to 50.
+    do_weight_bvals : bool , optional
+        If set, does a weighting on the bvalues in the gamma fit.
+    do_weight_pa : bool, optional
+        If set, does a powder averaging weighting in the gamma fit.
+    do_multiple_s0 : bool, optional
+        If set, takes into account multiple baseline signals.
+    nbr_processes : int, optional
+        The number of subprocesses to use.
+        Default: multiprocessing.cpu_count()
+
+    Returns
+    -------
+    fit_array : np.ndarray
+        Array containing the fit
+    """
+    data_shape = data.shape
+    if mask is None:
+        mask = np.sum(data, axis=3).astype(bool)
+
+    nbr_processes = multiprocessing.cpu_count() if nbr_processes is None \
+        or nbr_processes <= 0 else nbr_processes
+
+    # Ravel the first 3 dimensions while keeping the 4th intact, like a list of
+    # 1D time series voxels. Then separate it in chunks of len(nbr_processes).
+    data = data[mask].reshape((np.count_nonzero(mask), data_shape[3]))
+    chunks = np.array_split(data, nbr_processes)
+
+    chunk_len = np.cumsum([0] + [len(c) for c in chunks])
+    pool = multiprocessing.Pool(nbr_processes)
+    results = pool.map(fit_gamma_parallel,
+                       zip(chunks,
+                           itertools.repeat(gtab_infos),
+                           itertools.repeat(fit_iters),
+                           itertools.repeat(random_iters),
+                           itertools.repeat(do_weight_bvals),
+                           itertools.repeat(do_weight_pa),
+                           itertools.repeat(do_multiple_s0),
+                           np.arange(len(chunks))))
+    pool.close()
+    pool.join()
+
+    # Re-assemble the chunk together in the original shape.
+    fit_array = np.zeros((data_shape[0:3])+(4,))
+    tmp_fit_array = np.zeros((np.count_nonzero(mask), 4))
+    for i, fit in results:
+        tmp_fit_array[chunk_len[i]:chunk_len[i+1]] = fit
+
+    fit_array[mask] = tmp_fit_array
+
+    return fit_array
diff --git a/scripts/scil_compute_divide.py b/scripts/scil_compute_divide.py
new file mode 100644
index 000000000..0ef609a66
--- /dev/null
+++ b/scripts/scil_compute_divide.py
@@ -0,0 +1,226 @@
+#! /usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+Script to compute microstructure metrics using the DIVIDE method. In order to
+operate, the script needs at leats two different types of b-tensor encodings.
+Note that custom encodings are not yet supported, so that only the linear
+tensor encoding (LTE, b_delta = 1), the planar tensor encoding
+(PTE, b_delta = -0.5), the spherical tensor encoding (STE, b_delta = 0) and
+the cigar shape tensor encoding (b_delta = 0.5) are available. Moreover, all
+of `--in_dwis`, `--in_bvals`, `--in_bvecs` and `--in_bdeltas` must have the
+same number of arguments. Be sure to keep the same order of encodings
+throughout all these inputs and to set `--in_bdeltas` accordingly (IMPORTANT).
+
+By default, will output all possible files, using default names.
+Specific names can be specified using the file flags specified in the
+"File flags" section.
+
+If --not_all is set, only the files specified explicitly by the flags
+will be output.
+
+Based on Markus Nilsson, Filip Szczepankiewicz, Björn Lampinen, André Ahlgren,
+João P. de Almeida Martins, Samo Lasic, Carl-Fredrik Westin,
+and Daniel Topgaard. An open-source framework for analysis of multidimensional
+diffusion MRI data implemented in MATLAB.
+Proc. Intl. Soc. Mag. Reson. Med. (26), Paris, France, 2018.
+"""
+
+import argparse
+import logging
+
+import nibabel as nib
+import numpy as np
+
+from scilpy.image.utils import extract_affine
+from scilpy.io.image import get_data_as_mask
+from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
+                             assert_outputs_exist, add_force_b0_arg,
+                             add_processes_arg, add_verbose_arg)
+from scilpy.reconst.multi_processes import fit_gamma
+from scilpy.reconst.divide_fit import gamma_fit2metrics
+from scilpy.reconst.b_tensor_utils import generate_btensor_input
+
+
+def _build_arg_parser():
+    p = argparse.ArgumentParser(description=__doc__,
+                                formatter_class=argparse.RawTextHelpFormatter)
+
+    p.add_argument('--in_dwis', nargs='+', required=True,
+                   help='Path to the input diffusion volume for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bvals', nargs='+', required=True,
+                   help='Path to the bval file, in FSL format, for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bvecs', nargs='+', required=True,
+                   help='Path to the bvec file, in FSL format, for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bdeltas', nargs='+', type=float,
+                   choices=[0, 1, -0.5, 0.5], required=True,
+                   help='Value of b_delta for each b-tensor encoding type, '
+                        'in the same order as dwi, bval and bvec inputs.')
+
+    p.add_argument(
+        '--mask',
+        help='Path to a binary mask. Only the data inside the '
+             'mask will be used for computations and reconstruction.')
+    p.add_argument(
+        '--fa',
+        help='Path to a FA map. Needed for calculating the OP.')
+    p.add_argument(
+        '--tolerance', type=int, default=20,
+        help='The tolerated gap between the b-values to '
+             'extract\nand the current b-value. [%(default)s]')
+    p.add_argument(
+        '--fit_iters', type=int, default=1,
+        help='The number of time the gamma fit will be done [%(default)s]')
+    p.add_argument(
+        '--random_iters', type=int, default=50,
+        help='The number of iterations for the initial parameters search. '
+             '[%(default)s]')
+    p.add_argument(
+        '--do_weight_bvals', action='store_false',
+        help='If set, does not do a weighting on the bvalues in the gamma '
+             'fit.')
+    p.add_argument(
+        '--do_weight_pa', action='store_false',
+        help='If set, does not do a powder averaging weighting in the gamma '
+             'fit.')
+    p.add_argument(
+        '--do_multiple_s0', action='store_false',
+        help='If set, does not take into account multiple baseline signals.')
+
+    add_force_b0_arg(p)
+    add_processes_arg(p)
+    add_overwrite_arg(p)
+    add_verbose_arg(p)
+
+    p.add_argument(
+        '--not_all', action='store_true',
+        help='If set, only saves the files specified using the '
+             'file flags. (Default: False)')
+
+    g = p.add_argument_group(title='File flags')
+
+    g.add_argument(
+        '--md', metavar='file', default='',
+        help='Output filename for the MD.')
+    g.add_argument(
+        '--ufa', metavar='file', default='',
+        help='Output filename for the microscopic FA.')
+    g.add_argument(
+        '--op', metavar='file', default='',
+        help='Output filename for the order parameter.')
+    g.add_argument(
+        '--mk_i', metavar='file', default='',
+        help='Output filename for the isotropic mean kurtosis.')
+    g.add_argument(
+        '--mk_a', metavar='file', default='',
+        help='Output filename for the anisotropic mean kurtosis.')
+    g.add_argument(
+        '--mk_t', metavar='file', default='',
+        help='Output filename for the total mean kurtosis.')
+
+    return p
+
+
+def main():
+    parser = _build_arg_parser()
+    args = parser.parse_args()
+
+    if args.verbose:
+        logging.basicConfig(level=logging.DEBUG)
+    else:
+        logging.basicConfig(level=logging.INFO)
+
+    if not args.not_all:
+        args.md = args.md or 'md.nii.gz'
+        args.ufa = args.ufa or 'ufa.nii.gz'
+        args.op = args.op or 'op.nii.gz'
+        args.mk_i = args.mk_i or 'mk_i.nii.gz'
+        args.mk_a = args.mk_a or 'mk_a.nii.gz'
+        args.mk_t = args.mk_t or 'mk_t.nii.gz'
+
+    arglist = [args.md, args.ufa, args.op, args.mk_i, args.mk_a, args.mk_t]
+    if args.not_all and not any(arglist):
+        parser.error('When using --not_all, you need to specify at least ' +
+                     'one file to output.')
+
+    assert_inputs_exist(parser, [],
+                        optional=list(np.concatenate((args.in_dwis,
+                                                      args.in_bvals,
+                                                      args.in_bvecs))))
+    assert_outputs_exist(parser, args, arglist)
+
+    if not (len(args.in_dwis) == len(args.in_bvals)
+            == len(args.in_bvecs) == len(args.in_bdeltas)):
+        msg = """The number of given dwis, bvals, bvecs and bdeltas must be the
+              same. Please verify that all inputs were correctly inserted."""
+        raise ValueError(msg)
+
+    if len(np.unique(args.in_bdeltas)) < 2:
+        msg = """At least two different b-tensor shapes are needed for the
+              script to work properly."""
+        raise ValueError(msg)
+
+    affine = extract_affine(args.in_dwis)
+
+    tol = args.tolerance
+    force_b0_thr = args.force_b0_threshold
+
+    data, gtab_infos = generate_btensor_input(args.in_dwis,
+                                              args.in_bvals,
+                                              args.in_bvecs,
+                                              args.in_bdeltas,
+                                              force_b0_thr,
+                                              do_pa_signals=True,
+                                              tol=tol)
+
+    gtab_infos[0] *= 1e6  # getting bvalues to SI units
+
+    if args.mask is None:
+        mask = None
+    else:
+        mask = get_data_as_mask(nib.load(args.mask), dtype=bool)
+        if mask.shape != data.shape[:-1]:
+            raise ValueError("Mask is not the same shape as data.")
+
+    if args.fa is not None:
+        vol = nib.load(args.fa)
+        FA = vol.get_fdata(dtype=np.float32)
+
+    parameters = fit_gamma(data, gtab_infos, mask=mask,
+                           fit_iters=args.fit_iters,
+                           random_iters=args.random_iters,
+                           do_weight_bvals=args.do_weight_bvals,
+                           do_weight_pa=args.do_weight_pa,
+                           do_multiple_s0=args.do_multiple_s0,
+                           nbr_processes=args.nbr_processes)
+
+    microFA, MK_I, MK_A, MK_T = gamma_fit2metrics(parameters)
+    microFA[np.isnan(microFA)] = 0
+    microFA = np.clip(microFA, 0, 1)
+
+    if args.md:
+        nib.save(nib.Nifti1Image(parameters[..., 1].astype(np.float32),
+                                 affine), args.md)
+
+    if args.ufa:
+        nib.save(nib.Nifti1Image(microFA.astype(np.float32), affine), args.ufa)
+    if args.op:
+        if args.fa is not None:
+            OP = np.sqrt((3 * (microFA ** (-2)) - 2) / (3 * (FA ** (-2)) - 2))
+            OP[microFA < FA] = 0
+            nib.save(nib.Nifti1Image(OP.astype(np.float32), affine), args.op)
+        else:
+            logging.warning('The FA must be given in order to compute the OP.')
+    if args.mk_i:
+        nib.save(nib.Nifti1Image(MK_I.astype(np.float32), affine), args.mk_i)
+    if args.mk_a:
+        nib.save(nib.Nifti1Image(MK_A.astype(np.float32), affine), args.mk_a)
+    if args.mk_t:
+        nib.save(nib.Nifti1Image(MK_T.astype(np.float32), affine), args.mk_t)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/scil_compute_memsmt_fodf.py b/scripts/scil_compute_memsmt_fodf.py
new file mode 100644
index 000000000..fc0ef0928
--- /dev/null
+++ b/scripts/scil_compute_memsmt_fodf.py
@@ -0,0 +1,375 @@
+#! /usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+Script to compute multi-encoding multi-shell multi-tissue (memsmt)
+Constrained Spherical Deconvolution ODFs. In order to operate,
+the script only needs the data from one type of b-tensor encoding. However,
+giving only a spherical one will not produce good fODFs, as it only probes
+spherical shapes. As for planar encoding, it should technically
+work alone, but seems to be very sensitive to noise and is yet to be properly
+documented. We thus suggest to always use at least the linear encoding, which
+will be equivalent to standard multi-shell multi-tissue if used alone, in
+combinaison with other encodings. Note that custom encodings are not yet
+supported, so that only the linear tensor encoding (LTE, b_delta = 1), the
+planar tensor encoding (PTE, b_delta = -0.5), the spherical tensor encoding
+(STE, b_delta = 0) and the cigar shape tensor encoding (b_delta = 0.5) are
+available. Moreover, all of `--in_dwis`, `--in_bvals`, `--in_bvecs` and
+`--in_bdeltas` must have the same number of arguments. Be sure to keep the
+same order of encodings throughout all these inputs and to set `--in_bdeltas`
+accordingly (IMPORTANT).
+
+By default, will output all possible files, using default names.
+Specific names can be specified using the file flags specified in the
+"File flags" section.
+
+If --not_all is set, only the files specified explicitly by the flags
+will be output.
+
+Based on P. Karan et al., Bridging the gap between constrained spherical 
+deconvolution and diffusional variance decomposition via tensor-valued 
+diffusion MRI. Medical Image Analysis (2022)
+"""
+
+import argparse
+import logging
+
+from dipy.core.gradients import GradientTable
+from dipy.data import get_sphere, default_sphere
+from dipy.reconst import shm
+from dipy.reconst.mcsd import MultiShellResponse, MultiShellDeconvModel
+import nibabel as nib
+import numpy as np
+
+from scilpy.image.utils import extract_affine
+from scilpy.io.image import get_data_as_mask
+from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
+                             assert_outputs_exist, add_force_b0_arg,
+                             add_sh_basis_args, add_processes_arg,
+                             add_verbose_arg)
+from scilpy.reconst.multi_processes import fit_from_model, convert_sh_basis
+from scilpy.reconst.b_tensor_utils import generate_btensor_input
+
+
+def _build_arg_parser():
+    p = argparse.ArgumentParser(description=__doc__,
+                                formatter_class=argparse.RawTextHelpFormatter)
+
+    p.add_argument('in_wm_frf',
+                   help='Text file of WM response function.')
+    p.add_argument('in_gm_frf',
+                   help='Text file of GM response function.')
+    p.add_argument('in_csf_frf',
+                   help='Text file of CSF response function.')
+
+    p.add_argument('--in_dwis', nargs='+', required=True,
+                   help='Path to the input diffusion volume for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bvals', nargs='+', required=True,
+                   help='Path to the bval file, in FSL format, for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bvecs', nargs='+', required=True,
+                   help='Path to the bvec file, in FSL format, for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bdeltas', nargs='+', type=float,
+                   choices=[0, 1, -0.5, 0.5], required=True,
+                   help='Value of b_delta for each b-tensor encoding type, '
+                        'in the same order as dwi, bval and bvec inputs.')
+
+    p.add_argument(
+        '--sh_order', metavar='int', default=8, type=int,
+        help='SH order used for the CSD. (Default: 8)')
+    p.add_argument(
+        '--mask',
+        help='Path to a binary mask. Only the data inside the '
+             'mask will be used for computations and reconstruction.')
+    p.add_argument(
+        '--tolerance', type=int, default=20,
+        help='The tolerated gap between the b-values to '
+             'extract\nand the current b-value. [%(default)s]')
+
+    add_force_b0_arg(p)
+    add_sh_basis_args(p)
+    add_processes_arg(p)
+    add_overwrite_arg(p)
+    add_verbose_arg(p)
+
+    p.add_argument(
+        '--not_all', action='store_true',
+        help='If set, only saves the files specified using the '
+             'file flags. (Default: False)')
+
+    g = p.add_argument_group(title='File flags')
+
+    g.add_argument(
+        '--wm_out_fODF', metavar='file', default='',
+        help='Output filename for the WM fODF coefficients.')
+    g.add_argument(
+        '--gm_out_fODF', metavar='file', default='',
+        help='Output filename for the GM fODF coefficients.')
+    g.add_argument(
+        '--csf_out_fODF', metavar='file', default='',
+        help='Output filename for the CSF fODF coefficients.')
+    g.add_argument(
+        '--vf', metavar='file', default='',
+        help='Output filename for the volume fractions map.')
+    g.add_argument(
+        '--vf_rgb', metavar='file', default='',
+        help='Output filename for the volume fractions map in rgb.')
+
+    return p
+
+
+def single_tensor_btensor(gtab, evals, b_delta, S0=1):
+    # This function should be moved to Dipy at some point
+
+    if b_delta > 1 or b_delta < -0.5:
+        msg = """The value of b_delta must be between -0.5 and 1."""
+        raise ValueError(msg)
+
+    out_shape = gtab.bvecs.shape[:gtab.bvecs.ndim - 1]
+    gradients = gtab.bvecs.reshape(-1, 3)
+
+    evals = np.asarray(evals)
+    D_iso = np.sum(evals) / 3.
+    D_para = evals[np.argmax(abs(evals - D_iso))]
+    D_perp = evals[np.argmin(abs(evals - D_iso))]
+    D_delta = (D_para - D_perp) / (3 * D_iso)
+
+    S = np.zeros(len(gradients))
+    for (i, g) in enumerate(gradients):
+        theta = np.arctan2(np.sqrt(g[0] ** 2 + g[1] ** 2), g[2])
+        P_2 = (3 * np.cos(theta) ** 2 - 1) / 2.
+        b = gtab.bvals[i]
+        S[i] = S0 * np.exp(-b * D_iso * (1 + 2 * b_delta * D_delta * P_2))
+
+    return S.reshape(out_shape)
+
+
+def multi_shell_fiber_response(sh_order, bvals, wm_rf, gm_rf, csf_rf,
+                               b_deltas=None, sphere=None, tol=20):
+    # This function should be moved to Dipy at some point
+
+    bvals = np.array(bvals, copy=True)
+
+    n = np.arange(0, sh_order + 1, 2)
+    m = np.zeros_like(n)
+
+    if sphere is None:
+        sphere = default_sphere
+
+    big_sphere = sphere.subdivide()
+    theta, phi = big_sphere.theta, big_sphere.phi
+
+    B = shm.real_sh_descoteaux_from_index(m, n, theta[:, None], phi[:, None])
+    A = shm.real_sh_descoteaux_from_index(0, 0, 0, 0)
+
+    if b_deltas is None:
+        b_deltas = np.ones(len(bvals) - 1)
+
+    response = np.empty([len(bvals), len(n) + 2])
+
+    if bvals[0] < tol:
+        gtab = GradientTable(big_sphere.vertices * 0)
+        wm_response = single_tensor_btensor(gtab, wm_rf[0, :3], 1, wm_rf[0, 3])
+        response[0, 2:] = np.linalg.lstsq(B, wm_response, rcond=None)[0]
+
+        response[0, 1] = gm_rf[0, 3] / A
+        response[0, 0] = csf_rf[0, 3] / A
+        for i, bvalue in enumerate(bvals[1:]):
+            gtab = GradientTable(big_sphere.vertices * bvalue)
+            wm_response = single_tensor_btensor(gtab, wm_rf[i, :3],
+                                                b_deltas[i],
+                                                wm_rf[i, 3])
+            response[i+1, 2:] = np.linalg.lstsq(B, wm_response, rcond=None)[0]
+
+            response[i+1, 1] = gm_rf[i, 3] * np.exp(-bvalue * gm_rf[i, 0]) / A
+            response[i+1, 0] = csf_rf[i, 3] * np.exp(-bvalue
+                                                     * csf_rf[i, 0]) / A
+
+        S0 = [csf_rf[0, 3], gm_rf[0, 3], wm_rf[0, 3]]
+
+    else:
+        logging.warning('No b0 was given. Proceeding either way.')
+        for i, bvalue in enumerate(bvals):
+            gtab = GradientTable(big_sphere.vertices * bvalue)
+            wm_response = single_tensor_btensor(gtab, wm_rf[i, :3],
+                                                b_deltas[i],
+                                                wm_rf[i, 3])
+            response[i, 2:] = np.linalg.lstsq(B, wm_response, rcond=None)[0]
+
+            response[i, 1] = gm_rf[i, 3] * np.exp(-bvalue * gm_rf[i, 0]) / A
+            response[i, 0] = csf_rf[i, 3] * np.exp(-bvalue * csf_rf[i, 0]) / A
+
+        S0 = [csf_rf[0, 3], gm_rf[0, 3], wm_rf[0, 3]]
+
+    return MultiShellResponse(response, sh_order, bvals, S0=S0)
+
+
+def main():
+    parser = _build_arg_parser()
+    args = parser.parse_args()
+
+    if args.verbose:
+        logging.basicConfig(level=logging.DEBUG)
+    else:
+        logging.basicConfig(level=logging.INFO)
+
+    if not args.not_all:
+        args.wm_out_fODF = args.wm_out_fODF or 'wm_fodf.nii.gz'
+        args.gm_out_fODF = args.gm_out_fODF or 'gm_fodf.nii.gz'
+        args.csf_out_fODF = args.csf_out_fODF or 'csf_fodf.nii.gz'
+        args.vf = args.vf or 'vf.nii.gz'
+        args.vf_rgb = args.vf_rgb or 'vf_rgb.nii.gz'
+
+    arglist = [args.wm_out_fODF, args.gm_out_fODF, args.csf_out_fODF,
+               args.vf, args.vf_rgb]
+    if args.not_all and not any(arglist):
+        parser.error('When using --not_all, you need to specify at least ' +
+                     'one file to output.')
+
+    assert_inputs_exist(parser, [],
+                        optional=list(np.concatenate((args.in_dwis,
+                                                      args.in_bvals,
+                                                      args.in_bvecs))))
+    assert_outputs_exist(parser, args, arglist)
+
+    if not (len(args.in_dwis) == len(args.in_bvals)
+            == len(args.in_bvecs) == len(args.in_bdeltas)):
+        msg = """The number of given dwis, bvals, bvecs and bdeltas must be the
+              same. Please verify that all inputs were correctly inserted."""
+        raise ValueError(msg)
+
+    affine = extract_affine(args.in_dwis)
+
+    tol = args.tolerance
+    force_b0_thr = args.force_b0_threshold
+
+    wm_frf = np.loadtxt(args.in_wm_frf)
+    gm_frf = np.loadtxt(args.in_gm_frf)
+    csf_frf = np.loadtxt(args.in_csf_frf)
+
+    gtab, data, ubvals, ubdeltas = generate_btensor_input(args.in_dwis,
+                                                          args.in_bvals,
+                                                          args.in_bvecs,
+                                                          args.in_bdeltas,
+                                                          force_b0_thr,
+                                                          tol=tol)
+
+    # Checking mask
+    if args.mask is None:
+        mask = None
+    else:
+        mask = get_data_as_mask(nib.load(args.mask), dtype=bool)
+        if mask.shape != data.shape[:-1]:
+            raise ValueError("Mask is not the same shape as data.")
+
+    sh_order = args.sh_order
+
+    # Checking data and sh_order
+    if data.shape[-1] < (sh_order + 1) * (sh_order + 2) / 2:
+        logging.warning(
+            'We recommend having at least {} unique DWIs volumes, but you '
+            'currently have {} volumes. Try lowering the parameter --sh_order '
+            'in case of non convergence.'.format(
+                (sh_order + 1) * (sh_order + 2) / 2, data.shape[-1]))
+
+    # Checking response functions and computing msmt response function
+    if len(wm_frf.shape) == 1:
+        wm_frf = np.reshape(wm_frf, (1,) + wm_frf.shape)
+    if len(gm_frf.shape) == 1:
+        gm_frf = np.reshape(gm_frf, (1,) + gm_frf.shape)
+    if len(csf_frf.shape) == 1:
+        csf_frf = np.reshape(csf_frf, (1,) + csf_frf.shape)
+
+    if not wm_frf.shape[1] == 4:
+        raise ValueError('WM frf file did not contain 4 elements. '
+                         'Invalid or deprecated FRF format')
+    if not gm_frf.shape[1] == 4:
+        raise ValueError('GM frf file did not contain 4 elements. '
+                         'Invalid or deprecated FRF format')
+    if not csf_frf.shape[1] == 4:
+        raise ValueError('CSF frf file did not contain 4 elements. '
+                         'Invalid or deprecated FRF format')
+
+    memsmt_response = multi_shell_fiber_response(sh_order,
+                                                 ubvals,
+                                                 wm_frf, gm_frf, csf_frf,
+                                                 ubdeltas[1:],
+                                                 tol=tol)
+
+    reg_sphere = get_sphere('symmetric362')
+
+    # Computing memsmt-CSD
+    memsmt_model = MultiShellDeconvModel(gtab, memsmt_response,
+                                         reg_sphere=reg_sphere,
+                                         sh_order=sh_order)
+
+    # Computing memsmt-CSD fit
+    memsmt_fit = fit_from_model(memsmt_model, data,
+                                mask=mask, nbr_processes=args.nbr_processes)
+
+    shm_coeff = memsmt_fit.all_shm_coeff
+
+    nan_count = len(np.argwhere(np.isnan(shm_coeff[..., 0])))
+    voxel_count = np.prod(shm_coeff.shape[:-1])
+
+    if nan_count / voxel_count >= 0.05:
+        msg = """There are {} voxels out of {} that could not be solved by
+        the solver, reaching a critical amount of voxels. Make sure to tune the
+        response functions properly, as the solving process is very sensitive
+        to it. Proceeding to fill the problematic voxels by 0.
+        """
+        logging.warning(msg.format(nan_count, voxel_count))
+    elif nan_count > 0:
+        msg = """There are {} voxels out of {} that could not be solved by
+        the solver. Make sure to tune the response functions properly, as the
+        solving process is very sensitive to it. Proceeding to fill the
+        problematic voxels by 0.
+        """
+        logging.warning(msg.format(nan_count, voxel_count))
+
+    shm_coeff = np.where(np.isnan(shm_coeff), 0, shm_coeff)
+
+    vf = memsmt_fit.volume_fractions
+    vf = np.where(np.isnan(vf), 0, vf)
+
+    # Saving results
+    if args.wm_out_fODF:
+        wm_coeff = shm_coeff[..., 2:]
+        if args.sh_basis == 'tournier07':
+            wm_coeff = convert_sh_basis(wm_coeff, reg_sphere, mask=mask,
+                                        nbr_processes=args.nbr_processes)
+        nib.save(nib.Nifti1Image(wm_coeff.astype(np.float32),
+                                 affine), args.wm_out_fODF)
+
+    if args.gm_out_fODF:
+        gm_coeff = shm_coeff[..., 1]
+        if args.sh_basis == 'tournier07':
+            gm_coeff = gm_coeff.reshape(gm_coeff.shape + (1,))
+            gm_coeff = convert_sh_basis(gm_coeff, reg_sphere, mask=mask,
+                                        nbr_processes=args.nbr_processes)
+        nib.save(nib.Nifti1Image(gm_coeff.astype(np.float32),
+                                 affine), args.gm_out_fODF)
+
+    if args.csf_out_fODF:
+        csf_coeff = shm_coeff[..., 0]
+        if args.sh_basis == 'tournier07':
+            csf_coeff = csf_coeff.reshape(csf_coeff.shape + (1,))
+            csf_coeff = convert_sh_basis(csf_coeff, reg_sphere, mask=mask,
+                                         nbr_processes=args.nbr_processes)
+        nib.save(nib.Nifti1Image(csf_coeff.astype(np.float32),
+                                 affine), args.csf_out_fODF)
+
+    if args.vf:
+        nib.save(nib.Nifti1Image(vf.astype(np.float32), affine), args.vf)
+
+    if args.vf_rgb:
+        vf_rgb = vf / np.max(vf) * 255
+        vf_rgb = np.clip(vf_rgb, 0, 255)
+        nib.save(nib.Nifti1Image(vf_rgb.astype(np.uint8),
+                                 affine), args.vf_rgb)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/scil_compute_memsmt_frf.py b/scripts/scil_compute_memsmt_frf.py
new file mode 100644
index 000000000..a8051de31
--- /dev/null
+++ b/scripts/scil_compute_memsmt_frf.py
@@ -0,0 +1,305 @@
+#! /usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+Script to estimate response functions for multi-encoding multi-shell
+multi-tissue (memsmt) constrained spherical deconvolution. In order to operate,
+the script only needs the data from one type of b-tensor encoding. However,
+giving only a spherical one will not produce good fiber response functions, as
+it only probes spherical shapes. As for planar encoding, it should technically
+work alone, but seems to be very sensitive to noise and is yet to be properly
+documented. We thus suggest to always use at least the linear encoding, which
+will be equivalent to standard multi-shell multi-tissue if used alone, in
+combinaison with other encodings. Note that custom encodings are not yet
+supported, so that only the linear tensor encoding (LTE, b_delta = 1), the
+planar tensor encoding (PTE, b_delta = -0.5), the spherical tensor encoding
+(STE, b_delta = 0) and the cigar shape tensor encoding (b_delta = 0.5) are
+available. Moreover, all of `--in_dwis`, `--in_bvals`, `--in_bvecs` and
+`--in_bdeltas` must have the same number of arguments. Be sure to keep the
+same order of encodings throughout all these inputs and to set `--in_bdeltas`
+accordingly (IMPORTANT).
+
+The script computes a response function for white-matter (wm),
+gray-matter (gm), csf and the mean b=0.
+
+In the wm, we compute the response function in each voxels where
+the FA is superior at threshold_fa_wm.
+
+In the gm (or csf), we compute the response function in each voxels where
+the FA is below at threshold_fa_gm (or threshold_fa_csf) and where
+the MD is below threshold_md_gm (or threshold_md_csf).
+
+Based on P. Karan et al., Bridging the gap between constrained spherical
+deconvolution and diffusional variance decomposition via tensor-valued
+diffusion MRI. Medical Image Analysis (2022)
+"""
+
+import argparse
+import logging
+
+import nibabel as nib
+import numpy as np
+
+from scilpy.image.utils import extract_affine
+from scilpy.io.image import get_data_as_mask
+from scilpy.io.utils import (add_force_b0_arg,
+                             add_overwrite_arg, add_verbose_arg,
+                             assert_inputs_exist, assert_outputs_exist)
+from scilpy.reconst.frf import compute_msmt_frf
+from scilpy.utils.bvec_bval_tools import extract_dwi_shell
+from scilpy.reconst.b_tensor_utils import generate_btensor_input
+
+
+def buildArgsParser():
+
+    p = argparse.ArgumentParser(description=__doc__,
+                                formatter_class=argparse.RawTextHelpFormatter)
+
+    p.add_argument('out_wm_frf',
+                   help='Path to the output WM frf file, in .txt format.')
+    p.add_argument('out_gm_frf',
+                   help='Path to the output GM frf file, in .txt format.')
+    p.add_argument('out_csf_frf',
+                   help='Path to the output CSF frf file, in .txt format.')
+
+    p.add_argument('--in_dwis', nargs='+', required=True,
+                   help='Path to the input diffusion volume for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bvals', nargs='+', required=True,
+                   help='Path to the bval file, in FSL format, for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bvecs', nargs='+', required=True,
+                   help='Path to the bvec file, in FSL format, for each '
+                        'b-tensor encoding type.')
+    p.add_argument('--in_bdeltas', nargs='+', type=float,
+                   choices=[0, 1, -0.5, 0.5], required=True,
+                   help='Value of b_delta for each b-tensor encoding type, '
+                        'in the same order as \ndwi, bval and bvec inputs.')
+
+    p.add_argument('--mask',
+                   help='Path to a binary mask. Only the data inside the mask '
+                        'will be used for\ncomputations and reconstruction. '
+                        'Useful if no tissue masks are available.')
+    p.add_argument('--mask_wm',
+                   help='Path to the input WM mask file, used to improve the'
+                        ' final WM frf mask.')
+    p.add_argument('--mask_gm',
+                   help='Path to the input GM mask file, used to improve the '
+                        'final GM frf mask.')
+    p.add_argument('--mask_csf',
+                   help='Path to the input CSF mask file, used to improve the'
+                        ' final CSF frf mask.')
+
+    p.add_argument('--fa_thr_wm',
+                   default=0.7, type=float,
+                   help='If supplied, use this threshold to select single WM '
+                        'fiber voxels from \nthe FA inside the WM mask '
+                        'defined by mask_wm. \nEach voxel above this '
+                        'threshold will be selected. [%(default)s]')
+    p.add_argument('--fa_thr_gm',
+                   default=0.2, type=float,
+                   help='If supplied, use this threshold to select GM voxels '
+                        'from the FA inside \nthe GM mask defined by mask_gm. '
+                        '\nEach voxel below this threshold will be selected. '
+                        '[%(default)s]')
+    p.add_argument('--fa_thr_csf',
+                   default=0.1, type=float,
+                   help='If supplied, use this threshold to select CSF voxels '
+                        'from the FA inside \nthe CSF mask defined by '
+                        'mask_csf. \nEach voxel below this threshold will be '
+                        'selected. [%(default)s]')
+    p.add_argument('--md_thr_gm',
+                   default=0.0007, type=float,
+                   help='If supplied, use this threshold to select GM voxels '
+                        'from the MD inside \nthe GM mask defined by mask_gm. '
+                        '\nEach voxel below this threshold will be selected. '
+                        '[%(default)s]')
+    p.add_argument('--md_thr_csf',
+                   default=0.003, type=float,
+                   help='If supplied, use this threshold to select CSF '
+                        'voxels from the MD inside \nthe CSF mask defined by '
+                        'mask_csf. \nEach voxel below this threshold will be '
+                        'selected. [%(default)s]')
+
+    p.add_argument('--min_nvox',
+                   default=100, type=int,
+                   help='Minimal number of voxels needed for each tissue masks'
+                        ' in order to \nproceed to frf estimation. '
+                        '[%(default)s]')
+    p.add_argument('--tolerance',
+                   type=int, default=20,
+                   help='The tolerated gap between the b-values to '
+                        'extract and the current b-value. [%(default)s]')
+    p.add_argument('--dti_bval_limit',
+                   type=int, default=1200,
+                   help='The highest b-value taken for the DTI model. '
+                        '[%(default)s]')
+    p.add_argument('--roi_radii',
+                   default=[20], nargs='+', type=int,
+                   help='If supplied, use those radii to select a cuboid roi '
+                        'to estimate the \nresponse functions. The roi will '
+                        'be a cuboid spanning from the middle of \nthe volume '
+                        'in each direction with the different radii. The type '
+                        'is either \nan int (e.g. --roi_radii 10) or an '
+                        'array-like (3,) (e.g. --roi_radii 20 30 10). '
+                        '[%(default)s]')
+    p.add_argument('--roi_center',
+                   metavar='tuple(3)', nargs=3, type=int,
+                   help='If supplied, use this center to span the cuboid roi '
+                        'using roi_radii. \n[center of the 3D volume] '
+                        '(e.g. --roi_center 66 79 79)')
+
+    p.add_argument('--wm_frf_mask',
+                   metavar='file', default='',
+                   help='Path to the output WM frf mask file, the voxels used '
+                        'to compute the WM frf.')
+    p.add_argument('--gm_frf_mask',
+                   metavar='file', default='',
+                   help='Path to the output GM frf mask file, the voxels used '
+                        'to compute the GM frf.')
+    p.add_argument('--csf_frf_mask',
+                   metavar='file', default='',
+                   help='Path to the output CSF frf mask file, the voxels '
+                        'used to compute the CSF frf.')
+
+    p.add_argument('--frf_table',
+                   metavar='file', default='',
+                   help='Path to the output frf table file. Saves the frf for '
+                        'each b-value, in .txt format.')
+
+    add_force_b0_arg(p)
+    add_overwrite_arg(p)
+    add_verbose_arg(p)
+
+    return p
+
+
+def main():
+    parser = buildArgsParser()
+    args = parser.parse_args()
+
+    if args.verbose:
+        logging.basicConfig(level=logging.DEBUG)
+    else:
+        logging.basicConfig(level=logging.INFO)
+
+    assert_inputs_exist(parser, [],
+                        optional=list(np.concatenate((args.in_dwis,
+                                                      args.in_bvals,
+                                                      args.in_bvecs))))
+    assert_outputs_exist(parser, args, [args.out_wm_frf, args.out_gm_frf,
+                                        args.out_csf_frf])
+
+    if not (len(args.in_dwis) == len(args.in_bvals)
+            == len(args.in_bvecs) == len(args.in_bdeltas)):
+        msg = """The number of given dwis, bvals, bvecs and bdeltas must be the
+              same. Please verify that all inputs were correctly inserted."""
+        raise ValueError(msg)
+
+    affine = extract_affine(args.in_dwis)
+
+    if len(args.roi_radii) == 1:
+        roi_radii = args.roi_radii[0]
+    elif len(args.roi_radii) == 2:
+        parser.error('--roi_radii cannot be of size (2,).')
+    else:
+        roi_radii = args.roi_radii
+    roi_center = args.roi_center
+
+    tol = args.tolerance
+    dti_lim = args.dti_bval_limit
+    force_b0_thr = args.force_b0_threshold
+
+    gtab, data, ubvals, ubdeltas = generate_btensor_input(args.in_dwis,
+                                                          args.in_bvals,
+                                                          args.in_bvecs,
+                                                          args.in_bdeltas,
+                                                          force_b0_thr,
+                                                          tol=tol)
+
+    if not np.all(ubvals <= dti_lim):
+        if np.sum(ubdeltas == 1) > 0:
+            dti_ubvals = ubvals[ubdeltas == 1]
+        elif np.sum(ubdeltas == -0.5) > 0:
+            dti_ubvals = ubvals[ubdeltas == -0.5]
+        elif np.sum(ubdeltas == args.in_bdelta_custom) > 0:
+            dti_ubvals = ubvals[ubdeltas == args.in_bdelta_custom]
+        else:
+            raise ValueError("No encoding available for DTI.")
+        vol = nib.Nifti1Image(data, affine)
+        outputs = extract_dwi_shell(vol, gtab.bvals, gtab.bvecs,
+                                    dti_ubvals[dti_ubvals <= dti_lim],
+                                    tol=1)
+        indices_dti, data_dti, bvals_dti, bvecs_dti = outputs
+        # gtab_dti = gradient_table(np.squeeze(bvals_dti), bvecs_dti,
+        #                           btens=gtab.btens[indices_dti])
+        bvals_dti = np.squeeze(bvals_dti)
+        btens_dti = gtab.btens[indices_dti]
+    else:
+        data_dti = None
+        bvals_dti = None
+        bvecs_dti = None
+        btens_dti = None
+
+    mask = None
+    if args.mask is not None:
+        mask = get_data_as_mask(nib.load(args.mask), dtype=bool)
+        if mask.shape != data.shape[:-1]:
+            raise ValueError("Mask is not the same shape as data.")
+    mask_wm = None
+    mask_gm = None
+    mask_csf = None
+    if args.mask_wm:
+        mask_wm = get_data_as_mask(nib.load(args.mask_wm), dtype=bool)
+    if args.mask_gm:
+        mask_gm = get_data_as_mask(nib.load(args.mask_gm), dtype=bool)
+    if args.mask_csf:
+        mask_csf = get_data_as_mask(nib.load(args.mask_csf), dtype=bool)
+
+    responses, frf_masks = compute_msmt_frf(data, gtab.bvals, gtab.bvecs,
+                                            btens=gtab.btens,
+                                            data_dti=data_dti,
+                                            bvals_dti=bvals_dti,
+                                            bvecs_dti=bvecs_dti,
+                                            btens_dti=btens_dti,
+                                            mask=mask, mask_wm=mask_wm,
+                                            mask_gm=mask_gm, mask_csf=mask_csf,
+                                            fa_thr_wm=args.fa_thr_wm,
+                                            fa_thr_gm=args.fa_thr_gm,
+                                            fa_thr_csf=args.fa_thr_csf,
+                                            md_thr_gm=args.md_thr_gm,
+                                            md_thr_csf=args.md_thr_csf,
+                                            min_nvox=args.min_nvox,
+                                            roi_radii=roi_radii,
+                                            roi_center=roi_center,
+                                            tol=0,
+                                            force_b0_threshold=force_b0_thr)
+
+    masks_files = [args.wm_frf_mask, args.gm_frf_mask, args.csf_frf_mask]
+    for mask, mask_file in zip(frf_masks, masks_files):
+        if mask_file:
+            nib.save(nib.Nifti1Image(mask.astype(np.uint8), vol.affine),
+                     mask_file)
+
+    frf_out = [args.out_wm_frf, args.out_gm_frf, args.out_csf_frf]
+
+    for frf, response in zip(frf_out, responses):
+        np.savetxt(frf, response)
+
+    if args.frf_table:
+        if ubvals[0] < tol:
+            bvals = ubvals[1:]
+        else:
+            bvals = ubvals
+        response_csf = responses[2]
+        response_gm = responses[1]
+        response_wm = responses[0]
+        iso_responses = np.concatenate((response_csf[:, :3],
+                                        response_gm[:, :3]), axis=1)
+        responses = np.concatenate((iso_responses, response_wm[:, :3]), axis=1)
+        frf_table = np.vstack((bvals, responses.T)).T
+        np.savetxt(args.frf_table, frf_table)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/scil_compute_msmt_fodf.py b/scripts/scil_compute_msmt_fodf.py
index c8180edd8..dd3a97298 100755
--- a/scripts/scil_compute_msmt_fodf.py
+++ b/scripts/scil_compute_msmt_fodf.py
@@ -30,7 +30,8 @@
 from scilpy.io.image import get_data_as_mask
 from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
                              assert_outputs_exist, add_force_b0_arg,
-                             add_sh_basis_args, add_processes_arg)
+                             add_sh_basis_args, add_processes_arg,
+                             add_verbose_arg)
 from scilpy.reconst.multi_processes import fit_from_model, convert_sh_basis
 from scilpy.utils.bvec_bval_tools import (check_b0_threshold, normalize_bvecs,
                                           is_normalized_bvecs)
@@ -66,6 +67,7 @@ def _build_arg_parser():
     add_sh_basis_args(p)
     add_processes_arg(p)
     add_overwrite_arg(p)
+    add_verbose_arg(p)
 
     p.add_argument(
         '--not_all', action='store_true',
@@ -137,6 +139,7 @@ def main():
     # Checking data and sh_order
     b0_thr = check_b0_threshold(
         args.force_b0_threshold, bvals.min(), bvals.min())
+
     if data.shape[-1] < (sh_order + 1) * (sh_order + 2) / 2:
         logging.warning(
             'We recommend having at least {} unique DWIs volumes, but you '
diff --git a/scripts/scil_compute_ssst_fodf.py b/scripts/scil_compute_ssst_fodf.py
index 64fc3e5cd..70bc5c0d7 100755
--- a/scripts/scil_compute_ssst_fodf.py
+++ b/scripts/scil_compute_ssst_fodf.py
@@ -18,7 +18,8 @@
 import numpy as np
 
 from scilpy.io.image import get_data_as_mask
-from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
+from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg,
+                             assert_inputs_exist, add_verbose_arg,
                              assert_outputs_exist, add_force_b0_arg,
                              add_sh_basis_args, add_processes_arg)
 from scilpy.reconst.multi_processes import fit_from_model, convert_sh_basis
@@ -53,6 +54,7 @@ def _build_arg_parser():
     add_sh_basis_args(p)
     add_processes_arg(p)
     add_overwrite_arg(p)
+    add_verbose_arg(p)
 
     return p
 
diff --git a/scripts/tests/test_compute_divide.py b/scripts/tests/test_compute_divide.py
new file mode 100644
index 000000000..b4e34b298
--- /dev/null
+++ b/scripts/tests/test_compute_divide.py
@@ -0,0 +1,119 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import os
+import tempfile
+
+from scilpy.io.fetcher import fetch_data, get_home, get_testing_files_dict
+
+fetch_data(get_testing_files_dict(), keys=['btensor_testdata.zip'])
+tmp_dir = tempfile.TemporaryDirectory()
+
+
+def test_help_option(script_runner):
+    ret = script_runner.run('scil_compute_divide.py', '--help')
+    assert ret.success
+
+
+def test_nb_btensors_check(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_plan = os.path.join(get_home(), 'btensor_testdata',
+                               'dwi_planar.nii.gz')
+    in_bval_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvals')
+    in_bvec_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvecs')
+    fa = os.path.join(get_home(), 'btensor',
+                      'fa.nii.gz')
+
+    ret = script_runner.run('scil_compute_divide.py', '--in_dwis',
+                            in_dwi_lin, '--in_bvals', in_bval_lin,
+                            '--in_bvecs', in_bvec_lin, '--in_bdeltas', '1',
+                            '--fa', fa, '--do_weight_bvals',
+                            '--do_weight_pa', '--do_multiple_s0',
+                            '--processes', '1', '-f')
+    assert (not ret.success)
+
+    ret = script_runner.run('scil_compute_divide.py', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, '--in_bvals', in_bval_lin,
+                            in_bval_plan, '--in_bvecs', in_bvec_lin,
+                            in_bvec_plan, '--in_bdeltas', '1', '1',
+                            '--fa', fa, '--do_weight_bvals',
+                            '--do_weight_pa', '--do_multiple_s0',
+                            '--processes', '1', '-f')
+    assert (not ret.success)
+
+
+def test_inputs_check(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_plan = os.path.join(get_home(), 'btensor_testdata',
+                               'dwi_planar.nii.gz')
+    in_bval_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvals')
+    in_bvec_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvecs')
+    fa = os.path.join(get_home(), 'btensor',
+                      'fa.nii.gz')
+
+    ret = script_runner.run('scil_compute_divide.py', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, '--in_bvals', in_bval_lin,
+                            '--in_bvecs', in_bvec_lin, '--in_bdeltas', '1',
+                            '--fa', fa, '--do_weight_bvals',
+                            '--do_weight_pa', '--do_multiple_s0',
+                            '--processes', '1', '-f')
+    assert (not ret.success)
+
+    ret = script_runner.run('scil_compute_divide.py', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, '--in_bvals',
+                            in_bval_lin, in_bval_plan, '--in_bvecs',
+                            in_bvec_lin, in_bvec_plan, '--in_bdeltas', '1',
+                            '-0.5', '0', '--fa', fa, '--do_weight_bvals',
+                            '--do_weight_pa', '--do_multiple_s0',
+                            '--processes', '1', '-f')
+    assert (not ret.success)
+
+
+def test_execution_processing(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_plan = os.path.join(get_home(), 'btensor_testdata',
+                               'dwi_planar.nii.gz')
+    in_bval_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvals')
+    in_bvec_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvecs')
+    in_dwi_sph = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_spherical.nii.gz')
+    in_bval_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvals')
+    in_bvec_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvecs')
+    fa = os.path.join(get_home(), 'btensor_testdata',
+                      'fa.nii.gz')
+
+    ret = script_runner.run('scil_compute_divide.py', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, in_dwi_sph,
+                            '--in_bvals', in_bval_lin, in_bval_plan,
+                            in_bval_sph, '--in_bvecs', in_bvec_lin,
+                            in_bvec_plan, in_bvec_sph, '--in_bdeltas',
+                            '1', '-0.5', '0', '--fa', fa, '--do_weight_bvals',
+                            '--do_weight_pa', '--do_multiple_s0',
+                            '--processes', '1', '-f')
+    assert (ret.success)
diff --git a/scripts/tests/test_compute_memsmt_fodf.py b/scripts/tests/test_compute_memsmt_fodf.py
new file mode 100644
index 000000000..d64d74a3e
--- /dev/null
+++ b/scripts/tests/test_compute_memsmt_fodf.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import os
+import tempfile
+
+from scilpy.io.fetcher import fetch_data, get_home, get_testing_files_dict
+
+fetch_data(get_testing_files_dict(), keys=['btensor_testdata.zip'])
+tmp_dir = tempfile.TemporaryDirectory()
+
+
+def test_help_option(script_runner):
+    ret = script_runner.run('scil_compute_memsmt_fodf.py', '--help')
+    assert ret.success
+
+
+def test_inputs_check(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_plan = os.path.join(get_home(), 'btensor_testdata',
+                               'dwi_planar.nii.gz')
+    in_bval_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvals')
+    in_bvec_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvecs')
+    in_wm_frf = os.path.join(get_home(), 'btensor_testdata',
+                             'wm_frf.txt')
+    in_gm_frf = os.path.join(get_home(), 'btensor_testdata',
+                             'gm_frf.txt')
+    in_csf_frf = os.path.join(get_home(), 'btensor_testdata',
+                              'csf_frf.txt')
+
+    ret = script_runner.run('scil_compute_memsmt_fodf.py', in_wm_frf,
+                            in_gm_frf, in_csf_frf, '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, '--in_bvals',
+                            in_bval_lin, '--in_bvecs', in_bvec_lin,
+                            '--in_bdeltas', '1',
+                            '--wm_out_fODF', 'wm_fodf.nii.gz',
+                            '--gm_out_fODF', 'gm_fodf.nii.gz',
+                            '--csf_out_fODF', 'csf_fodf.nii.gz', '--vf',
+                            'vf.nii.gz', '--sh_order', '4', '--sh_basis',
+                            'tournier07', '--processes', '1', '-f')
+    assert (not ret.success)
+
+    ret = script_runner.run('scil_compute_memsmt_fodf.py', in_wm_frf,
+                            in_gm_frf, in_csf_frf, '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, '--in_bvals',
+                            in_bval_lin, in_bval_plan, '--in_bvecs',
+                            in_bvec_lin, in_bvec_plan, '--in_bdeltas',
+                            '1', '-0.5', '0',
+                            '--wm_out_fODF', 'wm_fodf.nii.gz',
+                            '--gm_out_fODF', 'gm_fodf.nii.gz',
+                            '--csf_out_fODF', 'csf_fodf.nii.gz', '--vf',
+                            'vf.nii.gz', '--sh_order', '4', '--sh_basis',
+                            'tournier07', '--processes', '1', '-f')
+    assert (not ret.success)
+
+
+def test_execution_processing(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_sph = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_spherical.nii.gz')
+    in_bval_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvals')
+    in_bvec_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvecs')
+    in_wm_frf = os.path.join(get_home(), 'btensor_testdata',
+                             'wm_frf.txt')
+    in_gm_frf = os.path.join(get_home(), 'btensor_testdata',
+                             'gm_frf.txt')
+    in_csf_frf = os.path.join(get_home(), 'btensor_testdata',
+                              'csf_frf.txt')
+
+    ret = script_runner.run('scil_compute_memsmt_fodf.py', in_wm_frf,
+                            in_gm_frf, in_csf_frf, '--in_dwis',
+                            in_dwi_lin, in_dwi_sph, '--in_bvals',
+                            in_bval_lin, in_bval_sph,
+                            '--in_bvecs', in_bvec_lin,
+                            in_bvec_sph, '--in_bdeltas', '1', '0',
+                            '--sh_order', '8', '--processes', '8', '-f')
+    assert ret.success
diff --git a/scripts/tests/test_compute_memsmt_frf.py b/scripts/tests/test_compute_memsmt_frf.py
new file mode 100644
index 000000000..9b8fc7fe7
--- /dev/null
+++ b/scripts/tests/test_compute_memsmt_frf.py
@@ -0,0 +1,159 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import os
+import tempfile
+
+from scilpy.io.fetcher import fetch_data, get_home, get_testing_files_dict
+
+fetch_data(get_testing_files_dict(), keys=['btensor_testdata.zip'])
+tmp_dir = tempfile.TemporaryDirectory()
+
+
+def test_help_option(script_runner):
+    ret = script_runner.run('scil_compute_memsmt_frf.py', '--help')
+    assert ret.success
+
+
+def test_roi_center_shape_parameter(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_plan = os.path.join(get_home(), 'btensor_testdata',
+                               'dwi_planar.nii.gz')
+    in_bval_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvals')
+    in_bvec_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvecs')
+    in_dwi_sph = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_spherical.nii.gz')
+    in_bval_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvals')
+    in_bvec_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvecs')
+
+    ret = script_runner.run('scil_compute_memsmt_frf.py', 'wm_frf.txt',
+                            'gm_frf.txt', 'csf_frf.txt', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, in_dwi_sph, '--in_bvals',
+                            in_bval_lin, in_bval_plan, in_bval_sph,
+                            '--in_bvecs', in_bvec_lin, in_bvec_plan,
+                            in_bvec_sph, '--in_bdeltas', '1', '-0.5', '0',
+                            '--roi_center', '1', '--min_nvox', '1', '-f')
+
+    assert (not ret.success)
+
+
+def test_roi_radii_shape_parameter(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_plan = os.path.join(get_home(), 'btensor_testdata',
+                               'dwi_planar.nii.gz')
+    in_bval_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvals')
+    in_bvec_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvecs')
+    in_dwi_sph = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_spherical.nii.gz')
+    in_bval_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvals')
+    in_bvec_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvecs')
+    ret = script_runner.run('scil_compute_memsmt_frf.py', 'wm_frf.txt',
+                            'gm_frf.txt', 'csf_frf.txt', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, in_dwi_sph, '--in_bvals',
+                            in_bval_lin, in_bval_plan, in_bval_sph,
+                            '--in_bvecs', in_bvec_lin, in_bvec_plan,
+                            in_bvec_sph, '--in_bdeltas', '1', '-0.5', '0',
+                            '--roi_radii', '37', '--min_nvox', '1', '-f')
+    assert ret.success
+
+    ret = script_runner.run('scil_compute_memsmt_frf.py', 'wm_frf.txt',
+                            'gm_frf.txt', 'csf_frf.txt', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, in_dwi_sph, '--in_bvals',
+                            in_bval_lin, in_bval_plan, in_bval_sph,
+                            '--in_bvecs', in_bvec_lin, in_bvec_plan,
+                            in_bvec_sph, '--in_bdeltas', '1', '-0.5', '0',
+                            '--roi_radii', '37', '37', '37',
+                            '--min_nvox', '1', '-f')
+    assert ret.success
+
+    ret = script_runner.run('scil_compute_memsmt_frf.py', 'wm_frf.txt',
+                            'gm_frf.txt', 'csf_frf.txt', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, in_dwi_sph, '--in_bvals',
+                            in_bval_lin, in_bval_plan, in_bval_sph,
+                            '--in_bvecs', in_bvec_lin, in_bvec_plan,
+                            in_bvec_sph, '--in_bdeltas', '1', '-0.5', '0',
+                            '--roi_radii', '37', '37', '37', '37', '37',
+                            '--min_nvox', '1', '-f')
+
+    assert (not ret.success)
+
+
+def test_inputs_check(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_plan = os.path.join(get_home(), 'btensor_testdata',
+                               'dwi_planar.nii.gz')
+    in_bval_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvals')
+    in_bvec_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvecs')
+
+    ret = script_runner.run('scil_compute_memsmt_frf.py', 'wm_frf.txt',
+                            'gm_frf.txt', 'csf_frf.txt', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, '--in_bvals',
+                            in_bval_lin, '--in_bvecs', in_bvec_lin,
+                            '--in_bdeltas', '1', '--min_nvox', '1', '-f')
+    assert (not ret.success)
+
+    ret = script_runner.run('scil_compute_memsmt_frf.py', 'wm_frf.txt',
+                            'gm_frf.txt', 'csf_frf.txt', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, '--in_bvals',
+                            in_bval_lin, in_bval_plan, '--in_bvecs',
+                            in_bvec_lin, in_bvec_plan, '--in_bdeltas',
+                            '1', '-0.5', '0', '--min_nvox', '1', '-f')
+    assert (not ret.success)
+
+
+def test_execution_processing(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi_lin = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_linear.nii.gz')
+    in_bval_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvals')
+    in_bvec_lin = os.path.join(get_home(), 'btensor_testdata',
+                               'linear.bvecs')
+    in_dwi_plan = os.path.join(get_home(), 'btensor_testdata',
+                               'dwi_planar.nii.gz')
+    in_bval_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvals')
+    in_bvec_plan = os.path.join(get_home(), 'btensor_testdata',
+                                'planar.bvecs')
+    in_dwi_sph = os.path.join(get_home(), 'btensor_testdata',
+                              'dwi_spherical.nii.gz')
+    in_bval_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvals')
+    in_bvec_sph = os.path.join(get_home(), 'btensor_testdata',
+                               'spherical.bvecs')
+    ret = script_runner.run('scil_compute_memsmt_frf.py', 'wm_frf.txt',
+                            'gm_frf.txt', 'csf_frf.txt', '--in_dwis',
+                            in_dwi_lin, in_dwi_plan, in_dwi_sph, '--in_bvals',
+                            in_bval_lin, in_bval_plan, in_bval_sph,
+                            '--in_bvecs', in_bvec_lin, in_bvec_plan,
+                            in_bvec_sph, '--in_bdeltas', '1', '-0.5', '0',
+                            '--min_nvox', '1', '-f')
+    assert ret.success