<style>
.nobr { white-space:nowrap; }
</style>

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 $\vec i_S$ to a "world coordinate" $\vec c_S = A_S\vec i_S$.

$$\begin{align}
\vec c_B &= A_B \vec i_B \\
\vec c_V^t &= A_V \vec i_V^t \\
\vec c_F &= A_F \vec i_F \\
\end{align}$$

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{align}
\mathcal T_{BF} &: \mathbf{B} \rightarrow \mathbf{F} \\
\mathcal T_{VB}^t &: \mathbf{V^t} \rightarrow \mathbf{B} \\
\vec c_F &= \mathcal T_{BF} \vec c_B \\
\vec c_B &= \mathcal T_{VB}^t \vec c_V^t
\end{align}$$

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 <span class="nobr">$\mathbf V_t$</span>, and what permits $\vec c_V^t$ and $\vec 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 $\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$:

$$\phi_F(\vec i_F)$$

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

$$\Phi(\vec 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:

$$\vec o_V = \begin{bmatrix}0\\-1\\0\end{bmatrix}$$

Phase-encoding directions may be `i`, `j` or `k` and may have positive or negative polarity. See [Algorithm](#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:

$$\vec r_V = \tau_{ro} \vec 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 $\vec 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{align}
\vec i_X &= A_X^{-1} \vec c_X \\
&= A_X^{-1} \mathcal T_{YX} \vec c_Y \\
&= A_X^{-1} \mathcal T_{YX} A_Y \vec i_Y
\end{align}$$

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$:

$$\vec i_V^t = A_V^{-1} \mathcal T_{BV}^t A_B \vec i_B$$

We would like to find an additional displacement $\Delta \vec 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:

$$
\Delta \vec i_V^t = \Phi(\mathcal{T}_{BF}\mathcal T_{VB}^tA_V\vec i_V^t)\vec 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{align}
\Delta \vec i_V^t &= \Phi(\mathcal{T}_{BF}\mathcal T_{VB}^tA_VA_V^{-1} \mathcal T_{BV}^t A_B \vec i_B)\tau_{ro}\vec o_V \\
&= \Phi(\mathcal{T}_{BF} A_B \vec i_B)\tau_{ro}\vec o_V \\
\end{align}$$

Hence, by resampling $\Phi$ once into $\mathbf B$,

$$\phi_B(\vec i_B) = \Phi(\mathcal{T}_{BF}A_B\vec i_B)$$

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

$$
\vec i_V^t = A_V^{-1} \mathcal T_{BV}^t A_B \vec i_B + \phi_B(\vec i_B)\tau_{ro}\vec 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:

* $\vec 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$
* $\phi_B$ is a precomputed fieldmap sampled in $\mathbf B$ and $\phi_B(\vec i_B)$ is a particular value in Hz
* $\tau_{ro}$ is the total readout time of the images in $\mathbf V^t$
* $\vec 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:

$$
\vec i_V^t = A_V^{-1} \mathcal T_{BV}^t A_B \vec i_B + \phi_B(\vec i_B)\tau_{ro}\vec 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 $\vec i_B$ and then applying the motion correction to yield $\vec i_V^t$ would yield:

$$
\vec i_V^t = A_V^{-1} \mathcal T_{BV}^t A_B (\vec i_B + \phi_B(\vec i_B)\tau_{ro}\vec o_V)
$$

The difference in source voxel selection is $(A_V^{-1} \mathcal T_{BV}^t A_B - I)\vec 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{align}
\vec i_X &= A_X^{-1} \vec c_X \\
&= A_X^{-1} \mathcal T_{YX}(\vec c_Y) \\
&= A_X^{-1} \mathcal T_{YX}(A_Y \vec i_Y)
\end{align}$$

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{align}
\mathcal T_{BF} &: \mathbf{B} \rightarrow \mathbf{F} \\
\mathcal T_{BV}^t &: \mathbf{B} \rightarrow \mathbf{V}^t \\
\mathcal T_{PB} &: \mathbf{P} \rightarrow \mathbf{B} \\
\mathcal T_{GP} &: \mathbf{G} \rightarrow \mathbf{P} \\
\end{align}$$

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

$$\begin{align}
\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{align}$$

To resample to $\mathbf G$ from each volume, we need:

$$\begin{align}
\vec i_V^t &= A_V^{-1} \mathcal T_{BV}^t(\mathcal T_{PB}(\mathcal T_{GP}
(A_G \vec i_G))) \\
&= A_V^{-1} \mathcal T_{GV}^t(A_G \vec i_G) \\
\end{align}$$

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

$$
\Delta \vec i_V^t = \Phi(\mathcal{T}_{BF}(\mathcal T_{VB}^t(A_V\vec i_V^t)))\vec r_V \\
$$

Substituting $\vec i_V^t$:

$$\begin{align}
\Delta \vec i_V^t &= \Phi(\mathcal{T}_{BF}(\mathcal T_{VB}^t(A_VA_V^{-1} \mathcal T_{BV}^t(\mathcal T_{PB}(\mathcal T_{GP}(A_G \vec i_G))))))\tau_{ro}\vec o_V \\
&= \Phi(\mathcal{T}_{BF}(\mathcal T_{PB}(\mathcal T_{GP}(A_G \vec i_G))))\tau_{ro}\vec o_V \\
&= \Phi(\mathcal{T}_{GF}(A_G \vec i_G))\tau_{ro}\vec o_V \\
\end{align}$$

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

$$\phi_G(\vec i_G) = \Phi(\mathcal{T}_{GF}(A_G \vec i_G))$$

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

$$
\vec i_V^t = A_V^{-1} \mathcal T_{GV}^t(A_G \vec i_G) + \phi_G(\vec i_G)\tau_{ro}\vec 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{align}{}
\vec i_V^t &= A_V^{-1} \mathcal T_{GV}^t(A_G \vec i_G) + \phi_G(\vec i_G)\tau_{ro}\vec o_V \\
&= A_V^{-1} \mathcal T_{GV}^t(A_G \vec i_G) + \phi_P(A_P^{-1} \mathcal T_{GP}(A_G \vec i_G))\tau_{ro}\vec o_V \\
&= A_V^{-1} \mathcal T_{GV}^t(A_G \vec i_G) + \phi_B(A_B^{-1} \mathcal T_{GB}(A_G \vec i_G))\tau_{ro}\vec o_V \\
&= A_V^{-1} \mathcal T_{GV}^t(A_G \vec i_G) + \phi_V(A_V^{-1} \mathcal T_{GV}(A_G \vec i_G))\tau_{ro}\vec o_V \\
\end{align}
$$

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$.

$$
\vec i_V^t = A_V^{-1} \mathcal T_{VV}^tA_V A_V^{-1} \mathcal T_{GV}(A_G \vec i_G) + \phi_V(A_V^{-1} \mathcal T_{GV}(A_G \vec i_G))\tau_{ro}\vec 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{align}{}
\vec i_V^t &= A_V^{-1} \mathcal T_{XV}^t(A_X \vec i_X) + \phi_X(\vec i_X)\tau_{ro}\vec o_V \\
&= A_V^{-1} \mathcal T_{VV}^tA_V A_V^{-1} \mathcal T_{XV}(A_X \vec i_X) + \phi_V(A_V^{-1} \mathcal T_{XV}(A_X \vec i_X))\tau_{ro}\vec o_V \\
&= A_V^{-1} \mathcal T_{VV}^tA_V A_V^{-1} \mathcal T_{WV}(\mathcal T_{XW}(A_X \vec i_X)) + \phi_V(A_V^{-1} \mathcal T_{WV}(\mathcal T_{XW}(A_X \vec i_X)))\tau_{ro}\vec o_V
\end{align}
$$

## Derivatives of transforms

The derivative of the affine transformation

dT_XY

Letting $T_{GV} = A_V^{-1} \mathcal T_{GV}(A_G \vec i_G)$ be the voxel-to-voxel warp,

$$
\begin{eqnarray}
\frac{d\vec i_V^t}{d\vec i_G} &= T_{GV}'(\vec i_G) + \phi_V'(T_{GV}(\vec i_G))T'_{GV}(\vec i_G)\tau_{ro}\vec o_V \\
&= T_{GV}'(\vec i_G)\left( 1 + \nabla \phi_V(T_{GV}(\vec i_G))\tau_{ro}\vec o_V\right) \\
&= T_{GV}'(\vec i_G)\left( 1 + \frac{\partial\phi_V}{\partial\vec o_V}(T_{GV}(\vec i_G))\tau_{ro}\right) \\
\end{eqnarray}
$$

Note that $\vec 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 $\vec 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{align}{}
\mathcal I_X(\vec i_X) &= \mathcal I_V^t(A_V^{-1}\mathcal T_{XV}^t(A_X \vec i_X) + \phi_X(\vec i_X)\tau_{ro}\vec o_V)
\end{align}
$$

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(\vec i_X) = \mathcal S_{XV}(\vec i_X)\mathcal I_V^t(A_V^{-1}\mathcal T_{XV}^t(A_X \vec i_X) + \phi_X(\vec i_X)\tau_{ro}\vec o_V)
$$

---

$$
\mathcal I_X(\vec i_X) = \mathcal S_{XV}(\vec i_X)\mathcal I_V^t(A_V^{-1}\mathcal T_{XV}^t(A_X \vec i_X) + \phi_V(A_V^{-1}\mathcal T_{XV}(A_X \vec i_X))\tau_{ro}\vec 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{align}{}
\mathcal S_{XV}(\vec i_X) &= |\mathbf J_{XV}(\vec i_X)|\cdot|\mathbf J_{\Phi}(A_V^{-1}\mathcal T_{XV}(A_X\vec i_X))| \\
& = |\mathbf J_{XV}(\vec i_X)|\left(1 + \tau_{ro}\frac{\partial\phi_V}{\partial\vec o_V}(A_V^{-1}\mathcal T_{XV}(A_X\vec i_X))\right)
\end{align}
$$

More generally, we can accumulate Jacobians

$$
|\mathbf J_{XY}(\vec i_X)| = \prod_{S=\mathbf X}^\mathbf Y|\mathbf J_{XS}(A_S^{-1}\mathcal T_{XS}A_X\vec 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{align}{}
\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{align}
$$

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

$$
\begin{align}{}
\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{align}
$$

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_\Phi(\vec i_B) = \phi_B(\vec i_B)\tau_{ro}\vec o_V $$

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

$$
\begin{align}{}
\mathcal D_\Phi(i, j, k) &= \phi_B(i, j, k)\tau_{ro}\vec o_V \\
&= \left[\begin{array}{}\phi_B(i, j, k)\tau_{ro}\\0\\0\end{array}\right]
\end{align}
$$

Therefore, the Jacobian determinant at each voxel index is

$$
\begin{align}{}
|\mathbf J_\Phi|
&= \left|
\begin{array}{}
1 + \tau_{ro}\frac{\partial}{\partial i}\phi_B &
0 &
0 \\
0 &
1 &
0 \\
0 &
0 &
1 \\
\end{array}
\right| \\
&= 1 + \tau_{ro}\frac{\partial\phi_B}{\partial i}
\end{align}
$$

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

$$
\begin{align}{}
\mathcal I_B(\vec i_B) &= |\mathbf J_\Phi(\vec i_B)|\mathcal I_V(\vec i_V^t) \\
&= |\mathbf J_\Phi(\vec i_B)|\mathcal I_V(A_V^{-1} \mathcal T_{BV}^tA_B \vec i_B + \phi_B(\vec i_B)\tau_{ro}\vec o_V)
\end{align}
$$

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