idvf: Iterative inversion of deformation vector field with adaptive bi-residual feedback control

Contents

• Algorithm overview
  • Two inverse consistency residuals
  • DVF inversion iteration with bi-residual feedback control
• References
• Getting started
  • Installation
  • Testing
• Software description
  • DVF inversion
  • Preprocessing & DVF characterization
  • Post-evaluation
  • GPU utilization
  • Eigenvalue calculations
• System environment
• License and community guidelines
• Contributors

Algorithm overview

This repository contains a set of MATLAB functions for fast and accurate inversion of 2D/3D deformation vector fields (DVFs), which are also known as flow or dense motion fields. The inversion algorithm framework and a unified analysis are described in References [1,2]; the code itself is published as Reference [3]. To our knowledge, this is the first DVF inversion framework that is supported by theoretical guarantees under mild and verifiable conditions. We give a brief overview of the DVF inversion algorithm which underlies the repository functions.

We assume that we are provided with a deformation vector field (DVF) denoted by u. The DVF describes a (non-linear) mapping from a reference image I_r onto a study image I_s via point-wise displacement, i.e.,

I_r(y) = I_s(y + u(y)),

at all locations y in the reference domain. We aim at computing the inverse DVF v such that

I_s(x) = I_r(x + v(x)),

at all x in the study domain.

Two inverse consistency residuals

If two DVFs u and v are inverse of each other, they must meet the inverse consistency (IC) condition. We define the study IC residual,

s(x) = v(x) + u(x + v(x)),

and the reference IC residual,

r(y) = u(y) + v(y + u(y)).

The DVFs u and v are inverse-consistent if both residuals are zero. The residuals are related to the inversion error, e(x) = v(x) − v*(x), where v* is the true inverse DVF. The IC residuals are computationally available, while the inversion error is unknown. We use the IC residuals as feedback in iterative DVF inversion, and exercise adaptive control over the feedback for global convergence and local acceleration.

DVF inversion iteration with bi-residual feedback control

Due to the non-linear nature of the DVF mapping, one resorts to iterative procedures for DVF inversion. Our procedure has two stages: preprocessing and adaptive inversion iteration.

In the preprocessing stage, we compute the DVF Jacobians and their eigenvalues over the spatial domain. We use this spectral information to determine data-adaptive parameters for feedback control.

The inversion iteration proceeds in two phases. During phase one, at step k+1, we update the current inverse estimate v_k using the study IC residual s_k as feedback:

v_k+1(x) = v_k(x) − (1 − μ_k(x)) s_k(x),

where μ is the feedback control parameter. With the adaptive control schemes provided in this repository, the iteration is guaranteed to converge globally, under certain mild conditions [1].

Once the errors are made sufficiently small, the iteration is switched to phase two:

v_k+1(x) = v_k(x) − r_k(x + v_k(x)),

where the reference IC residual r_k at displaced locations is used as feedback. Inversion errors can be estimated using the IC residuals and the DVF spectral information. The local convergence rate is quadratic during phase two. We may refer to this iteration as implicit Newton iteration. However, phase-two iteration steps do not entail explicit formation and inversion of the Newton matrix; the explicit Newton step suffers from several numerical problems.

Phase transition and integration in the inversion iteration is facilitated by a multi-resolution scheme, elaborated in Reference [2].

References

[1] A. Dubey*, A.S. Iliopoulos*, X. Sun, F.F. Yin, and L. Ren, "Iterative inversion of deformation vector fields with feedback control", Medical Physics, 45(7):3147-3160, 2018. (arXiv)

[2] A. Dubey, Symmetric completion of deformable registration via bi-residual inversion, PhD thesis, Duke University, Durham, NC, USA, 2018.

[3] A.S. Iliopoulos, A. Dubey, and X. Sun, "idvf: Iterative inversion of deformation vector field with adaptive bi-residual feedback control", Journal of Open Source Software, 4(35):1076, 2019.

Getting started

Installation

To use idvf, simply add its top-level directory to the MATLAB path. All functions are organized in packages.

Testing

To test idvf, run test_idvf (which simply runs the included demo scripts demo_inversion_2d, demo_inversion_3d_z0, and demo_inversion_3d_zsin). Each demo performs inversion of a synthetic forward DVF using 8 different feedback control settings, and produces the following visualizations:

• Deformation of a synthetic grid-like image by the forward DVF.
• Spectral measure maps of the forward DVF (preprocessing).
• IC residual magnitudes (post-evaluation) with each control setting: percentiles at each step of the iteration, and spatial maps at the end of the iteration.
• Image-space difference maps in reference image recovery with each control setting.

The adaptive control schemes are guaranteed to converge; that is, the corresponding IC residuals tend towards zero, at a rate that depends on the control setting. The spectral measure maps provide useful information regarding the invertibility and controllability of the forward DVF over its spatial domain.

Software description

DVF inversion

The main function is dvf.inversion. It takes a forward DVF u and returns the inverse DVF v. A 3D DVF is represented by a 4D array of size N_x x N_y x N_z x 3, in single or double precision. Specifically, U(i,j,k,1:3) is the 3D forward displacement vector u(x) at voxel x = (i,j,k). Two-dimensional DVFs are represented analogously.

We provide a number of options in function dvf.inversion, which are input as name-value pairs. These options pertain to the choice of feedback control schemes and iteration parameters. The function documentation details the supported options.

The memory requirement of dvf.inversion is approximately 5 times the amount of memory for the input DVF data array.

Preprocessing & DVF characterization

Preprocessing functions are invoked automatically by the main function when adaptive feedback control is chosen. Functions dvf.jacobian and dvf.eigJacobian compute the DVF Jacobians and their eigenvalues, respectively, at all pixels/voxels. Function dvf.feedbackControlVal takes the eigenvalues and returns adaptive feedback control parameter values per the chosen scheme.

Function dvf.ntdcMeasures calculates three characteristic spectral measures of a DVF. Specifically, it returns the algebraic control index map, the non-translational component spectral radius map, and the determinant map. See Reference [1] for a description of the spectral measures and their relation to the inversion iteration.

Post-evaluation

The two IC residuals enable post-evaluation of inverse DVF estimates in the absence of the true inverse DVF. Function dvf.inversion returns IC residual magnitude percentiles at each iteration step as optional output. Function dvf.icResidual returns an IC residual field given a DVF pair (reference or study IC residual, depending on the order of the input DVFs).

Certain regions of the study domain are invalid for IC residual calculations. These regions cover points that are either mapped outside the domain by the forward DVF, or lie outside the deformed reference domain. Function dvf.maskDomain calculates the valid domain.

GPU utilization

All dvf.* functions are GPU-enabled, supported by the MATLAB Parallel Computing Toolbox. If a compatible GPU is available, GPU computation is invoked by casting the input data arrays into gpuArray objects; for example:

V = dvf.inversion(gpuArray(U));

The output data arrays are returned as gpuArray objects.

Eigenvalue calculations

Eigenvalue calculations in dvf.eigJacobian are currently slow with large DVFs. They are implemented as a loop that calls the MATLAB eig function for each Jacobian over the spatial domain. This approach is inefficient but has the benefit of numerical stability. In our experience with respiratory DVFs over thoracic and abdominal CT domains, we have not encountered any major issues with efficient but numerically unstable alternatives such as closed-form cubic root calculation. Nevertheless, we do not include such alternatives in idvf, in the interest of stability and simplicity over efficiency.

System environment

The idvf code was developed and tested on MATLAB R2018b. Dependencies to MATLAB toolboxes include the following functions:

• imresize, imresize3, imdilate, imerode, and imclose (Image Processing Toolbox, tested with version 10.3);
• prctile (Statistics and Machine Learning Toolbox, tested with version 11.4); and
• gpuArray (Parallel Computing Toolbox, tested with version 6.13).

GPU computations were tested on an NVIDIA Quantum TXR113-1000R and CUDA 10.0 drivers.

License and community guidelines

The idvf code is licensed under the GNU general public license v3.0. If you wish to contribute to idvf or report any bugs/issues, please see our contribution guidelines and code of conduct.

Contributors

• Design, development, testing:
  Alexandros-Stavros Iliopoulos, Abhishek Dubey, and Xiaobai Sun
  Department of Computer Science, Duke University

• Consulting on medical CT image analysis and applications:
  Lei Ren and Fang-Fang Yin
  Department of Radiation Oncology, Duke University Medical School

• Review, suggestions, bug reports:
  @Kevin-Mattheus-Moerman