2.1

Learned Iterative CT Reconstruction

The forward model of a CT acquisition can be written as

where x ∈ ℝ^M is the volume, p ∈ ℝ^N is the projection data, A ∈ ℝ^N×M is the forward operator defined by the imaging geometry (cone-beam in this case), and ϵ ∈ ℝ^N is additive noise. Typically, recovering the volume x from the measured data p is an ill-posed inverse problem meaning that A is not square and there are multiple solutions for x which are consistent with the measured data p. Hence, the reconstruction is formulated as a regularized optimization problem

Here, D : ℝ^M × ℝ^N → ℝ is a function which measures the consistency of volume x and measured data p. We choose the data consistency term in a least squares sense given as . R : ℝ^M → ℝ is a regularizer which helps finding a favorable solution x* and is weighted against the data consistency term by a scalar μ ∈ ℝ. Using a simple gradient descent optimization scheme, the resulting update formula is

where A^T ∈ ℝ^M×N is the adjoint operator of A, λ ∈ ℝ is a sufficiently small step size, and n is the iteration index. To learn a flexible regularizer from data, we replace its gradient by a network N_θ : ℝ^M → ℝ^M with free parameters θ which allows to fit the regularizing component directly from data. The final update formula is given as

If x represents a 2D image, this equation could inspire an unrolled network architecture (Figure 1a) and the parameters θ can be trained from pairs of input projection data and ground truth reconstruction for a fixed number of unrolled iterations.¹ This is infeasible for 3D cone-beam data due to extremely high GPU memory requirements.

Figure 1:

Comparison of an unrolled network (a) and a neural ODE-based version (b). While the unrolled version explicitly repeats the same network block for a fixed number of steps, the ODE solver receives only one network block parameterizing the temporal derivative of the volume. The solver computes the numerical integration internally. All trainable parts of our architecture are highlighted in yellow, the fixed parts are gray.

2.2

Neural Ordinary Differential Equations

Equation 4 has a residual form of the type

with f_θ(x, p) = –λ(A^T(Ax – p)+μN_θ(x)). The function f_θ describes how the volume xⁿ changes incrementally. Chen et al.¹² proposed to regard such residual neural network architectures as the numerical integration of some underlying continuous ordinary differential equation. Here, the continuous differential equation would be . Following that idea, a full unrolled iterative reconstruction is similar to the solution of an initial value problem of the given ODE starting from an initial condition x⁰ integrated until some end time T

There exist a number of different numerical solvers to approximate solutions of such initial value problems. In this work, we use a fixed step-size Runge-Kutta solver of order 4 which integrates Eq. 6 by dividing the interval [0,., T] into a fixed number of steps S to solve the integral numerically. This highlights the analogy to residual networks of depth S.

To be able to combine this ODE-based problem formulation with a trainable network architecture, we need to compute a loss that is based on the output of the ODE solver and use its gradient to update the weights θ contained in f_θ. As demonstrated by Chen et al.,¹² this gradient can be obtained without backpropagating through the internals of the solver. Instead, the solver is regarded as a black box and the gradient with respect to θ is computed by solving another ODE backward in time (adjoint sensitivity method). This allows to obtain gradients with a memory cost that is independent of the number of steps S taken by the solver. Coming back to the analogy with residual networks, we can unroll the reconstruction problem in many steps using neural ODEs without further increasing the memory cost.

2.3

Network Architecture

In the neural ODE setting, a neural network defines the dynamics of the system, i.e., its temporal derivative. Following the classical problem formulation in Eq. 4, we design this network using two branches: (1) A data consistency branch with no trainable parameters incorporating the system’s forward and backward model as known operator and (2) a regularization branch which is trained from data. Figure 1b illustrates the proposed network architecture. The data consistency branch implements the term A^T(Ax – p). Operators A and A^T are the CT forward and backprojection under the correct cone-beam geometry, respectively. We use the differentiable version of these operators described by Syben et al.¹⁴ which computes analytical gradients and allows for a direct embedding of these operators in neural networks. For the regularization branch, we use a standard 3D U-Net¹⁵ with depth 4 and 8 feature maps on the first level which are doubled in each stage, ReLU activation function, and instance normalization. In total, this leads to a network with 255 000 free parameters. We further introduce an additional single trainable parameter γ (referred to as data consistency weight) which is multiplied to the output of the data consistency branch before adding the output of both branches together in order to enable a data-optimal weighting between the data consistency and regularization component. This network together with the initialization x⁰ is passed to the ODE solver.

3.1

Data

The data set consists of 42 walnuts scanned under cone-beam CT geometry (cone angle: 40°).¹⁶ It contains the raw projection data and a corresponding ground truth reconstruction. The projection images are acquired on three circular trajectories on different heights along the walnuts’ long axes with a full rotation of 360° divided into 1200 angular steps each. The ground truth reconstruction is computed iteratively from the full set of acquired projections. As input to our algorithm, we use projection data from only the central one of the three trajectories and downsample the projections by a factor of 10 in angular direction and 2 in the spatial directions. This results in 120 projection images of size 384 × 486 pixels with an angular increment of 3° covering a full rotation of 360°. Hence, the algorithm has much less projection data available than has been used for computing the ground truth reconstruction which serves as learning target and is downsampled by a factor of 2 resulting in volumes of size 251³. We use walnut number 1 for validation and walnut number 2 for testing. The rest is used for training. The used data, its preprocessing and the train-test-splits are identical to Rudzusika et al.¹⁰

3.2

Training Details

We use the neural ODE solver provided by Chen et al.¹² Integration of Eq. 6 is performed from 0 to T = 1 with a fixed step size of 0.05. This results in 20 steps and 80 evaluations of the network per forward and backward pass as the fourth order Runge-Kutta solver takes four evaluations per step. The initial volume x⁰ is a Feldmann-Davis-Kress (FDK) reconstruction of the projection data. The last layer of the U-Net is initialized with zeros such that the regularizer has no influence upon initialization and the data consistency weight is initialized with γ = 0.01. The parameters are optimized by an Adam optimizer with learning rate 1 × 10^–4 for the U-Net and 1 × 10^–2 for the data consistency weight. Training is performed with batch size 1 for 80 epochs and an L1-Loss evaluated only inside the cylindrical scan field of view (FOV) captured in each projection. The model weights corresponding to the lowest validation loss during training are selected for further evaluation.

3.3

Reference Methods

The reconstruction of the proposed method is compared to (1) an FDK reconstruction, (2) a SIRT reconstruction with non-negativity constraint (500 iterations)¹⁷ and (3) a total variation (TV) regularized reconstruction using gradient updates (300 iterations).¹⁴