Here we consider the problem of applying susceptibility distortion correction to EPI imaging data. We derive an efficient algorithm for simultaneously correcting for head motion and susceptibility distortions.

Background

Consider a collection of MRI scans of a human brain. When collected, each image is in the "scanner space", which is to say that the elements of the data arrays (voxels) correspond to locations inside an MR scanner. When comparing images for the same subject, this space is generally unhelpful, as subjects may move within the scanner, and the scanner operator may select a different location within the scanner to consider the origin

$(0,0,0)$.

Therefore images must be brought into register with one another. Generally, one uses successive approximations of transformations, such as rotations and translations, to minimize a cost function. The resulting transformation may be used to transform coordinates in one image to the equivalent coordinates in another.

Images contain an internal affine transformation that translates from indices in the data matrix into world coordinates, conventionally in millimeters (mm) right, anterior and superior of some origin. Because the world coordinate space is abstracted from the details of the data array, it is convenient to represent transformations between images as transformations between their world coordinates, not their array indices.

Head motion

During a sequence of scans, such as functional or diffusion MRI (fMRI, dMRI), subjects' brains generally move slightly due to a combination of breathing or other movements of the head and neck, as well as pulsing of the vasculature of the brain. Head motion correction is a series of registrations from each volume in a series to a single reference volume, which may be the first volume, a separately acquired reference volume, or any other convenient target.

Head motion correction is typically modeled as rigid transformations, with three translational and three rotational degrees of freedom per volume.

Susceptibility distortion

fMRI and dMRI sequences generally take advantage of magnetic susceptibility differences in tissues, such as oxygenated and deoxygenated blood, to measure the signals of interest. As a result, the sequences also suffer from distortion in macroscopic susceptibility differences, most notably in air passages. This distortion appears as stretching and compression of the image along the phase-encoding direction, and it scales with a parameter called the total readout time and the inhomogeneity of the primary magnetic field

$B_0$.

Susceptibility distortions are an interaction of the tissues, a property of the subject, and the magnetic field of the immobile scanner. Therefore, a subject in motion can produce susceptibility artifacts that are constant in direction relative to the scanner while different in direction relative to the principal axes of the brain.

It is therefore potentially problematic to correct for head motion and then susceptibility distortion, even with concatenated transformations. The susceptibility distortion must be calculated with respect to the rotation of the brain in the magnetic field

Terminology and conventions

In this section I introduce notation and terminology. The concepts, if not the specific notation should be familiar to those who have previously read a mathematical treatment of susceptibility distortion correction. The only (possibly) new concept is that of a readout field.

Spaces, indices, coordinates

Here a space is defined by an image, which contains a 3-dimensional data array and an affine matrix that map indices into world coordinates. The world coordinates and the space may sometimes be used interchangeably; the indices into the data array will always be specified when intended.

We have the following spaces:

• $\mathbf{F}$ – Fieldmap space
• $\mathbf{B}$ – BOLD reference space
• ${\mathbf{V}}^t$ – BOLD volume space for the volume at time
  $t$

The problem at hand is to resample the voxel data in

${\mathbf{V}}^t$ to

$\mathbf{B}$, taking into account distortions that are known to scale with a parameter sampled in

$\mathbf{F}$.

For a given space

$\mathbf{S}$, there is an affine matrix that transforms a discrete index

${\overrightarrow{i}}_S$ to a "world coordinate"

${\overrightarrow{c}}_S=A_S{\overrightarrow{i}}_S$.

$$\begin{array}{rl}{\overrightarrow{c}}_B& =A_B{\overrightarrow{i}}_B\\ {\overrightarrow{c}}_V^t& =A_V{\overrightarrow{i}}_V^t\\ {\overrightarrow{c}}_F& =A_F{\overrightarrow{i}}_F\end{array}$$

Without loss of generality, the BOLD volumes share a common affine, but not a common space. If a term with the

$V$ subscript omits the

$t$ superscript, the value of the term in that context does not change with

$t$.

Affine transformations

By assumption, the subject moves between volumes, so volumes have been registered to one another. In this case, the BOLD reference space

$\mathbf{B}$ acts as a common registration target. Only affine registrations are considered here, and for simplicity they may be assumed to be rigid (rotation and translation) transforms.

A transform

${\mathcal{T}}_{XY}$ is the affine transform from coordinates in space

$\mathbf{X}$ to

$\mathbf{Y}$. Assume we have invertible transforms:

$$\begin{array}{rl}{\mathcal{T}}_{BF}& :\mathbf{B}\to \mathbf{F}\\ {\mathcal{T}}_{VB}^t& :{\mathbf{V}}^{\mathbf{t}}\to \mathbf{B}\\ {\overrightarrow{c}}_F& ={\mathcal{T}}_{BF}{\overrightarrow{c}}_B\\ {\overrightarrow{c}}_B& ={\mathcal{T}}_{VB}^t{\overrightarrow{c}}_V^t\end{array}$$

To invert a transform, we may either use the

${}^{-1}$ superscript or swap the subscripts to indicate from/to.

${\mathcal{T}}_{VB}^t$ and its inverse

${\mathcal{T}}_{BV}^t$ vary by

$t$. This is in fact the only source of variation between volumes

${\mathbf{V}}_t$, and what permits

${\overrightarrow{c}}_V^t$ and

${\overrightarrow{i}}_V^t$ to be meaningfully superscripted.

$B_0$ inhomogeneity maps ("fieldmaps")

MRI acquisitions in general are sensitive to to inhomogeneities in the scanner magnetic field (

$B_0$), which are induced by inhomogeneous materials, such as the human head. A fieldmap is a measure of

$\mathrm{\Delta }{B}_0$, with the unit of Hz (radians/sec), and when we estimate this value (by any technique), we can describe the result as a discretely sampled field in the space

$\mathbf{F}$:

$${\varphi }_F({\overrightarrow{i}}_F)$$

More usefully, we will consider the function over the world coordinates, which will be interpolated over the discrete field:

$$\mathrm{\Phi }({\overrightarrow{c}}_F)$$

Readout vector

Echo-planar images are acquired with phase-encoding and frequency-encoding directions. Importantly the phase-encoding direction has a polarity, indicating which end of the data matrix is collected at an early or late phase.

BOLD volumes share a phase-encoding direction that corresponds to an (oriented) axis in the data matrix, for example j- indicates that readout will occur in the negative direction in the second dimension, which we can indicate with an orientation vector field:

$${\overrightarrow{o}}_V=\left[\begin{array}{c}0\\ -1\\ 0\end{array}\right]$$

Phase-encoding directions may be i, j or k and may have positive or negative polarity. See Algorithm for a Python implementation.

The readout field is a vector field over the voxel space of the BOLD image, and it scales with the total readout time (

${\tau }_{ro}$) of the BOLD image:

$${\overrightarrow{r}}_V={\tau }_{ro}{\overrightarrow{o}}_V$$

This readout field is an orientation vector in voxel-seconds. If scaled by a fieldmap in Hz, it would provide a displacement map in voxels along the direction of the phase encoding axis.

Per-volume distortion correction

The problem is to correct a bold volume in

${\mathbf{V}}^t$ simultaneously accounting for motion and distortion to resample in

$\mathbf{B}$.

First, let's review resampling. When resampling from a source space

$\mathbf{X}$ to a target space

$\mathbf{Y}$, we select a grid of indices

${\overrightarrow{i}}_Y$ in the target space, and find values in

$\mathbf{X}$ to place in each array element (voxel). Assuming we have an affine transform

${\mathcal{T}}_{YX}$ from

$\mathbf{Y}$ to

$\mathbf{X}$, we can map our new indices into our source data:

$$\begin{array}{rl}{\overrightarrow{i}}_X& =A_X^{-1}{\overrightarrow{c}}_X\\ & =A_X^{-1}{\mathcal{T}}_{YX}{\overrightarrow{c}}_Y\\ & =A_X^{-1}{\mathcal{T}}_{YX}{A}_Y{\overrightarrow{i}}_Y\end{array}$$

We are looking to resample a volume

${\mathbf{V}}^t$ into the the BOLD reference volume

$\mathbf{B}$, given a transform

${\mathcal{T}}_{BV}^t$:

$${\overrightarrow{i}}_V^t=A_V^{-1}{\mathcal{T}}_{BV}^t{A}_B{\overrightarrow{i}}_B$$

We would like to find an additional displacement

$\mathrm{\Delta }{\overrightarrow{i}}_V^t$ that would account for distortion induced shift. At each index

$i_V^t$, we can find a corresponding shift by mapping into the volume world space, then reference, and finally fieldmap space, and multiplying the scalar result by our readout vector:

$$\mathrm{\Delta }{\overrightarrow{i}}_V^t=\mathrm{\Phi }({\mathcal{T}}_{BF}{\mathcal{T}}_{VB}^t{A}_V{\overrightarrow{i}}_V^t){\overrightarrow{r}}_V$$

Taken over all indices, this is the fieldmap, resampled into the volume space

${\mathbf{V}}^t$, multiplied the readout vector to get a voxel shift map. Using our earlier calculations, we can decompose this into the fieldmap resampled into the reference space

$\mathbf{B}$ and the orientation vector:

$$\begin{array}{rl}\mathrm{\Delta }{\overrightarrow{i}}_V^t& =\mathrm{\Phi }({\mathcal{T}}_{BF}{\mathcal{T}}_{VB}^t{A}_V{A}_V^{-1}{\mathcal{T}}_{BV}^t{A}_B{\overrightarrow{i}}_B){\tau }_{ro}{\overrightarrow{o}}_V\\ & =\mathrm{\Phi }({\mathcal{T}}_{BF}{A}_B{\overrightarrow{i}}_B){\tau }_{ro}{\overrightarrow{o}}_V\end{array}$$

Hence, by resampling

$\mathrm{\Phi }$ once into

$\mathbf{B}$,

$${\varphi }_B({\overrightarrow{i}}_B)=\mathrm{\Phi }({\mathcal{T}}_{BF}{A}_B{\overrightarrow{i}}_B)$$

We are able to fully calculate the index-to-index mapping:

$${\overrightarrow{i}}_V^t=A_V^{-1}{\mathcal{T}}_{BV}^t{A}_B{\overrightarrow{i}}_B+{\varphi }_B({\overrightarrow{i}}_B){\tau }_{ro}{\overrightarrow{o}}_V$$

Somewhat surprisingly, the fieldmap does not need to be resampled into the individual

${\mathbf{V}}^t$ spaces in order to fully account for the effect of distortions, because it is a shift in indices based on the geometry of the target space

$\mathbf{B}$.

What is novel is that the shift applies to indices in

${\mathbf{V}}^t$ and not in

$\mathbf{B}$ where they seem to naturally occur.

Let recap the components of this final equation:

• ${\overrightarrow{i}}_B$ represents any index in the target image we are trying to compute
• $A_V^{-1}{\mathcal{T}}_{BV}^t{A}_B$ is an affine transform from the indices of an image in the reference space
  $\mathbf{B}$ to the indices of our source image in
  ${\mathbf{V}}^t$
• ${\varphi }_B$ is a precomputed fieldmap sampled in
  $\mathbf{B}$ and
  ${\varphi }_B({\overrightarrow{i}}_B)$ is a particular value in Hz
• ${\tau }_{ro}$ is the total readout time of the images in
  ${\mathbf{V}}^t$
• ${\overrightarrow{o}}_V$ is a constant orientation vector that encodes phase-encoding direction of
  ${\mathbf{V}}^t$

Note that there is no requirement for

$\mathbf{B}$ to be the BOLD reference space, nor must the transformations be affine. The critical features are that the fieldmap be resampled into the target space and used to calculate voxel shifts in the source space.

Comparison with naive composition

The final equation we derived was:

$${\overrightarrow{i}}_V^t=A_V^{-1}{\mathcal{T}}_{BV}^t{A}_B{\overrightarrow{i}}_B+{\varphi }_B({\overrightarrow{i}}_B){\tau }_{ro}{\overrightarrow{o}}_V$$

The standard way of combining head motion and distortion corrections is to concatenate transforms, which composes index mapping functions. Here, applying the voxel shift map to

${\overrightarrow{i}}_B$ and then applying the motion correction to yield

${\overrightarrow{i}}_V^t$ would yield:

$${\overrightarrow{i}}_V^t=A_V^{-1}{\mathcal{T}}_{BV}^t{A}_B({\overrightarrow{i}}_B+{\varphi }_B({\overrightarrow{i}}_B){\tau }_{ro}{\overrightarrow{o}}_V)$$

The difference in source voxel selection is

$(A_V^{-1}{\mathcal{T}}_{BV}^t{A}_B-I){\overrightarrow{o}}_V$, which intuitively measures the degree of rotation of phase-encoding axis encoded in the

${\mathcal{T}}_{BV}^t$. For small rotations or rotations about the phase-encoding axis, the difference is minimal.

Extension to nonlinear transforms

In the previous section, we assumed transformations were affine for conceptual and notational simplicity.
We will now show that invertibility is sufficient, and demonstrate that any target space may be used, and not merely an intermediate BOLD reference space.

Let

${\mathcal{T}}_{XY}$ be an invertible mapping from

$\mathbf{X}$ to

$\mathbf{Y}$:

$$\begin{array}{rl}{\overrightarrow{i}}_X& =A_X^{-1}{\overrightarrow{c}}_X\\ & =A_X^{-1}{\mathcal{T}}_{YX}({\overrightarrow{c}}_Y)\\ & =A_X^{-1}{\mathcal{T}}_{YX}(A_Y{\overrightarrow{i}}_Y)\end{array}$$

Let's consider the following common spaces:

• $\mathbf{F}$ – Fieldmap space
• $\mathbf{B}$ – BOLD reference space
• ${\mathbf{V}}^t$ – BOLD volume space for the volume at time
  $t$
• $\mathbf{P}$ – "Physical" space, as measured by a T1-weighted image, or similar
• $\mathbf{G}$ – Group space, for comparing data across subjects, such as MNI152

The following transforms (and their inverses) have been calculated:

$$\begin{array}{rl}{\mathcal{T}}_{BF}& :\mathbf{B}\to \mathbf{F}\\ {\mathcal{T}}_{BV}^t& :\mathbf{B}\to {\mathbf{V}}^t\\ {\mathcal{T}}_{PB}& :\mathbf{P}\to \mathbf{B}\\ {\mathcal{T}}_{GP}& :\mathbf{G}\to \mathbf{P}\end{array}$$

For notational convenience, the following compositions are defined, which map group coordinates into the fieldmap and individual BOLD coordinate spaces:

$$\begin{array}{rl}{\mathcal{T}}_{GV}^t& ={\mathcal{T}}_{BV}^t\circ {\mathcal{T}}_{PB}\circ {\mathcal{T}}_{GP}\\ {\mathcal{T}}_{GF}& ={\mathcal{T}}_{BF}\circ {\mathcal{T}}_{PB}\circ {\mathcal{T}}_{GP}\end{array}$$

To resample to

$\mathbf{G}$ from each volume, we need:

$$\begin{array}{rl}{\overrightarrow{i}}_V^t& =A_V^{-1}{\mathcal{T}}_{BV}^t({\mathcal{T}}_{PB}({\mathcal{T}}_{GP}(A_G{\overrightarrow{i}}_G)))\\ & =A_V^{-1}{\mathcal{T}}_{GV}^t(A_G{\overrightarrow{i}}_G)\end{array}$$

We again seek a displacement based on the fieldmap, which has been registered to the BOLD reference space, as above:

$$\mathrm{\Delta }{\overrightarrow{i}}_V^t=\mathrm{\Phi }({\mathcal{T}}_{BF}({\mathcal{T}}_{VB}^t(A_V{\overrightarrow{i}}_V^t))){\overrightarrow{r}}_V$$

Substituting

${\overrightarrow{i}}_V^t$:

$$\begin{array}{rl}\mathrm{\Delta }{\overrightarrow{i}}_V^t& =\mathrm{\Phi }({\mathcal{T}}_{BF}({\mathcal{T}}_{VB}^t(A_V{A}_V^{-1}{\mathcal{T}}_{BV}^t({\mathcal{T}}_{PB}({\mathcal{T}}_{GP}(A_G{\overrightarrow{i}}_G)))))){\tau }_{ro}{\overrightarrow{o}}_V\\ & =\mathrm{\Phi }({\mathcal{T}}_{BF}({\mathcal{T}}_{PB}({\mathcal{T}}_{GP}(A_G{\overrightarrow{i}}_G)))){\tau }_{ro}{\overrightarrow{o}}_V\\ & =\mathrm{\Phi }({\mathcal{T}}_{GF}(A_G{\overrightarrow{i}}_G)){\tau }_{ro}{\overrightarrow{o}}_V\end{array}$$

The first term resamples the fieldmap, in Hz, onto the target space:

$${\varphi }_G({\overrightarrow{i}}_G)=\mathrm{\Phi }({\mathcal{T}}_{GF}(A_G{\overrightarrow{i}}_G))$$

Hence, the final index-to-index mapping is:

$${\overrightarrow{i}}_V^t=A_V^{-1}{\mathcal{T}}_{GV}^t(A_G{\overrightarrow{i}}_G)+{\varphi }_G({\overrightarrow{i}}_G){\tau }_{ro}{\overrightarrow{o}}_V$$

We have now shown that all that is required to use this algorithm is some sequence of transforms that can map both the fieldmap and the individual BOLD volumes to a target space. The BOLD reference volume is a convenient space, but any other intermediate spaces may be used.

For example, if it is simpler to register the fieldmap to the physical space, it remains possible to find a target space. However, typically fieldmap correction is used to improve the anatomical fidelity of the EPI prior to registration to the anatomical scan. The above analysis suggests that a fieldmap registered to the anatomical space could be used when computing the cost function during coregistration.

A note on the interchangeability of spaces

The above analysis establishes that we may resample a fieldmap into
a convenient space, and we use the target space as the most clearly
convenient space. However, the exact same logic can produce any of
the following mappings:

$$\begin{array}{rl}{\overrightarrow{i}}_V^t& =A_V^{-1}{\mathcal{T}}_{GV}^t(A_G{\overrightarrow{i}}_G)+{\varphi }_G({\overrightarrow{i}}_G){\tau }_{ro}{\overrightarrow{o}}_V\\ & =A_V^{-1}{\mathcal{T}}_{GV}^t(A_G{\overrightarrow{i}}_G)+{\varphi }_P(A_P^{-1}{\mathcal{T}}_{GP}(A_G{\overrightarrow{i}}_G)){\tau }_{ro}{\overrightarrow{o}}_V\\ & =A_V^{-1}{\mathcal{T}}_{GV}^t(A_G{\overrightarrow{i}}_G)+{\varphi }_B(A_B^{-1}{\mathcal{T}}_{GB}(A_G{\overrightarrow{i}}_G)){\tau }_{ro}{\overrightarrow{o}}_V\\ & =A_V^{-1}{\mathcal{T}}_{GV}^t(A_G{\overrightarrow{i}}_G)+{\varphi }_V(A_V^{-1}{\mathcal{T}}_{GV}(A_G{\overrightarrow{i}}_G)){\tau }_{ro}{\overrightarrow{o}}_V\end{array}$$

In all cases, we are mapping indices in

$\mathbf{G}$ onto indices in

${\mathbf{V}}^t$. The only condition is that the fieldmap must be
interpolated at a high enough resolution to make these terms equivalent.

Here we make explicit the notion of a BOLD reference space with an
identical affine and index grid as each individual volume.
This is typically true of the space

$\mathbf{B}$, but rather than
require

$A_B=A_V$, we will simply use

$A_V$.

$${\overrightarrow{i}}_V^t=A_V^{-1}{\mathcal{T}}_{VV}^t{A}_V{A}_V^{-1}{\mathcal{T}}_{GV}(A_G{\overrightarrow{i}}_G)+{\varphi }_V(A_V^{-1}{\mathcal{T}}_{GV}(A_G{\overrightarrow{i}}_G)){\tau }_{ro}{\overrightarrow{o}}_V$$

Finally, let us abstract from any particular target space. We will
consider

$\mathbf{V}$,

$\mathbf{X}$ and if an intermediate space is needed
we will use

$\mathbf{W}$.

Therefore:

$$\begin{array}{rl}{\overrightarrow{i}}_V^t& =A_V^{-1}{\mathcal{T}}_{XV}^t(A_X{\overrightarrow{i}}_X)+{\varphi }_X({\overrightarrow{i}}_X){\tau }_{ro}{\overrightarrow{o}}_V\\ & =A_V^{-1}{\mathcal{T}}_{VV}^t{A}_V{A}_V^{-1}{\mathcal{T}}_{XV}(A_X{\overrightarrow{i}}_X)+{\varphi }_V(A_V^{-1}{\mathcal{T}}_{XV}(A_X{\overrightarrow{i}}_X)){\tau }_{ro}{\overrightarrow{o}}_V\\ & =A_V^{-1}{\mathcal{T}}_{VV}^t{A}_V{A}_V^{-1}{\mathcal{T}}_{WV}({\mathcal{T}}_{XW}(A_X{\overrightarrow{i}}_X))+{\varphi }_V(A_V^{-1}{\mathcal{T}}_{WV}({\mathcal{T}}_{XW}(A_X{\overrightarrow{i}}_X))){\tau }_{ro}{\overrightarrow{o}}_V\end{array}$$

Derivatives of transforms

The derivative of the affine transformation

dT_XY

Letting

$T_{GV}=A_V^{-1}{\mathcal{T}}_{GV}(A_G{\overrightarrow{i}}_G)$ be the voxel-to-voxel warp,

$$\begin{array}{rc}\frac{d{\overrightarrow{i}}_V^t}{d{\overrightarrow{i}}_G}& =T_{GV}^{\prime }({\overrightarrow{i}}_G)+{\varphi }_V^{\prime }(T_{GV}({\overrightarrow{i}}_G))T_{GV}^{\prime }({\overrightarrow{i}}_G){\tau }_{ro}{\overrightarrow{o}}_V\\ & =T_{GV}^{\prime }({\overrightarrow{i}}_G)(1+\mathrm{\nabla }{\varphi }_V(T_{GV}({\overrightarrow{i}}_G)){\tau }_{ro}{\overrightarrow{o}}_V)\\ & =T_{GV}^{\prime }({\overrightarrow{i}}_G)(1+\frac{\partial {\varphi }_V}{\partial {\overrightarrow{o}}_V}(T_{GV}({\overrightarrow{i}}_G)){\tau }_{ro})\end{array}$$

Note that

${\overrightarrow{o}}_V$ is a positive or negative unit vector
along one of the dimensions of

$\mathbf{V}$,
so the gradient simplifies to the partial derivative in the direction of

${\overrightarrow{o}}_V$.

We also neglect

${\mathcal{T}}_{VV}^t$

Resampling and intensity rescaling

To this point, we have considered the problem of mapping coordinates
from a target grid onto a source grid.
Resampling an image, however, uses mapped coordinates to interpolate
intensity values from the source image into the target grid.

We will use the symbol

${\mathcal{I}}_X^t$ to indicate the image intensity
at a location in space

$\mathbf{X}$ at time

$t$.

$$\begin{array}{rl}{\mathcal{I}}_X({\overrightarrow{i}}_X)& ={\mathcal{I}}_V^t(A_V^{-1}{\mathcal{T}}_{XV}^t(A_X{\overrightarrow{i}}_X)+{\varphi }_X({\overrightarrow{i}}_X){\tau }_{ro}{\overrightarrow{o}}_V)\end{array}$$

For any neighborhood in

$\mathbf{X}$, we would like the average intensity
of

${\mathcal{I}}_X$ to match the average intensity of the
corresponding region of

${\mathcal{I}}_V^t$.
So we need some per-voxel scaling parameter

${\mathcal{S}}_{XV}$, which
we would apply after resampling:

$${\mathcal{I}}_X({\overrightarrow{i}}_X)={\mathcal{S}}_{XV}({\overrightarrow{i}}_X){\mathcal{I}}_V^t(A_V^{-1}{\mathcal{T}}_{XV}^t(A_X{\overrightarrow{i}}_X)+{\varphi }_X({\overrightarrow{i}}_X){\tau }_{ro}{\overrightarrow{o}}_V)$$

---

$${\mathcal{I}}_X({\overrightarrow{i}}_X)={\mathcal{S}}_{XV}({\overrightarrow{i}}_X){\mathcal{I}}_V^t(A_V^{-1}{\mathcal{T}}_{XV}^t(A_X{\overrightarrow{i}}_X)+{\varphi }_V(A_V^{-1}{\mathcal{T}}_{XV}(A_X{\overrightarrow{i}}_X)){\tau }_{ro}{\overrightarrow{o}}_V)$$

---

We will show that the scale factor is the product of the Jacobian
determinants of the

${\mathcal{T}}_{XV}$ transformation and the
transformation described by the fieldmap.

$$\begin{array}{rl}{\mathcal{S}}_{XV}({\overrightarrow{i}}_X)& =|{\mathbf{J}}_{XV}({\overrightarrow{i}}_X)|\cdot |{\mathbf{J}}_{\mathrm{\Phi }}(A_V^{-1}{\mathcal{T}}_{XV}(A_X{\overrightarrow{i}}_X))|\\ & =|{\mathbf{J}}_{XV}({\overrightarrow{i}}_X)|(1+{\tau }_{ro}\frac{\partial {\varphi }_V}{\partial {\overrightarrow{o}}_V}(A_V^{-1}{\mathcal{T}}_{XV}(A_X{\overrightarrow{i}}_X)))\end{array}$$

More generally, we can accumulate Jacobians

$$|{\mathbf{J}}_{XY}({\overrightarrow{i}}_X)|=\prod _{S=\mathbf{X}}^{\mathbf{Y}}|{\mathbf{J}}_{XS}(A_S^{-1}{\mathcal{T}}_{XS}{A}_X{\overrightarrow{i}}_X)|$$

Jacobian determinant modulation

Consider a differentiable transform

$\mathcal{T}:{\mathbb{R}}^3\to {\mathbb{R}}^3$, and let

$\mathcal{D}$ be the displacement field of

$\mathcal{T}$, such that:

$$\begin{array}{rl}\mathcal{T}(x,y,z)& =\left[\begin{array}{}x\\ y\\ z\end{array}\right]+\mathcal{D}(x,y,z)\\ & =\left[\begin{array}{}x+{\mathcal{D}}_x(x,y,z)\\ y+{\mathcal{D}}_y(x,y,z)\\ z+{\mathcal{D}}_z(x,y,z)\end{array}\right]\end{array}$$

At each point we can compute the Jacobian of the transformation:

$$\begin{array}{rl}{\mathbf{J}}_{\mathcal{T}}& =\left[\begin{array}{}\frac{\partial {\mathcal{T}}_x}{\partial x}& \frac{\partial {\mathcal{T}}_x}{\partial y}& \frac{\partial {\mathcal{T}}_x}{\partial z}\\ \frac{\partial {\mathcal{T}}_y}{\partial x}& \frac{\partial {\mathcal{T}}_y}{\partial y}& \frac{\partial {\mathcal{T}}_y}{\partial z}\\ \frac{\partial {\mathcal{T}}_z}{\partial x}& \frac{\partial {\mathcal{T}}_z}{\partial y}& \frac{\partial {\mathcal{T}}_z}{\partial z}\end{array}\right]\\ & =\left[\begin{array}{}1+\frac{\partial {\mathcal{D}}_x}{\partial x}& \frac{\partial {\mathcal{D}}_x}{\partial y}& \frac{\partial {\mathcal{D}}_x}{\partial z}\\ \frac{\partial {\mathcal{D}}_y}{\partial x}& 1+\frac{\partial {\mathcal{D}}_y}{\partial y}& \frac{\partial {\mathcal{D}}_y}{\partial z}\\ \frac{\partial {\mathcal{D}}_z}{\partial x}& \frac{\partial {\mathcal{D}}_z}{\partial y}& 1+\frac{\partial {\mathcal{D}}_z}{\partial z}\end{array}\right]\end{array}$$

The Jacobian determinant describes the proportional change in volume of
a unit cube in the source coordinate space when transformed into the
target coordinate space.

We now consider the specific case of the voxel shift map in the BOLD
reference space

$\mathbf{B}$. This is a displacement field where the
coordinates are voxel indices:

$${\mathcal{D}}_{\mathrm{\Phi }}({\overrightarrow{i}}_B)={\varphi }_B({\overrightarrow{i}}_B){\tau }_{ro}{\overrightarrow{o}}_V$$

Without loss of generality, suppose that the phase-encoding direction
is i, and therefore:

$$\begin{array}{rl}{\mathcal{D}}_{\mathrm{\Phi }}(i,j,k)& ={\varphi }_B(i,j,k){\tau }_{ro}{\overrightarrow{o}}_V\\ & =\left[\begin{array}{}{\varphi }_B(i,j,k){\tau }_{ro}\\ 0\\ 0\end{array}\right]\end{array}$$

Therefore, the Jacobian determinant at each voxel index is

$$\begin{array}{rl}|{\mathbf{J}}_{\mathrm{\Phi }}|& =\left|\begin{array}{}1+{\tau }_{ro}\frac{\partial }{\partial i}{\varphi }_B& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right|\\ & =1+{\tau }_{ro}\frac{\partial {\varphi }_B}{\partial i}\end{array}$$

This is simply one plus the gradient of the voxel shift map in the phase-encoding direction. Resampling with scaling becomes:

$$\begin{array}{rl}{\mathcal{I}}_B({\overrightarrow{i}}_B)& =|{\mathbf{J}}_{\mathrm{\Phi }}({\overrightarrow{i}}_B)|{\mathcal{I}}_V({\overrightarrow{i}}_V^t)\\ & =|{\mathbf{J}}_{\mathrm{\Phi }}({\overrightarrow{i}}_B)|{\mathcal{I}}_V(A_V^{-1}{\mathcal{T}}_{BV}^t{A}_B{\overrightarrow{i}}_B+{\varphi }_B({\overrightarrow{i}}_B){\tau }_{ro}{\overrightarrow{o}}_V)\end{array}$$

Note that

${\mathcal{T}}_{BV}^t$ is affine and will thus scale all voxels equally. We are here concerned with relative changes in brightness, and can ignore any scaling effects of this term.

Note also that this Jacobian does not simplify to this trivial form as
soon as

Algorithm

Assume an interpolation function interpolate3D:

def interpolate3D(
    data: Array,
    index0: int,
    index1: int,
    index2: int,
) -> float: ...

We define a resampling function that accepts a transformation from
voxel indices to voxel indices in addition to a 4D displacement map,
such that the value at a 3D voxel index is a 3-vector.
Because the transformed indices may not be integers, we use
interpolate3D to interpolate the source data:

def resample(
    source: Array,
    target: Array,
    transformTS: Affine,
    voxel_shift_map: Array,
) -> Array:
    I, J, K = target.shape
    output = Array((I, J, K))
    for iT in range(I):
        for jT in range(J):
            for kT in range(K):
                iS, jS, kS, _ = transformTS @ [iT, jT, kT, 1]
                d_iS, d_jS, d_kS = voxel_shift_map[iT, jT, kT]
                output[i, j, k] = interpolate3D(
                    source, iS + d_iS, jS + d_jS, kS + d_kS
                )
    return output

import numpy as np

def orientation_vector(
    phase_encoding_direction: OneOf["i", "j", "k", "i-", "j-", "k-"],
) -> Vector:
    vector = np.zeros((3,))
    polarity = 1 if phase_encoding_direction[1:] != '-' else -1
    vector['ijk'.index(phase_encoding_direction[0])] = polarity
    return vector

The final algorithm to correct the volume is as follows:

def correct_volume(
    imageV: (Array, Affine),
    imageB: (Array, Affine),
    imageF: (Array, Affine),
    transformVB: Affine,
    transformBF: Affine,
    total_readout_time: float,
    phase_encoding_direction: OneOf["i", "j", "k", "i-", "j-", "k-"],
) -> Array:
    reference_fieldmap = resample(
        source=imageF.array,
        target=imageB.array,
        transformTS=inv(imageF.affine) @ transformBF @ imageB.affine,
        voxel_shift_map=Array(imageB.shape + (1,)) * 0,
    )
    orientV = orientation_vector(phase_encoding_direction)
    delta_i_V = (reference_fieldmap * total_readout_time) @ orientV
    corrected = resample(
        source=imageV.array,
        target=imageB.array,
        transformTS=inv(imageV.affine) @ inv(transformVB) @ imageB.affine,
        voxel_shift_map=delta_i_V,
    )
    return corrected