<commit_before>.. _response_fn_estimation:

Response function estimation
============

A compulsory step in spherical deconvolution is deriving the 'response
function (RF)', which is used as the kernel during the deconvolution
step. For the white matter, this is the signal expected for a voxel
containing a single, coherently-oriented fibre bundle. While some groups
prefer to define this function using some *ad-hoc* template function
(e.g. a diffusion tensor with empirical diffusivities), the MRtrix
contributors are in preference of deriving this function directly from
the image data, typically by averaging the diffusion signal from a set
of empirically-determined 'single-fibre (SF)' voxels.

The process of estimating this function from the data is however
non-trivial; there is no single unambiguous way in which this should be
done. Earlier in the beta version of MRtrix3, we provided a command
``dwi2response`` that advertised automated determination of the response
function, based on a `published
method <http://www.sciencedirect.com/science/article/pii/S1053811913008367>`__
with a few additional enhancements. Unfortunately user testing showed
that this algorithm would not produce the desired result in a number of
circumstances, and the available command-line options for altering its
behaviour were not intuitive.

As a result, we are now instead providing ``dwi2response`` as a
*script*. This was done for a few reasons. Firstly, it means that we can
provide multiple different mechanisms / algorithms for response function
estimation, all accessible within the one script, allowing users to
experiment with different approaches. Secondly, because these Python
scripts are more accessible to most users than C++ code, the algorithms
themselves can be modified, or some may even choose to try devising
their own heiristics for response function estimation. Thirdly, it
reinforces the fact that there is unfortunately *not* a black-box,
one-size-fits-all solution to this problem.

Here I will discuss some of the technical aspects of response function
estimation, and describe the mechanisms by which the currently provided
algorithms work. If however you are not interested in the nitty-gritty
of this process, feel free to scroll to the bottom of the page.

Necessary steps
---------------

Looking at the process of response function estimation in full detail,
there are four crucial steps. For each of these, I will also briefly
mention the typical process used.

1. Select those image voxels that are to be used when determining the
   response function - the 'single-fibre mask'. *Typical*: Varies.

2. Estimate the direction of the underlying fibres in each voxel.
   *Typical*: Often the diffusion tensor fit is still used for this
   purpose; though CSD itself can also be used as long as an initial
   response function estimate is available.

3. Rotate the signal measured in each single-fibre voxel in such a way
   that the estimated fibre direction coincides with the z-axis.
   *Typical*: This may be done by rotating the diffusion gradient table
   according to the estimated fibre direction; or if the diffusion
   signal is converted to spherical harmonics, then a spherical
   convolution can be used.

4. Combine these signals to produce a single response function.
   *Typical*: The ``m=0`` terms of the spherical harmonic series (which
   are rotationally symmetric about the z-axis) are simply averaged
   across single-fibre voxels.

Of these steps, the first is the one that has caused the greatest
difficulty, and is also the principle mechanism where the provided
response function estimation algorithms vary. It will therefore be the
primary focus of this document, though note that the other aspects of
this process may also change with ongoing development.

``dwi2response`` algorithms
---------------------------

``fa``
^^^^^^

In the previous version of MRtrix ('0.2'), the following heuristic was
suggested in the documentation for deriving the response function:

-  Erode a brain mask by a few voxels, to omit any voxels near the edge
   of the brain;

-  Select those voxels within the mask that have a Fractional Anisotropy
   (FA) of 0.7 or greater;

-  The ``estimate_response`` command would then be used to generate a
   response function, which would internally perform diffusion tensor
   estimation to get the fibre directions as well as the gradient
   reorientation.

Rather than this series of commands, ``dwi2response`` now provides a
similar heuristic in-built as the ``fa`` algorithm. The primary
difference is that by default, it will instead select the 300 voxels
with the highest FA (though this can be modified at the command-line).

This algorithm is provided partly for nostalgic purposes, but it also
highlights the range of possibilities for single-fbre voxel selection.
One of the problems associated with this approach (over and above the
feeling of uncleanliness from using the tensor model) is that in white
matter regions close to CSF, Gibbs ringing can make the signal in *b=0*
images erroneously low, which causes an artificial increase in FA, and
therefore such voxels get included in the single-fibre mask.

``manual``
^^^^^^^^^^

This algorithm is provided for cases where none of the available
algorithms give adequate results, for deriving multi-shell multi-tissue
response functions in cases where the voxel mask for each tissue must be
defined manually, or for anyone who may find it useful if trying to
devise their own mechanism for response function estimation. It requires
manual definition of both the single-fibre voxel mask (or just a voxel
mask for isotropic tissues); the fibre directions can also be provided
manually if necessary (otherwise a tensor fit will be used).

``msmt_5tt``
^^^^^^^^^^^^

This algorithm is intended for deriving multi-shell, multi-tissue
response functions that are compatible with the new Multi-Shell
Multi-Tissue (MSMT) CSD algorithm. The response function estimation
algorithm is identical to that described in `the
manuscript <http://linkinghub.elsevier.com/retrieve/pii/S1053-8119(14)00644-2>`__:
As long as EPI inhomogeneity field correction has been performed, and a
tissue-segmented anatomical image (prepared in the 5TT format for
:ref:`ACT <act>`) is provided with good
prior rigid-body alignment to the diffusion images, then these
high-resolution tissue segmentations can be used to identify
single-tissue voxels in the diffusion images. This algorithm is
hard-wired to provide response functions for the most typical use case
for MSMT CSD: An isotropic grey matter response, an anisotropic white
matter response, and an isotropic CSF response; the output response
functions are provided in the format expected by the ``msdwi2fod``
command. Those wishing to experiment with different multi-tissue
response function configurations will need to use the ``manual``
algorithm (which will provide a multi-shell response function if the
input DWI contains such data).

For reference, this algorithm operates as follows:

1. Resample the 5TT segmented image to diffusion image space.

2. For each of the three tissues (GM, WM, CSF), select those voxels that
   obey the following criteria:

-  The tissue partial volume fraction must be at least 0.95.

-  For GM and CSF, the FA must be no larger than 0.2.

3. For WM, use the mask derived from step 2 as the initialisation to the
   ``tournier`` algorithm, to select single-fibre voxels.

4. Derive a multi-shell response for each tissue for each of these three
   tissues. For GM and CSF, use *lmax=0* for all shells.

``tax``
^^^^^^^

This algorithm is a fairly accurate reimplementation of the approach
proposed by `Tax et
al. <http://www.sciencedirect.com/science/article/pii/S1053811913008367>`__.
The operation of the algorithm can be summarized as follows:

1. Initialise the response function using a relatively 'fat' profile,
   and the single-fibre mask using all brain voxels.

2. Perform CSD in all single-fibre voxels.

3. Exclude from the single-fibre voxel mask those voxels where the
   resulting FOD detects more than one discrete fibre population, e.g.
   using the ratio of the amplitudes of the first and second tallest
   peaks.

4. Re-calculate the response function using the updated single-fibre<commit_msg>Response function estimation concept update to reflect WM-GM-CSF ordering of dwi2response<commit_after>.. _response_fn_estimation:

Response function estimation
============

A compulsory step in spherical deconvolution is deriving the 'response
function (RF)', which is used as the kernel during the deconvolution
step. For the white matter, this is the signal expected for a voxel
containing a single, coherently-oriented fibre bundle. While some groups
prefer to define this function using some *ad-hoc* template function
(e.g. a diffusion tensor with empirical diffusivities), the MRtrix
contributors are in preference of deriving this function directly from
the image data, typically by averaging the diffusion signal from a set
of empirically-determined 'single-fibre (SF)' voxels.

The process of estimating this function from the data is however
non-trivial; there is no single unambiguous way in which this should be
done. Earlier in the beta version of MRtrix3, we provided a command
``dwi2response`` that advertised automated determination of the response
function, based on a `published
method <http://www.sciencedirect.com/science/article/pii/S1053811913008367>`__
with a few additional enhancements. Unfortunately user testing showed
that this algorithm would not produce the desired result in a number of
circumstances, and the available command-line options for altering its
behaviour were not intuitive.

As a result, we are now instead providing ``dwi2response`` as a
*script*. This was done for a few reasons. Firstly, it means that we can
provide multiple different mechanisms / algorithms for response function
estimation, all accessible within the one script, allowing users to
experiment with different approaches. Secondly, because these Python
scripts are more accessible to most users than C++ code, the algorithms
themselves can be modified, or some may even choose to try devising
their own heiristics for response function estimation. Thirdly, it
reinforces the fact that there is unfortunately *not* a black-box,
one-size-fits-all solution to this problem.

Here I will discuss some of the technical aspects of response function
estimation, and describe the mechanisms by which the currently provided
algorithms work. If however you are not interested in the nitty-gritty
of this process, feel free to scroll to the bottom of the page.

Necessary steps
---------------

Looking at the process of response function estimation in full detail,
there are four crucial steps. For each of these, I will also briefly
mention the typical process used.

1. Select those image voxels that are to be used when determining the
   response function - the 'single-fibre mask'. *Typical*: Varies.

2. Estimate the direction of the underlying fibres in each voxel.
   *Typical*: Often the diffusion tensor fit is still used for this
   purpose; though CSD itself can also be used as long as an initial
   response function estimate is available.

3. Rotate the signal measured in each single-fibre voxel in such a way
   that the estimated fibre direction coincides with the z-axis.
   *Typical*: This may be done by rotating the diffusion gradient table
   according to the estimated fibre direction; or if the diffusion
   signal is converted to spherical harmonics, then a spherical
   convolution can be used.

4. Combine these signals to produce a single response function.
   *Typical*: The ``m=0`` terms of the spherical harmonic series (which
   are rotationally symmetric about the z-axis) are simply averaged
   across single-fibre voxels.

Of these steps, the first is the one that has caused the greatest
difficulty, and is also the principle mechanism where the provided
response function estimation algorithms vary. It will therefore be the
primary focus of this document, though note that the other aspects of
this process may also change with ongoing development.

``dwi2response`` algorithms
---------------------------

``fa``
^^^^^^

In the previous version of MRtrix ('0.2'), the following heuristic was
suggested in the documentation for deriving the response function:

-  Erode a brain mask by a few voxels, to omit any voxels near the edge
   of the brain;

-  Select those voxels within the mask that have a Fractional Anisotropy
   (FA) of 0.7 or greater;

-  The ``estimate_response`` command would then be used to generate a
   response function, which would internally perform diffusion tensor
   estimation to get the fibre directions as well as the gradient
   reorientation.

Rather than this series of commands, ``dwi2response`` now provides a
similar heuristic in-built as the ``fa`` algorithm. The primary
difference is that by default, it will instead select the 300 voxels
with the highest FA (though this can be modified at the command-line).

This algorithm is provided partly for nostalgic purposes, but it also
highlights the range of possibilities for single-fbre voxel selection.
One of the problems associated with this approach (over and above the
feeling of uncleanliness from using the tensor model) is that in white
matter regions close to CSF, Gibbs ringing can make the signal in *b=0*
images erroneously low, which causes an artificial increase in FA, and
therefore such voxels get included in the single-fibre mask.

``manual``
^^^^^^^^^^

This algorithm is provided for cases where none of the available
algorithms give adequate results, for deriving multi-shell multi-tissue
response functions in cases where the voxel mask for each tissue must be
defined manually, or for anyone who may find it useful if trying to
devise their own mechanism for response function estimation. It requires
manual definition of both the single-fibre voxel mask (or just a voxel
mask for isotropic tissues); the fibre directions can also be provided
manually if necessary (otherwise a tensor fit will be used).

``msmt_5tt``
^^^^^^^^^^^^

This algorithm is intended for deriving multi-shell, multi-tissue
response functions that are compatible with the new Multi-Shell
Multi-Tissue (MSMT) CSD algorithm. The response function estimation
algorithm is identical to that described in `the
manuscript <http://linkinghub.elsevier.com/retrieve/pii/S1053-8119(14)00644-2>`__:
As long as EPI inhomogeneity field correction has been performed, and a
tissue-segmented anatomical image (prepared in the 5TT format for
:ref:`ACT <act>`) is provided with good
prior rigid-body alignment to the diffusion images, then these
high-resolution tissue segmentations can be used to identify
single-tissue voxels in the diffusion images. This algorithm is
hard-wired to provide response functions for the most typical use case
for MSMT CSD: An isotropic grey matter response, an anisotropic white
matter response, and an isotropic CSF response; the output response
functions are provided in the format expected by the ``msdwi2fod``
command. Those wishing to experiment with different multi-tissue
response function configurations will need to use the ``manual``
algorithm (which will provide a multi-shell response function if the
input DWI contains such data).

For reference, this algorithm operates as follows:

1. Resample the 5TT segmented image to diffusion image space.

2. For each of the three tissues (WM, GM, CSF), select those voxels that
   obey the following criteria:

-  The tissue partial volume fraction must be at least 0.95.

-  For GM and CSF, the FA must be no larger than 0.2.

3. For WM, use the mask derived from step 2 as the initialisation to the
   ``tournier`` algorithm, to select single-fibre voxels.

4. Derive a multi-shell response for each tissue for each of these three
   tissues. For GM and CSF, use *lmax=0* for all shells.

``tax``
^^^^^^^

This algorithm is a fairly accurate reimplementation of the approach
proposed by `Tax et
al. <http://www.sciencedirect.com/science/article/pii/S1053811913008367>`__.
The operation of the algorithm can be summarized as follows:

1. Initialise the response function using a relatively 'fat' profile,
   and the single-fibre mask using all brain voxels.

2. Perform CSD in all single-fibre voxels.

3. Exclude from the single-fibre voxel mask those voxels where the
   resulting FOD detects more than one discrete fibre population, e.g.
   using the ratio of the amplitudes of the first and second tallest
   peaks.

4. Re-calculate the response function using the updated single-fibre
{'fraction_non_alphanumeric': 0.05315455524961092, 'fraction_numerical': 0.010235843409553453, 'mean_word_length': 4.465161923454367, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4618209', 'n_tokens_mistral': 4389, 'n_tokens_neox': 3995, 'n_words': 2431}