diff --git a/.zenodo.json b/.zenodo.json
index c28789127b..6c8db536d2 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -711,6 +711,11 @@
     {
       "name": "Marina, Ana"
     },
+    {
+      "affiliation": "University of Sydney",
+      "name": "Close, Thomas",
+      "orcid": "0000-0002-4160-2134"
+    },
     {
       "name": "Davison, Andrew"
     },
diff --git a/nipype/interfaces/mrtrix3/__init__.py b/nipype/interfaces/mrtrix3/__init__.py
index 53af56ef65..9ea9850b1b 100644
--- a/nipype/interfaces/mrtrix3/__init__.py
+++ b/nipype/interfaces/mrtrix3/__init__.py
@@ -13,6 +13,8 @@
     MRConvert,
     MRResize,
     DWIExtract,
+    SHConv,
+    SH2Amp,
 )
 from .preprocess import (
     ResponseSD,
diff --git a/nipype/interfaces/mrtrix3/reconst.py b/nipype/interfaces/mrtrix3/reconst.py
index 3f21e9ad54..be0832a3f8 100644
--- a/nipype/interfaces/mrtrix3/reconst.py
+++ b/nipype/interfaces/mrtrix3/reconst.py
@@ -4,7 +4,7 @@
 
 import os.path as op
 
-from ..base import traits, TraitedSpec, File, Undefined, InputMultiObject
+from ..base import traits, TraitedSpec, File, InputMultiObject, isdefined
 from .base import MRTrix3BaseInputSpec, MRTrix3Base
 
 
@@ -50,10 +50,18 @@ class FitTensorInputSpec(MRTrix3BaseInputSpec):
             "only applies to the non-linear methods"
         ),
     )
+    predicted_signal = File(
+        argstr="-predicted_signal %s",
+        desc=(
+            "specify a file to contain the predicted signal from the tensor "
+            "fits. This can be used to calculate the residual signal"
+        ),
+    )
 
 
 class FitTensorOutputSpec(TraitedSpec):
     out_file = File(exists=True, desc="the output DTI file")
+    predicted_signal = File(desc="Predicted signal from fitted tensors")
 
 
 class FitTensor(MRTrix3Base):
@@ -81,6 +89,8 @@ class FitTensor(MRTrix3Base):
     def _list_outputs(self):
         outputs = self.output_spec().get()
         outputs["out_file"] = op.abspath(self.inputs.out_file)
+        if isdefined(self.inputs.predicted_signal):
+            outputs["predicted_signal"] = op.abspath(self.inputs.predicted_signal)
         return outputs
 
 
@@ -144,12 +154,23 @@ class EstimateFODInputSpec(MRTrix3BaseInputSpec):
             "[ az el ] pairs for the directions."
         ),
     )
+    predicted_signal = File(
+        argstr="-predicted_signal %s",
+        desc=(
+            "specify a file to contain the predicted signal from the FOD "
+            "estimates. This can be used to calculate the residual signal."
+            "Note that this is only valid if algorithm == 'msmt_csd'. "
+            "For single shell reconstructions use a combination of SHConv "
+            "and SH2Amp instead."
+        ),
+    )
 
 
 class EstimateFODOutputSpec(TraitedSpec):
     wm_odf = File(argstr="%s", desc="output WM ODF")
     gm_odf = File(argstr="%s", desc="output GM ODF")
     csf_odf = File(argstr="%s", desc="output CSF ODF")
+    predicted_signal = File(desc="output predicted signal")
 
 
 class EstimateFOD(MRTrix3Base):
@@ -183,10 +204,17 @@ class EstimateFOD(MRTrix3Base):
     def _list_outputs(self):
         outputs = self.output_spec().get()
         outputs["wm_odf"] = op.abspath(self.inputs.wm_odf)
-        if self.inputs.gm_odf != Undefined:
+        if isdefined(self.inputs.gm_odf):
             outputs["gm_odf"] = op.abspath(self.inputs.gm_odf)
-        if self.inputs.csf_odf != Undefined:
+        if isdefined(self.inputs.csf_odf):
             outputs["csf_odf"] = op.abspath(self.inputs.csf_odf)
+        if isdefined(self.inputs.predicted_signal):
+            if self.inputs.algorithm != "msmt_csd":
+                raise Exception(
+                    "'predicted_signal' option can only be used with "
+                    "the 'msmt_csd' algorithm"
+                )
+            outputs["predicted_signal"] = op.abspath(self.inputs.predicted_signal)
         return outputs
 
 
diff --git a/nipype/interfaces/mrtrix3/tests/test_auto_ConstrainedSphericalDeconvolution.py b/nipype/interfaces/mrtrix3/tests/test_auto_ConstrainedSphericalDeconvolution.py
index 1e87ebc377..c395f0d1c8 100644
--- a/nipype/interfaces/mrtrix3/tests/test_auto_ConstrainedSphericalDeconvolution.py
+++ b/nipype/interfaces/mrtrix3/tests/test_auto_ConstrainedSphericalDeconvolution.py
@@ -77,6 +77,10 @@ def test_ConstrainedSphericalDeconvolution_inputs():
             argstr="-nthreads %d",
             nohash=True,
         ),
+        predicted_signal=dict(
+            argstr="-predicted_signal %s",
+            extensions=None,
+        ),
         shell=dict(
             argstr="-shell %s",
             sep=",",
@@ -112,6 +116,9 @@ def test_ConstrainedSphericalDeconvolution_outputs():
             argstr="%s",
             extensions=None,
         ),
+        predicted_signal=dict(
+            extensions=None,
+        ),
         wm_odf=dict(
             argstr="%s",
             extensions=None,
diff --git a/nipype/interfaces/mrtrix3/tests/test_auto_EstimateFOD.py b/nipype/interfaces/mrtrix3/tests/test_auto_EstimateFOD.py
index d91d270c4d..2d15207571 100644
--- a/nipype/interfaces/mrtrix3/tests/test_auto_EstimateFOD.py
+++ b/nipype/interfaces/mrtrix3/tests/test_auto_EstimateFOD.py
@@ -80,6 +80,10 @@ def test_EstimateFOD_inputs():
             argstr="-nthreads %d",
             nohash=True,
         ),
+        predicted_signal=dict(
+            argstr="-predicted_signal %s",
+            extensions=None,
+        ),
         shell=dict(
             argstr="-shell %s",
             sep=",",
@@ -115,6 +119,9 @@ def test_EstimateFOD_outputs():
             argstr="%s",
             extensions=None,
         ),
+        predicted_signal=dict(
+            extensions=None,
+        ),
         wm_odf=dict(
             argstr="%s",
             extensions=None,
diff --git a/nipype/interfaces/mrtrix3/tests/test_auto_FitTensor.py b/nipype/interfaces/mrtrix3/tests/test_auto_FitTensor.py
index cbadd78fa6..7cf38faf8c 100644
--- a/nipype/interfaces/mrtrix3/tests/test_auto_FitTensor.py
+++ b/nipype/interfaces/mrtrix3/tests/test_auto_FitTensor.py
@@ -54,6 +54,10 @@ def test_FitTensor_inputs():
             position=-1,
             usedefault=True,
         ),
+        predicted_signal=dict(
+            argstr="-predicted_signal %s",
+            extensions=None,
+        ),
         reg_term=dict(
             argstr="-regularisation %f",
             max_ver="0.3.13",
@@ -71,6 +75,9 @@ def test_FitTensor_outputs():
         out_file=dict(
             extensions=None,
         ),
+        predicted_signal=dict(
+            extensions=None,
+        ),
     )
     outputs = FitTensor.output_spec()
 
diff --git a/nipype/interfaces/mrtrix3/tests/test_auto_SH2Amp.py b/nipype/interfaces/mrtrix3/tests/test_auto_SH2Amp.py
new file mode 100644
index 0000000000..ab75fc1f8a
--- /dev/null
+++ b/nipype/interfaces/mrtrix3/tests/test_auto_SH2Amp.py
@@ -0,0 +1,55 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from ..utils import SH2Amp
+
+
+def test_SH2Amp_inputs():
+    input_map = dict(
+        args=dict(
+            argstr="%s",
+        ),
+        directions=dict(
+            argstr="%s",
+            extensions=None,
+            mandatory=True,
+            position=-2,
+        ),
+        environ=dict(
+            nohash=True,
+            usedefault=True,
+        ),
+        in_file=dict(
+            argstr="%s",
+            extensions=None,
+            mandatory=True,
+            position=-3,
+        ),
+        nonnegative=dict(
+            argstr="-nonnegative",
+        ),
+        out_file=dict(
+            argstr="%s",
+            extensions=None,
+            name_source=["in_file"],
+            name_template="%s_amp.mif",
+            position=-1,
+            usedefault=True,
+        ),
+    )
+    inputs = SH2Amp.input_spec()
+
+    for key, metadata in list(input_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(inputs.traits()[key], metakey) == value
+
+
+def test_SH2Amp_outputs():
+    output_map = dict(
+        out_file=dict(
+            extensions=None,
+        ),
+    )
+    outputs = SH2Amp.output_spec()
+
+    for key, metadata in list(output_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(outputs.traits()[key], metakey) == value
diff --git a/nipype/interfaces/mrtrix3/tests/test_auto_SHConv.py b/nipype/interfaces/mrtrix3/tests/test_auto_SHConv.py
new file mode 100644
index 0000000000..992e6984a8
--- /dev/null
+++ b/nipype/interfaces/mrtrix3/tests/test_auto_SHConv.py
@@ -0,0 +1,52 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from ..utils import SHConv
+
+
+def test_SHConv_inputs():
+    input_map = dict(
+        args=dict(
+            argstr="%s",
+        ),
+        environ=dict(
+            nohash=True,
+            usedefault=True,
+        ),
+        in_file=dict(
+            argstr="%s",
+            extensions=None,
+            mandatory=True,
+            position=-3,
+        ),
+        out_file=dict(
+            argstr="%s",
+            extensions=None,
+            name_source=["in_file"],
+            name_template="%s_shconv.mif",
+            position=-1,
+            usedefault=True,
+        ),
+        response=dict(
+            argstr="%s",
+            extensions=None,
+            mandatory=True,
+            position=-2,
+        ),
+    )
+    inputs = SHConv.input_spec()
+
+    for key, metadata in list(input_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(inputs.traits()[key], metakey) == value
+
+
+def test_SHConv_outputs():
+    output_map = dict(
+        out_file=dict(
+            extensions=None,
+        ),
+    )
+    outputs = SHConv.output_spec()
+
+    for key, metadata in list(output_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(outputs.traits()[key], metakey) == value
diff --git a/nipype/interfaces/mrtrix3/utils.py b/nipype/interfaces/mrtrix3/utils.py
index d13b5d0ce7..f999db7001 100644
--- a/nipype/interfaces/mrtrix3/utils.py
+++ b/nipype/interfaces/mrtrix3/utils.py
@@ -765,3 +765,125 @@ class MRResize(MRTrix3Base):
     _cmd = "mrresize"
     input_spec = MRResizeInputSpec
     output_spec = MRResizeOutputSpec
+
+
+class SHConvInputSpec(CommandLineInputSpec):
+    in_file = File(
+        exists=True,
+        argstr="%s",
+        mandatory=True,
+        position=-3,
+        desc="input ODF image",
+    )
+    # General options
+    response = File(
+        exists=True,
+        mandatory=True,
+        argstr="%s",
+        position=-2,
+        desc=("The response function"),
+    )
+    out_file = File(
+        name_template="%s_shconv.mif",
+        name_source=["in_file"],
+        argstr="%s",
+        position=-1,
+        usedefault=True,
+        desc="the output spherical harmonics",
+    )
+
+
+class SHConvOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc="the output convoluted spherical harmonics file")
+
+
+class SHConv(CommandLine):
+    """
+    Convolve spherical harmonics with a tissue response function. Useful for
+    checking residuals of ODF estimates.
+
+
+    Example
+    -------
+
+    >>> import nipype.interfaces.mrtrix3 as mrt
+    >>> sh = mrt.SHConv()
+    >>> sh.inputs.in_file = 'csd.mif'
+    >>> sh.inputs.response = 'response.txt'
+    >>> sh.cmdline
+    'shconv csd.mif response.txt csd_shconv.mif'
+    >>> sh.run()                                 # doctest: +SKIP
+    """
+
+    _cmd = "shconv"
+    input_spec = SHConvInputSpec
+    output_spec = SHConvOutputSpec
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs["out_file"] = op.abspath(self.inputs.out_file)
+        return outputs
+
+
+class SH2AmpInputSpec(CommandLineInputSpec):
+    in_file = File(
+        exists=True,
+        argstr="%s",
+        mandatory=True,
+        position=-3,
+        desc="input ODF image",
+    )
+    # General options
+    directions = File(
+        exists=True,
+        mandatory=True,
+        argstr="%s",
+        position=-2,
+        desc=(
+            "The gradient directions along which to sample the spherical "
+            "harmonics MRtrix format"
+        ),
+    )
+    out_file = File(
+        name_template="%s_amp.mif",
+        name_source=["in_file"],
+        argstr="%s",
+        position=-1,
+        usedefault=True,
+        desc="the output spherical harmonics",
+    )
+    nonnegative = traits.Bool(
+        argstr="-nonnegative", desc="cap all negative amplitudes to zero"
+    )
+
+
+class SH2AmpOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc="the output convoluted spherical harmonics file")
+
+
+class SH2Amp(CommandLine):
+    """
+    Sample spherical harmonics on a set of gradient orientations.  Useful for
+    checking residuals of ODF estimates.
+
+
+    Example
+    -------
+
+    >>> import nipype.interfaces.mrtrix3 as mrt
+    >>> sh = mrt.SH2Amp()
+    >>> sh.inputs.in_file = 'sh.mif'
+    >>> sh.inputs.directions = 'grads.txt'
+    >>> sh.cmdline
+    'sh2amp sh.mif grads.txt sh_amp.mif'
+    >>> sh.run()                                 # doctest: +SKIP
+    """
+
+    _cmd = "sh2amp"
+    input_spec = SH2AmpInputSpec
+    output_spec = SH2AmpOutputSpec
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs["out_file"] = op.abspath(self.inputs.out_file)
+        return outputs
diff --git a/nipype/testing/data/grads.txt b/nipype/testing/data/grads.txt
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/nipype/testing/data/sh.mif b/nipype/testing/data/sh.mif
new file mode 100644
index 0000000000..e69de29bb2