Deformable image registration for adaptive radiotherapy with guaranteed local rigidity constraints

Background Deformable image registration (DIR) is a key component in many radiotherapy applications. However, often resulting deformations are not satisfying, since varying deformation properties of different anatomical regions are not considered. To improve the plausibility of DIR in adaptive radiotherapy in the male pelvic area, this work integrates a local rigidity deformation model into a DIR algorithm. Methods A DIR framework is extended by constraints, enforcing locally rigid deformation behavior for arbitrary delineated structures. The approach restricts those structures to rigid deformations, while surrounding tissue is still allowed to deform elastically. The algorithm is tested on ten CT/CBCT male pelvis datasets with active rigidity constraints on bones and prostate and compared to the Varian SmartAdapt deformable registration (VSA) on delineations of bladder, prostate and bones. Results The approach with no rigid structures (REG0) obtains an average dice similarity coefficient (DSC) of 0.87 ± 0.06 and a Hausdorff-Distance (HD) of 8.74 ± 5.95 mm. The new approach with rigid bones (REG1) yields a DSC of 0.87 ± 0.07, HD 8.91 ± 5.89 mm. Rigid deformation of bones and prostate (REG2) obtains 0.87 ± 0.06, HD 8.73 ± 6.01 mm, while VSA yields a DSC of 0.86 ± 0.07, HD 10.22 ± 6.62 mm. No deformation grid foldings are observed for REG0 and REG1 in 7 of 10 cases; for REG2 in 8 of 10 cases, with no grid foldings in prostate, an average of 0.08 % in bladder (REG2: no foldings) and 0.01 % inside the body contour. VSA exhibits grid foldings in each case, with an average percentage of 1.81 % for prostate, 1.74 % for bladder and 0.12 % for the body contour. While REG1 and REG2 keep bones rigid, elastic bone deformations are observed with REG0 and VSA. An average runtime of 26.2 s was achieved with REG1; 31.1 s with REG2, compared to 10.5 s with REG0 and 10.7 s with VMS. Conclusions With accuracy in the range of VSA, the new approach with constraints delivers physically more plausible deformations in the pelvic area with guaranteed rigidity of arbitrary structures. Although the algorithm uses an advanced deformation model, clinically feasible runtimes are achieved. Electronic supplementary material The online version of this article (doi:10.1186/s13014-016-0697-4) contains supplementary material, which is available to authorized users.


A A.1 Registration Framework
Let R : R 3 → R denote the fixed reference image and T : R 3 → R the moving template image with compact support in domain Ω ⊂ R 3 . The goal of image registration is to find a plausible transformation y : Ω → R 3 that encodes the spatial correspondence between the two images R and T . In variational approaches, this is modeled by an objective function J which typically consists of a distance measure D describing image similarity and a regularizer S which penalizes implausible deformations [1]. Here, the transformation y is determined by solving the optimization problem where α > 0 balances between image similarity and deformation regularity of y and T (y) := T • y denotes the deformed template image.
Distance Measure Given the CBCT to CT registration problem, the normalized gradient fields distance measure D is chosen which is suitable for multimodal registration. This distance measure is based on the assumption that in two given images, regardless of modality, object edges and boundaries are always represented by intensity changes, i.e. image gradients. Furthermore it is easy to interpret and can be implemented efficiently and fully parallelized [2,3]. With · ε := ·, · + ε 2 , we define the NGF distance measure as where τ, ≥ 0 are so-called edge parameters, that allow for image modality specific filtering between noise and relevant image gradients.
Regularizer The second term in the objective function (1) is a regularizer.
Here, curvature regularization is used [4]. With the decomposition y(x) = x + u(x), the curvature regularizer is given by where "∆" denotes the Laplace operator and u x , u y , u z are the components of the displacement in x-, y-and z-direction, respectively.

A.2 Local Rigidity
The framework described in Section A.1 yields a non-linear deformation y that optimizes image similarity based on image gradients and deformation regularity based on second order derivatives. However, these deformations do not incorporate any additional knowledge about special anatomical deformation properties. Therefore, as proposed in [5] an additional constraint is added to the optimization problem in (1). Given a segmentation of a stiff structure in the planning CT, we define the domain of the segmentation as Σ ⊂ Ω. To model the stiff behavior of this structure, we incorporate a constraint which enforces rigid (i.e. rotational and translational) deformations on this sub-domain. Writing the local rigidity constraint as C(y) = 0 and assuming that it is active on Σ leads to s.t. C(y)(x) = 0 for all x ∈ Σ.
As there can be several rigid structures acting independently of each other, in general there will be multiple disjoint sub-domains Σ k , k = 1, . . . , M . To restrict structures to rigid deformations only, the constraint term C(y) is formulated as for all x ∈ Σ k and k = 1, . . . , M, where Q defines a three-dimensional rotation matrix, depending on the rotation angles θ = (α, β, γ) and the components of the translation vector b = (b 1 , b 2 , b 3 ) .
Constraints elimination Using the formulation in (3), constraints in the optimization problem (2) can be eliminated by substituting the deformation y. Instead of directly optimizing y, all the points that are in some Σ k are deformed in terms of a rigid deformation and its parameters, yielding where y 0 describes the deformation of all points that are outside of the Σ k 's and thus should deform non-rigidly, while y 1 , . . . , y M describe rigid deformations of points that are inside Σ 1 , . . . , Σ M , respectively. Substituting this in the objective function in (2), the optimization problem becomes unconstrained, and can be written as Here, the optimization is now no longer performed with respect to the deformation, but with respect to a new parameter set [y 0 , (θ 1 , b 1 ), . . . , (θ M , b M )] where each rigid region Σ k is controlled by six parameters in θ k , b k . Since each Σ k is a subset of Ω, the desired rigid areas are defined in the reference image space.

A.3 Discretization and implementation
The optimization problem described in A.1 and A.2 is solved using a discretizethen-optimize strategy [1]. The derived continuous model is first discretized on a grid and the obtained function is then optimized using well-established techniques from numerical optimization. To solve the discretized optimization problem, an iterative L-BFGS algorithm [6] is used for minimization. The L-BFGS algorithm is then embedded in a multi-level optimization scheme, where the optimization problem is consecutively solved on coarser to finer grids.
Here, a special strategy is used to prolongate the result obtained on a coarser level to a starting guess on a corresponding fine level. First, the full deformation field is prolongated as it would be using a standard deformable approach. After prolongation a new match of grid points to rigid areas is determined, generating a new set of unconstrained grid points while keeping the same number of rigid parameters θ k , b k , k = 1, . . . , M , since the number of rigid regions does not change. The result from the coarser level is then used as a starting guess for those rigid parameters.