Inelastic Deformable Image Registration (i-DIR): Capturing Sliding Motion through Automatic Detection of Discontinuities

Abstract

Deformable image registration (DIR) is an image-analysis method with a broad range of applications in biomedical sciences. Current applications of DIR on computed-tomography (CT) images of the lung and other organs under deformation suffer from large errors and artifacts due to the inability of standard DIR methods to capture sliding between interfaces, as standard transformation models cannot adequately handle discontinuities. In this work, we aim at creating a novel inelastic deformable image registration (i-DIR) method that automatically detects sliding surfaces and that is capable of handling sliding discontinuous motion. Our method relies on the introduction of an inelastic regularization term in the DIR formulation, where sliding is characterized as an inelastic shear strain. We validate the i-DIR by studying synthetic image datasets with strong sliding motion, and compare its results against two other elastic DIR formulations using landmark analysis. Further, we demonstrate the applicability of the i-DIR method to medical CT images by registering lung CT images. Our results show that the i-DIR method delivers accurate estimates of a local lung strain that are similar to fields reported in the literature, and that do not exhibit spurious oscillatory patterns typically observed in elastic DIR methods. We conclude that the i-DIR method automatically locates regions of sliding that arise in the dorsal pleural cavity, delivering significantly smaller errors than traditional elastic DIR methods.

1. Introduction

Deformable image registration (DIR) is an image-analysis technique used to determine the optimal transformation that establishes the spatial correspondence of a point between two images. When constructing a DIR method, three key elements need to be defined: (i) the transformation model, (ii) the regularizer, and (iii) the similarity measure [1]. These elements allow for the classification of DIR methods, and the reader is referred to [2] for a complete review. In particular, in this work we are concerned with the ability of the method to capture large displacements in the optimal transformation between medical images. From this perspective, transformation models can be divided into continuous-displacement transformations [3], which are suitable for small-deformation problems, and incremental diffeomorphic transformations based on the integration of flow equations [4,5], which can capture large deformations in DIR problems. While diffeomorphic methods have proven advantageous in capturing the large-displacement kinematics in DIR, continuous displacement models have been preferred in the field of medical imaging, as they provide a simple and efficient computational framework to DIR [6].

DIR has essential applications in radiology, such as the fusion of an anatomical image with a functional image [7], image-guided radiotherapy [8], and in treatment and surgery planning [9]. DIR has proven fundamental in the study of the deformation mechanisms that take place at a regional level in human lungs, where the primary inputs for determining regional deformation are deformation measures based on the Jacobian matrix of the optimal transformation resulting from DIR of lung computed-tomography (CT) images [6,10]. DIR-based biomechanical analysis has revealed significant spatial differences in the magnitude, anisotropy, and heterogeneity of regional deformation in the lung of normal human subjects [10,11], measured in terms of volumetric change expressed either as a Jacobian determinant [12] or as regional volumetric strain and deformation invariants that quantify linear and surface changes [13]. Further, spatial patterns of regional deformation obtained from DIR have been found to significantly differ from normal lungs in asthmatic patients [14] and patients with chronic obstructive pulmonary disease [15], highlighting the potential of biomechanical analysis in understanding, and potentially detecting disease progression. Estimates of volumetric strain have been correlated with lung inflammation and injury in mechanically-ventilated lungs, suggesting that regional deformation obtained from DIR can be useful in the prevention of ventilation-induced lung injury in critical-care patients [16,17,18]. A fundamental limitation of current DIR techniques and libraries is the poor performance experienced when images display motion discontinuities such as contact and organ sliding. Sliding typically occur inside the human body due to the existing lubricated interfaces between internal structures in the thoracic cage (e.g., lungs, chest wall, heart) [10,19,20]; and between the liver and other abdominal organs (e.g., kidney, diaphragm) [21,22]. Interestingly, sliding in the lung fissures has been detected from computed-tomography (CT) images of the lungs using DIR methods, where supraphysiological levels of shearing deformation colocalize with the fissures [23]. While useful for anatomy detection purposes, no regional tissue distortion is expected to occur in sliding regions, which invalidates the accuracy of regional deformation estimates from traditional DIR methods in the lungs in regions close to fissures and the pleural cavity. The main responsible for such spurious deformation levels in sliding is the transformation model that most DIR methods assume, typically constructed using interpolation schemes that deliver globally continuous and smooth transformation mappings [24]. As a consequence, traditional DIR methods cannot capture material interface or motion discontinuities, thus hindering the accuracy of the image registration and the associated biomechanical analysis [11].

Traditional DIR techniques have been modified to capture sliding either by using alternative regularization terms, as well as enhanced transformation models [24,25,26,27]. One example of the former is the diffusion-based approach [28], where the normal component of the displacement field near the sliding boundaries is continuous, and a direction-dependent regularization term is assumed such that it penalizes jumps in the normal direction but allows for a discontinuous displacement field in the tangential direction [25]. This direction-dependent registration model shows good registration accuracy but underperforms when the intensity contrast near the boundaries is low, which can be the case of lobar fissures in CT images of the lung. A similar approach that employs local weighting and direction-dependent anisotropic diffusion smoothing resulted in more realistic displacement fields than methods using global smoothing regularization [26]. Alternatively, sliding motion in DIR has been approached by using novel transformation models that allow for discontinuities at predefined boundaries. One such example is the use of a linear combination of multiple B-spline functions and a sliding constraint [27]. This enhanced formulation of the classical free-form deformation (FFD) model [3] delivered accurate estimations of the displacement deformation field in 16 patients with lung cancer. A step further is the extension of the FFD free-form deformation method, which consists in enhancing B-spline basis functions with discontinuous functions that have jumps defined at the discontinuity surface [22,24], a formulation that has been termed extended FFD (XFFD). XFFD has shown to deliver high accuracy when registering synthetic images with strong sliding discontinuities, as well as lung and liver images where high levels of sliding are present. While incorporating additional constraints that account for organ sliding results in better deformation estimates, several limitations preclude them from their direct applicability in the analysis of regional strain deformation in the lungs. A key limitation presented by XFFD models, which is also shared by other B-spline methods [27,29] and diffusion-based methods [25,26], is the fact that they rely on the definition of the sliding boundaries prior to the DIR analysis, which is typically done using semi-automatic segmentation methods, and largely depends on expert knowledge of the spatial location of the discontinuity boundaries. The a-priori definition of the sliding boundaries and the need for fine grids to capture curved surfaces where sliding occurs largely limit the applicability of current DIR methods in registering lung images, where the pleural cavity and fissures have intricate surface geometries and may not be easy to detect by the non-expert user.

The scientific question that motivates this work reads as follows: Is it possible to accurately capture sliding motion in DIR without the predefined knowledge of the sliding boundaries? To answer this question, in this work we aim at proposing and validating an inelastic DIR (i-DIR) method that allows for the automatic detection of sliding boundaries and that can handle discontinuous sliding motion on such surfaces. 2. Materials and Methods 2.1. Deformable Image Registration Elastic Formulation

In the following we adopt a variational framework for DIR problems [2,30], which will be the starting point of the i-DIR formulation. LetΩ⊂^Rnbe a domain of interest (image support),R:Ω→Rbe the reference image andT:Ω→Rbe the target image. The DIR problem aims to establish an optimal transformationu:Ω→^Rnthat best aligns the reference and target images. To this end, we consider the functional spaceV:=^H1(Ω,^Rn)and define a similarity functionalD:V→Rthat penalizes differences between the reference image R and the resampled target imageT∘(id+u) . A popular choice for the similarity measure in mono-modal applications of DIR [31,32] is the sum of squared-differences

D[w]:=12_∫Ω ^{[T(x+w(x))−R(x)]2}dΩ,∀w∈V,

which we will consider throughout this work. We remark that other choices of image similarity models such as those based on cross-correlation and mutual information measures can also be included in this formulation [30,33]. Further, we define a regularizerS:V→Rthat provides smoothness to the optimal transformation as well as it avoids ill-posedness of the DIR problem. A popular choice due to its physical meaning is the elastic regularizer

S[w]:=_∫Ω ^We(∇w)dΩ,

where the elastic energy density takes the form

^We(∇w):=μ^{∥∇w+∇wT∥2}+λ2^(divw)2,

with ∇ the gradient operator, div the divergence operator, andλandμthe Lamé constants. With these elements, we define the elastic registration functional as

Π[w]:=αD[w]+S[w],

whereα>0is a weighting parameter. Then, the optimal transformationuis the minimizer of the elastic registration functional, and the DIR problem is formulated as the following variational problem: Findusuch that

Π[u]=minw∈VΠ[w].

We note that the optimal transformationu can be interpreted as a displacement field that maps a point between its locations in the reference and target images. Moreover, the choice of the elastic deformation energy as a regularization term confers the DIR problem the physical interpretation of an elasticity problem [34], which has been widely exploited in the literature [1]. Further, and based upon this physical interpretation of the optimal transformationu, we define the strain tensor operator

ε(u):=12∇u+∇u^T,

and we note that the elastic energy density defined in (3) can be rewritten as

^We(ε)=με:ε+λ2^(traceε)2,

where : signifies the tensor scalar (inner) product and trace represents the trace operator. We further define the stress tensor associated to this elastic energy by

σ(ε):=∂^We∂ε=2με+λtrace(ε)I,

whereIis the identity tensor.

2.2. The Inelastic Deformable Image Registration (I-Dir) Method

As discussed in the introduction, the elastic regularizer is not suited to handle discontinuities in the displacement field.As a result, sliding motion is not captured by traditional DIR methods. To address this limitation, here we draw ideas from the mechanics of inelastic solids, which aims at modeling inelastic deformation processes that result in localized softening in a solid. In the following, we briefly summarize the main ingredients of a traditional von Mises plasticity model, for a comprehensive review of the theory of inelastic solids we refer the reader to [35]. Inelastic deformation in metals is driven by shearing deformation mechanisms, where sliding in the plane of maximum shearing occurs when the shear stress in that plane overcomes a critical yield stress, resembling frictional sliding motion.

We note that the mechanical behavior of an inelastic solid is path dependent, which we represent through a time-dependence of the associated displacement field and strain tensors. To model inelastic deformation processes, we adopt the standard additive decomposition of the strain tensor

ε=ε^e+ε^p

whereε^ecorresponds to the elastic strain tensor which is assumed to disappear as the load is removed, andε^p is the inelastic strain tensor which captures permanent deformations (i.e., sliding) that will remain in the solid after the load is removed. A sketch of this traditional decomposition of deformations is included in Figure 1. For the purposes of image registration, we note the inelastic strain tensor will capture sliding that does not generate deformation in a tissue, and therefore we will quantify regional deformation solely based on the elastic strain tensor. The additive decomposition carries onto the instantaneous evolution of strain components, and we note that (9) implies that

ε˙=^ε˙e+^ε˙p,

where()˙indicates partial derivatives with respect to time. To reflect the path-dependent nature of inelastic solids, we consider the effective inelastic strainq∈M:=^L2(Ω,R) as the hardening internal variable. Following a thermodynamic formalism, we assume a free energy density function (for rate-independent plasticity) that extends the elastic energy density (7) and takes the form

A(ε,ε^p,q)=^We(ε−ε^p)+^Wp(q),

where we assume that the stored plastic energy takes the form^Wp=12H^q2, with H being the hardening modulus. Then, the elastic constitutive relation reads

σ(ε^e)=∂W∂ε^e=2με^e+λtrace(ε^e)I,

and the relation between the hardening internal variable and its thermodynamic conjugate stress is

_σc(q)=∂A∂q=Hq.

The inelastic behavior of a solid is modeled by defining a inelastic potential, also known as the yield function, which in the case of a von Mises solid takes the form

Φ(σ,_σc):=32_{sij sij}−_σc

where_sijare the components of the deviatoric stress tensorsdefined as

s(σ):=σ−12traceσI.

Following an associative plasticity framework [35], the evolution of the inelastic strain tensor is governed by the flow rule

^ε˙p=γ˙∂Φ∂σ(σ,_σc)=γ˙32s∥s∥

whereγ˙is the inelastic multiplier and the norm defined as∥(·)∥=(·):(·). The evolution of internal variables is also dictated by the gradients of the inelastic potential

q˙=−γ˙∂Φ∂_σc(σ,_σc)=γ˙.

Finally, an inelastic model must comply with the loading/unloading complimentary conditions

Φ(σ,_σc)≤0,γ˙≥0,γ˙Φ(σ,_σc)=0.

Complimentary conditions (18) are interpreted as follows: Inelastic deformations will occur (γ˙>0) only when the current stress state reaches the yield surfaceΦ(σ,_σc)=0 . Otherwise, the inelastic evolution must be null to comply with (18), i.e.,γ˙=0 , which in turn implies that the change in total deformation will correspond only to changes in elastic deformation, as governed by (10).

In general, the relationship between^ε˙pand q is not independent. Letσ¯be the effective stress defined as

σ¯=32_{sij sij}.

Then, the relation between^ε˙pand q is assumed to follow the Prandtl-Reuss flow rule that reads

_ε˙ijp=^ε¯˙p32_sijσ¯,

where^ε¯˙p stands for the evolution of the accumulated plastic strain. In view of the plastic flow rule, the accumulated plastic strain is equivalent to q, so by integration of (20) and assuming isotropic hardening, we have the following relation,

^ε¯p=_∫0t23_{ε˙ijp ε˙ijp}dt≡q

(see Appendix A.1 for details).

2.3. Time and Space Discretization

Givenε˙ , the set of equations (10), (12), (13), (16) and (18) constitutes an inelastic constitutive initial value problem, which has been traditionally solved using a return-mapping algorithm based on an implicit backward-Euler temporal discretization. To this end, the time variable is discretized in generic subintervals[_tn,_tn+1]. Then, a series of incremental problems are obtained, where the main variables of the inelasticity model are assumed to be known at the timet=_tn, and need to be solved fort=_tn+1 , giving rise to classical return-mapping algorithms. The details about the numerical discretization of return mapping algorithms can be found elsewhere [36]. Conveniently, the elastoplastic incremental problems can be reformulated as incremental variational (minimization) problems, which gives rise to the theory of variational updates in the computational solid mechanics community [37,38,39]. In the following, we draw ideas from the theory of variational updates in plasticity to formulate the inelastic DIR model. The general framework consists in formulating the evolution of an elastoplastic solid as a sequence of incremental variational minimization problems. To this end, we first integrate the flow rule (16) using a Backward-Euler scheme to obtain

ε_n+1p−ε_np=(_qn+1−_qn)32s(ε_n+1−ε_n+1p)∥s(ε_n+1−ε_n+1p)∥

Solving forε_n+1p from the non-linear equation (22) delivers an incremental update for the inelastic strain tensor, which we express as

ε_n+1p=ε_n+1∗p(ε_n+1,_qn+1)

which depends solely ofε_n+1and_qn+1. Based on this flow-rule update, we define the effective incremental energy density fort=_tn as [38],

_Wn(ε)=inf_qn+1gn(ε,_qn+1)

where,

_gn(ε,_qn+1)=A(ε,ε_n+1∗p(ε,_qn+1),_qn+1)−_An+Δt·^ψ∗|_qn+1−_qn|Δt

whereΔt=_tn+1−_tn, and^ψ∗ stands for the dual dissipation potential [38] that governs the time evolution of the hardening variable, which in our case is defined as^ψ∗=_σy|Δq|, withΔq=_qn+1−_qn . The minimization problem involved in the definition (25) is equivalent to the stationary condition

0∈∂A∂_qn+1+∂^ψ∗|_qn+1−_qn|Δt

which, for the rate-independent case reads (see details in Appendix A.2),

_σ¯n+1pre−3μΔq−_σc(_qn+1)=∂^ψ∗|_qn+1−_qn|Δt

where_σ¯n+1preis the elastic predictor for_σ¯n+1 . Substituting (13) into (27) we obtain

_σ¯n+1pre−3μΔq−H(_qn+Δq)−∂^ψ∗|Δq|Δt=0

where the sub-differential of^ψ∗is defined as:

∂^ψ∗=[−_σy,_σy]ifΔq=0_σyifΔq>0−_σyifΔq<0

The solution ofΔq from (28) involves two mutually exclusive steps, giving rise to a return-mapping algorithm which involves an elastic predictor and a plastic corrector steps, see Appendix A.3.

With the definition of the effective incremental energy density, we now postulate the inelastic DIR formulation as a sequence of effective variational problems. For a generic time step, the displacement fieldu_nis assumed to be known, and we find the displacement fieldu_n+1by solving the problem

_Πneff[u_n+1]=minw∈V_Πneff[w],

where the inelastic DIR functional reads

_Πneff[w]=_Sneff[w]+αD[w],

and the inelastic regularizing term takes the form

_Sneff[w]=_{∫Ω Wn}(ε(w)).

To solve the minimization problem (30) we consider the stationary condition

_Rn[w;v]:=ddϵ_{Πneff[w+ϵv]ϵ=0}=0,∀v∈V.

The residual in (33) takes the form

_Rn[w;v]:=α_∫Ωv·(T(w)−R)∇T(w)+_∫Ωε(v):_σn+1(ε(w)),

where

σ_n+1(ε):=∂_Wn∂ε(ε),

represents the stress tensor update [38]. The residual Equation (33) constitutes a nonlinear problem, which we approach by means of linearization. To this end, we consider the Gauteaux differential defined as

TR[w,v;Δw]:=α_∫Ωv·∇T(w)⊗∇T(w)+(T(w)−R)∇∇T(w)·Δw+_∫Ωε(v):_Dn+1epε(Δw),

where

_Dn+1ep(ε):=^∂2 _Wn∂ε²,

is the consistent tangent tensor, see Appendix A.4. Thus, the linearized version of the residual problem reads: Given an initial guessw∈V, find the incrementΔwsuch that

_Rn[w;v]+TR[w,v;Δw]=0∀v∈V,

and we iterate over this linearized problem until a convergence criterion is reached.

To solve the continuous linear variational problem defined in (38) we adopt a Ritz-Galerkin finite-element approach. To this end, we construct the finite-element space

^Vh=^vh:^Ωh→^{Rn vh}:=∑A=1m_{NA vA},with_vA∈^Rn⊂V,

where{_N1,…,_Nm} is the set of basis functions. Using this finite-element space, we approximately solve the variational problem (38), i.e., we solve the problem: Given an initial guessu^h, find the incrementΔu^hsuch that

_Rn[w^h;v^h]+_TRn[u^h,v^h;Δu^h]=0∀v^h∈^Vh.

Using standard arguments (e.g., see [40]) we can show that (40) is equivalent to solving the linear system of equations

_KnΔu=_Fn,

whereΔuis a vector with the nodal values of the incrementΔu^h,_Knis the tangent matrix and_Rn is the residual for the previous guess, all of which are defined in Appendix A.5. After convergence is reached for the Newton step, the internal variables q at_tn+1 are updated and stored at the element level. We further note that, in order to provide stability and unisolvence of the problem, we adopt the approach set forth in [30], where we impose orthogonality conditions to the displacement fields and assume Neumann boundary conditions.

2.4. Performance Assessment and Metrics

The i-DIR method was implemented using an in-house Python code. In order to contrast the results of the i-DIR method with other DIR methods, we considered the open source Nifty Reg library [41] which efficiently implements the FFD method [3] with elastic regularization. To understand the effect of the inelastic regularization term over the purely elastic counterpart for a FE method, we also consider the comparison with an Elastic FEM registration. To study the performance of the three methods considered here (FFD, Elastic FEM, and i-DIR), we constructed synthetic reference and target images that simulated planar sliding over a chessboard-like image studied by Rua and co-workers [22], which we refer to as the synthetic dataset with sliding motion, see Figure 2. The synthetic images have a resolution of80×80 pixels, and the target image is constructed in such a way that it resembled a dislocation or sliding motion, with a known displacement of 5 pixels. We remark that, as the motion corresponds to a uniform vertical displacement (i.e., rigid body motion in both blocks) of the right-hand side of the image, the exact strain field is equal to zero, as rigid motions do not generate strain. To assess the method’s performance on anatomical images, we considered sagittal planes of CT thorax images of a normal volunteer under spontaneous breathing at total lung capacity (reference image) and functional residual capacity (target image), see Figure 3. The images were randomly selected from a small CT lung dataset of normal subjects employed in a previous study [13], and the sagittal planes were arbitrarily chosen so that large deformations were explicitly depicted to capture sliding.

To quantitatively evaluate the performance of the DIR methods, we considered the traditional residual sum of squared differences (RSS) between the reference and resampled images, defined as:

RSS=∑i=1m∑j=1n^{(_Rij−_{(T∘(id+u))ij})2}.

In addition, the normalized target registration error (TRE) was also computed. The TRE is defined as

TRE=_∑i=1N^(_pi−_qi)2N

where_pi(x,y)and_qi(x,y) are the ith landmark in the target image (fixed landmark) and the moving landmark, respectively, and N is the total number of landmarks. For the synthetic dataset, 375 landmarks were positioned around the discontinuity surface, as shown in Figure 2. In the case of the lung dataset, two sets of landmarks were considered, to analyze the effect of landmark selection and position regarding the sliding surface. The first case considered 9 landmarks positioned entirely inside the lung (Figure 3, top row). The second case considered 12 landmarks placed in the ribs of the dorsal region (Figure 3, bottom row).

In addition to the RSS and TRE metrics, the resampled image (T∘(id+u)), difference image (R−T∘(id+u)), and warped reference image ((φ,R)) are reported for all cases. To study the mechanical performance of all the methods studied, we constructed images of the elastic volumetric strain, defined as

_εvole:=trace(ε^e),

and images of the elastic von Mises strain, which takes the form

_εvme:=23ε^e:ε^e.

We note that for purely elastic methods (FFD, Elastic FEM), the elastic strain tensorε^eis replaced by the total strain tensorε in (44) and (45). Finally, a sensitivity analysis is conducted for the i-DIR method on the synthetic dataset to understand the effect of the initial yield stress on its registration performance.

2.5. Parameter Settings

For the FEM models (elastic FEM and i-DIR), we established an incremental approach for the weighting parameterα. We set an initial value ofα=0.01, and once a convergence tolerance was exceeded, we systematically increasedα, until we reach a value ofα=1400. Values for the Lamé constants (μ,λ), the initial yield limit (_σ0 ) and the hardening modulus (H) are included in Table 1.

Following the pyramidal approach described in [41], the FFD model consisted of a global registration and three consecutive local registration processes. The penalty term associated with the FFD model, also known as bending energy (BE), was set to:B_E1=1×¹⁰⁻⁹,B_E2=2.5×¹⁰⁻⁶andB_E3=1×¹⁰⁻⁴, for the three local registration respectively.

In terms of numerical discretization, both the elastic FEM and i-DIR models, employed structured triangular finite element meshes. As shown in Figure 4 the synthetic dataset used a mesh of size 5618 elements and for the lung dataset a mesh of size 18,860 elements. For visualization purposes, we further refined our results into structured meshes of size 64,800 and 28,800 elements for the synthetic and lung dataset, respectively. As for the FFD model, we projected the deformation mapping field (output) and compute the mechanical measures into a refined structured triangular finite element mesh of the same size as the FEM DIR models, in order to have a fair comparison.

3. Results 3.1. Synthetic Dataset with Sliding Motion

The performance of each registration model using the synthetic dataset with sliding motion is reported in terms of resampled and difference images in Figure 5. The i-DIR method accurately captures the vertical sliding and delivers the best resampled image, when compared to the other elastic methods (Figure 5, top row). Most of the errors in the resampled images are located in a small neighborhood of the line where sliding takes place. When analyzing the difference images (Figure 5, bottom row), the FFD method results in considerable voxel-wise differences at the boundaries of the squares that propagate from the sliding line throughout the checkerboard domain. In contrast, small differences are observed around the sliding line in the elastic FEM case. No visible differences are observed for the i-DIR case when compared to the other two methods.

Warped reference images showing the resulting displacement field for each method are reported in Figure 6. A close-up around the sliding region shows a continuous displacement field with a vortex-like pattern over the sliding line that slowly dissipates to the right for the FFD case. A similar displacement field pattern is observed for the elastic FEM case, but with an attenuated vortex pattern. In contrast, the i-DIR method delivers a uniformly vertical displacement field on the region to the right of the sliding line, and zero displacements to the left of the sliding line, being able to identify the discontinuity surface as well as capturing the discontinuous displacement field.

The RSS and TRE metrics for all three methods are shown in Table 2. The i-DIR method delivers the lowest values for these performance metrics, followed by the Elastic FEM method.

Figure 7 shows the elastic deformation fields associated to the three registration methods. The elastic volumetric strain displayed by the FFD model, shows an erratic pattern throughout the entire image, with high values of both compressive (peak value of −1.22) and expansive (peak value of 0.83) deformation near the vicinity of the discontinuity surface. The case of the elastic FEM model delivers a volumetric strain field with localized strain concentrations around the sliding surface with peak values of −0.75 and 1.27. In contrast, the i-DIR model delivers a volumetric strain field that is zero in the majority of the region of analysis, with small concentrations around the sliding surface with peak values of −0.39 and 0.63. The resulting elastic von Mises strain field, which characterizes shear distortions, is shown in the bottom row of Figure 7. Similarly to the case of volumetric strain, the FFD method results in a highly oscillating field that take on non-zero values everywhere in the image domain, reaching peak values of 1.72. The Elastic FEM method displays high strain concentrations around the sliding plane with peak von Mises strain values that are similar to the FFD case (2.38), but the strain field rapidly dissipate away from the discontinuity plane. In contrast, the i-DIR results in a narrow region around the sliding plane with low values (peak value of 0.38), with the rest of the image domain resulting in zero von Mises strain.

The sensitivity of the RSS error to the value of the initial yield stress in the i-DIR method is shown in Figure 8. Yield stress values smaller that0.1result in RSS errors that do no change considerably, delivering the highest accuracy observed for all three methods. In contrast, yield stress values above1.0 deliver a much higher RSS error, which approaches that of the Elastic FEM method, see Table 2.

3.2. Registration of Lung CT Images

Resampled and difference images for the lung CT dataset are shown in Figure 9, top row. All three methods deliver similar results of the resampled image. We note however, that resampled images from both the FFD and Elastic FE methods show distorted rib cuts in the dorsal region, while the ribs are accurately resampled in the case of the i-DIR method. The misalignment of the ribs is also observed in the difference images of the FFD and Elastic FEM cases, see Figure 9, bottom row. In contrast, the i-DIR case reports zero difference values in the regions where ribs are located.

Warped reference images are shown in Figure 10, where a close-up shows the displacement fields around the sliding pleural cavity. A continuous, and almost uniform upward displacement field is observed for the case of the FFD and Elastic FE methods. In contrast, the i-DIR method delivers an upward displacement field inside the lung, right next to a region comprising the ribs with null displacement, with the jump in displacement magnitude located on the sliding pleural cavity.

Performance metrics for the lung dataset are included in Table 3. We note that all three methods result in similar values for the case of RSS and TRE using inside-lung landmarks. However, the i-DIR method shows a remarkable advantage over the other methods for the case of TRE using rib landmarks.

The elastic volumetric strain distribution resulting from the registration of lung images are shown in Figure 11, top row. The FFD model delivers a highly oscillating field that results in excessive strain values with peaks as high as −1.68 and 1.02, located both inside and outside the lung domain. In contrast, the Elastic FEM model displays a more uniform volumetric strain distribution inside the lung, with a smooth pattern of strain. However, high strain localizations are observed outside the lung in the dorsal region where the ribs are located, with oscillating values. The i-DIR model delivers a smooth distribution of volumetric strain inside the lung, that quickly transitions to small levels os strain immediatly outside the lung. Further, the largest strain levels are found in the regions near the diaphragm. Outside the lung, we mostly observe zero volumetric deformation throughout the remaining image domain. The von Mises strain fields are shown in Figure 11, bottom row. Similary to the case of volumetric strain, the FFD model delivers a highly oscillating field with a peak value in the order of 1.7 both outside and inside the lung. The Elastic FEM method results in a distribution with smaller strain magnitudes, which in some parts of the lung boundary are rapidly reduced to zero. In the case of the i-DIR method, a smooth distribution of non-zero strain is observed inside the lung with the highest values close to the diaphragm and dorsal region. The von Mises strain distribution sharply decays to zero in the regions outside the lung.

4. Discussion

The results for the synthetic dataset with sliding motion show that among the three DIR methods studied, the i-DIR delivers the best resampled image, accurately accomodating the sliding motion, see Figure 5, top row. We note that the elastic DIR methods suffer from spurious displacements around the sliding line that result in distorted resampled images, see Figure 5, bottom row. Further, the i-DIR method is capable of capturing the discontinuous displacement field imposed by the sliding motion, while elastic DIR methods fail to capture the jump in displacements and result in spurious displacement fields, see Figure 6. Further, we have shown that for this example, the i-DIR method consistently delivers RSS and TRE metrics that confirm the superior performance of the i-DIR method when compared to the Elastic FEM and FFD methods, see Table 2. Other methods proposed in the literature have also shown a remarkable performance in the registration of the synthetic dataset with sliding motion. In particular, the XFFD method has shown to be capable of accurately capture the sliding motion by introducing a discontinuous transformation model that delivers an optimal resampling [22]. To this end, the XFFD method necessitates the definition of the sliding surface a priori in order to deliver accurate and efficient results. Here, we have shown that the i-DIR method does not require a priori information about the sliding surface. Further, the sliding plane did not coincide with any element edges in the discretization. This feature represents an important advantage over existing methods, as the i-DIR is capable of detecting sliding discontinuities in an automatic way, lending itself to the registration of images with arbitrary sliding discontinuities.

From the perspective of quantifying local deformation by means of DIR, we remark that the sliding mechanism present in the synthetic dataset corresponds to a rigid (sliding) motion between two adjacent blocks, and therefore no deformation is expected to occur in any of the blocks after sliding. Figure 7 shows that the FFD method induces spuriously high levels of both volumetric and deviatoric deformations around the sliding plane, which is consistent with previous findings for the synthetic dataset with sliding motion [22]. Similarly to the case of the displacement field, the strong warping and deformation propagates throughout the image domain, creating nonphysical high volumetric and deviatoric strain levels away from the discontinuity. The error in the strain predictions is strongly attenuated by the elastic FE method, which still concentrates high values in a neighborhood of the sliding plane. In contrast, the i-DIR method reports low levels of deviatoric deformation on a narrow band around the discontinuity surface, and negligible errors in the estimation of volumetric strain, see Figure 7. This result shows that the i-DIR method not only automatically captures the sliding motion with accuracy, but also delivers precise estimates of the strain fields, even in the present of strong discontinuities.

The sensitivity of the iDIR model to the yield stress parameter shows that for parameter values_σy≤¹⁰⁻²no appreciable improvement is obtained in terms of the RSS error. Further, we note that high values of yield stress deliver errors that are equal to those reported by the Elastic FEM method. These results show that tuning the yield stress parameter is essential for obtaining accurate results from the registration process.

The i-DIR method was also assessed in the analysis of medical CT images of the lung, where sliding is expected to occur when registering images from resting states to maximal inspiration effort [23]. When comparing resampled images, we showed that the i-DIR method delivered errors in registering the domains inside the lung that are comparable to those found in elastic DIR methods, see Figure 9. This conclusion is supported by the performance metrics RSS and TRE for the case of inside-lung landmarks reported in Table 3, where no marked differences are observed among the FFD, Elastic FEM and i-DIR methods. However, when assessing anatomical structures that are outside the lung, i.e., ribs, we observe that the i-DIR accurately resamples them to the correct location, whereas elastic DIR methods fail to achieve a reasonable result, see Figure 9, bottom row. This observation is confirmed by the results obtained in the TRE when using rib landmarks, where the i-DIR delivers errors that are one order of magnitude smaller than the error provided by elastic DIR methods, see Table 3. Once again, we attribute the good performance of the i-DIR method to its ability to handle discontinuous sliding motion by considering inelastic deformations in those regions, a feature not displayed by elastic methods, see Figure 10.

The evaluation of the elastic strain fields from registering lung CT images results in conclusions similar to those obtained in the case of the synthetic image: elastic DIR methods introduce highly oscillating fields both for the volumetric and deviatoric components of the elastic strain tensor, see Figure 11. We remark here that previous works on lung image registration have reported oscillatory strain fields when using elastic DIR methods that employ B-splines or other smooth and continuous basis functions for the construction of the deformation model [22,42]. However, these oscillations have been shown to hinder the accuracy of the estimations of local pulmonary deformation, as they arise due to the inability of the deformation model to capture discontinuities in the displacement field [11]. Notably, the i-DIR method delivers a smooth distribution of elastic strain inside the lung domain, with a sharp decay outside that approaches a state of no deformation. Further, the spatial patterns of volumetric strain delivered by the i-DIR method are in good agreement with those reported in the literature for normal human lungs that have been analyzed by isolating the lung domain [13], where larger volumetric strains are observed in the dorsal (dependent) and basal regions of the lung.

In conclusion, we have introduced a novel inelastic model for DIR that automatically captures sliding without a priori knowledge or assumptions about the spatial location of discontinuities or the need of a segmentation to denote the slipping domain. We note that other DIR formulations have been proposed in the literature to handle motion discontinuities without the need of segmenting the sliding region before the analysis [25]. In order to do so, these methods consider some assumptions related to the specific physiological behavior of the organ. For the case of the lungs, slipping motion is restricted to the edges of the image, and slippage occurs along the edge of the image. We note that these assumptions are not required by the i-DIR method, and therefore it represents a truly automatic technique for the detection of discontinuities in DIR of arbitrary images. The key ingredient to achieve this performance is the introduction of an inelastic energy term, which automatically locates regions of high shearing deformation associated to sliding and locally modifies the effective mechanical properties, allowing for higher levels of shear deformation in localized domains. We remark that, while inelastic formulations are standard in the field of computational mechanics [38], the inclusion of inelastic energy regularizers is novel in the field of image analysis, and, to the best of our knowledge, has not been pursued in the past in in the field of DIR. For the application of lung images, it is worth mentioning that we aim at automatically capture sliding, particularly between the lung boundary and the ribs (clearly identified within the images), and not necessarily improve the registration accuracy inside the lung. The above is supported by the results obtained when measuring the RSS error, which demonstrates that our i-DIR model holds a comparable performance with traditional DIR methods, especially in areas with no sliding. However, the inelastic model is considerably superior in capturing slippage at the lung edges, which is again substantiated by a better performance when measuring the TRE using rib-landmarks.

The present work can be extended in several directions. One limitation of the current computer implementation of the i-DIR method is the large wall-clock time required to solve the optimization problem, which can take up to 40 times the time required by optimized elastic DIR methods. This limitation may be alleviated by implementations that leverage the power of GPUs in DIR libraries [41]. In addition, due to the high computational demands, the current version of the i-DIR method has only been applied to 2D images. We remark that the motion and deformation analysis based on 2D images of the thorax constitutes an important limitation of this work in the biomechanical characterization of the lung. However, we also remark that under normal conditions, the dominant orientation of displacements in the lung is in the apico-basal direction [10], which is included in the sagittal images considered in this study. Future extensions should focus on DIR implementations for 3D CT thoracic images, based on which a complete biomechanical study can be performed to fully understand the 3D nature of deformations in the lung. In addition, we note that due to the large level of strains experienced by the lung under full inspiration, the elastic energy component employed in the i-DIR formulation may not be suitable, as it corresponds to the elastic deformation energy for small strain levels [35]. To overcome this limitation, hyperelastic warping fomulations have been proposed which employ elastic energy terms that are compatible with large deformation [43]. The use of hyperelastic energy terms in the future versions of the i-DIR method constitutes a promising avenue of research.

Generate image description

Figure 1. Schematics comparing the behaviour between an elastic and an inelastic approach.

Generate image description

Figure 2. Synthetic dataset with sliding motion: Reference (left) and target (right) images. Landmarks used for computing the target registration error (TRE) are shown in red for the reference image and in green for the target image.

Generate image description

Figure 3. Lung dataset: Reference (left) and target (right) images. (Row A): TRE analysis using landmarks inside the lung. (Row B): TRE analysis using using landmarks on the dorsal ribs. Red and green marks indicate landmarks in the reference and target images, respectively.

Generate image description

Figure 4. Numerical discretization of the synthetic and lung datasets. (Column A) Structured triangular meshes and (Column B) refined structured triangular meshes.

Generate image description

Figure 5. Registration of synthetic dataset with sliding motion. (Top row) resampled images using free-form deformation (FFD), Elastic finite element method and inelastic deformable image registration (iDIR) methods and reference image, (bottom row) difference images. Colorbar indicates the absolute intensity difference between images.

Generate image description

Figure 6. Warped reference image and displacement field for the synthetic dataset with sliding motion. Red arrows show the displacement field in a neighborhood of the sliding plane.

Generate image description

Figure 7. Elastic deformation fields for the synthetic dataset with sliding motion resulting from the different registration methods: elastic volumetric strain (top row, colorbar displays strain magnitude), and elastic von Mises strain (bottom row, colorbar displays strain magnitude).

Generate image description

Figure 8. Sensitivity of the i-DIR method, measured in terms of residual sum of squared differences (RSS) error, to the choice of initial yield stress value.

Generate image description

Figure 9. Registration of the lung dataset and comparison between methods. (Top row) resampled images, (bottom row) difference images. Colorbar indicates the absolute intensity difference between images. Reference image is included for comparison purposes.

Generate image description

Figure 10. Warped reference image and displacement field for the lung dataset. Red arrows show the displacement field in a neighborhood of the sliding surface.

Generate image description

Figure 11. Elastic deformation fields for the lung dataset resulting from the different registration methods: elastic volumetric strain (top row, colorbar displays strain magnitude), and elastic von Mises strain (bottom row, colorbar displays strain magnitude).

Parameter	Value
μ	1.36
λ	0.34
σ0	0.1
H	0.3

Model	RSS	TRE
FFD	16.32	1.17
Elastic FEM	3.19	0.46
i-DIR	0.28	0.22

Model	RSS	TRE	TRE
		(Inside-Lung Landmarks)	(Rib Landmarks)
FFD	11.64	6.82	13.98
Elastic FEM	13.25	6.74	13.68
i-DIR	12.24	6.99	0.77

Author Contributions

Conceptualization, D.E.H.; methodology, C.I.A. and D.E.H.; software, C.I.A.; validation, C.I.A. and D.E.H.; formal analysis, C.I.A. and D.E.H.; investigation, C.I.A. and D.E.H.; resources, D.E.H.; writing-original draft preparation, C.I.A. and D.E.H.; writing-review and editing, D.E.H.; visualization, C.I.A.; supervision, D.E.H.; project administration, D.E.H.; funding acquisition, D.E.H. Both authors have read and agreed to the published version of the manuscript.

Funding

This work received financial support from the Chilean National Agency for Research and Development (ANID) through grant FONDECYT Regular#1180832awarded to D.E. and by the Millennium Science Initiative Program-NCN17_129. C.I.A. acknowledges the support of ANID through the CONICYT Doctoral fellowship.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

i-DIR	Inelastic deformable image registration
FEM	Finite element method
FFD	Free form deformation
TLC	Total lung capacity
FRC	Functional residual capacity
CT	Computed tomography
RSS	Residual sum of squared differences
TRE	Target registration error

Appendix A. Mathematical Definitions and Demonstrations

Appendix A.1. Relation between the Rate of Plastic Strain and the Internal Variables

In general, the relationship between^ε˙pand q is not independent and is given by the Prandtl-Reuss flow rule:

_ε˙ijp=^ε¯˙p32_sijσ¯

where^ε¯˙pstands for the evolution of the accumulated plastic strain andσ¯for the effective stress, which for the von Mises model takes the form

σ¯=32_{sij sij}.

For multi-axial plasticity we can rewrite the flow rule (20) as,

_ε˙ijp=^ε¯˙p _Mij

whereM≡_Mij=32_sijσ¯ stands for the instantaneous direction of plastic flow. Similarly, the effective stress (A2) can be expressed as

σ¯=_{σij Mij}.

Going a step back and clearing^ε¯˙p from (A1) we have,

^ε¯˙p=23_{ε˙ijp ε˙ijp}

which in view of the plastic flow rule, is analogous to

^ε¯˙p=q˙

Then by integration of (A5) and again, assuming isotropic hardening, we have the following relation,

^ε¯p=_∫0t23_{ε˙ijp ε˙ijp}dt≡q

Finally, we compute the driving forces for q as,

y=-∂A∂q=∂^We∂_εijeMij-∂^Wp∂q(q)=_{σij Mij}-∂^Wp∂q(q)

y=σ¯-_σc.

Appendix A.2. Incremental Flow Rule Update

Following an incremental flow rule of the type,

ε_n+1p=ε_np+ΔqM

=ε_np+Δq32s_n+1σ¯n+1

let

_gn(ε_n+1,_qn+1)=A(ε_n+1,ε_n+1p(_qn+1),_qn+1)-_An+Δt·^ψ∗|_qn+1-_qn|Δt

Then we seek to minimize (A12) with respect to_qn+1, such that,

∂_gn∂_qn+1=0⟹inf_qn+1gn(ε_n+1,_qn+1)

we solve the above in a sub-differential way, such that

0∈∂A∂_qn+1+∂^ψ∗|_qn+1-_qn|Δt

0∈∂^We∂ε_n+1e·∂ε_n+1e∂_qn+1(ε_n+1-ε_np-(_qn+1-_qn)M)+∂^Wp∂q(_qn+1)+∂^ψ∗|_qn+1-_qn|Δt

0∈-σ_n+1·M+∂^Wp∂q(_qn+Δq)+∂^ψ∗|Δq|Δt

0∈-_σ¯n+1+∂^Wp∂q(_qn+Δq)+∂^ψ∗|Δq|Δt

Recalling that the von Mises flow vector is purely deviatoric [36], we have that

s_n+1=2με_n+1e=2μ(ε_n+1-ε_n+1p)=2με_n+1-ε_np-Δq32s_n+1σ¯n+1=2με_n+1-ε_np-2μΔq32s_n+1σ¯n+1

s_n+1=s_n+1pre-3μΔqs_n+1σ¯n+1

Since the predictive and updated deviatoric stress are co-linear (s_n+1‖s_n+1pre), we can state that,

M=32s_n+1σ¯n+1=32s_{n+1preσ¯n+1pre}

⟹s_n+1=s_{n+1preσ¯n+1preσ¯n+1}

From the above we can re-write (A18) as,

s_n+1=s_n+1pre-3μΔqs_{n+1preσ¯n+1pre}

s_n+1=1-3μΔq_σ¯n+1pres_n+1pre

Now replacing (A20) in (A22), we have,

s_{n+1preσ¯n+1preσ¯n+1}=1-3μΔq_σ¯n+1pres_n+1pre

⟹_σ¯n+1=_σ¯n+1pre-3μΔq

Going a step back to (A17) we can explicitly define,

_σ¯n+1pre-3μΔq∈∂^Wp∂q(_qn+Δq)+∂^ψ∗|Δq|Δt

Finally, assuming that^Wp=12H_qn+12and^ψ∗=_σy|Δq| , we can rewrite (A25) as:

_σ¯n+1pre-3μΔq-H(_qn+Δq)-_σy=0

which eventually delivers,

Δq=_σ¯n+1pre-H_qn-_σy3μ+H

Appendix A.3. Return Mapping Algorithm

The solution ofΔq in (A27) involves a two mutually exclusive steps:

(i)

An elastic predictor, such that,

Δq=0_qn+1pre=_qnσ¯n+1pre={2μ(ε_n+1-ε_n+1p(_qn+1pre))+λtrace(ε_n+1-ε_n+1p(_qn+1pre))I}·M=32s_n+1pre·s_n+1pre

replacing in (A26), we have

_σ¯n+1pre-H_qn-[-_σy,_σy]=0

⟹_σ¯n+1pre≤H_qn+_{σyσ¯n+1pre}≥-H_qn-_σy

and

(ii)

A plastic corrector, where we have two possible cases:

(a) if the elastic trial lies within the elastic domain

Φ(_σ¯n+1pre)≤0⟹_σ¯n+1pre∈[-_σy,_σy]

there is no plastic evolution within the time interval(_tn,_tn+1), and therefore we update our variables:

_(·)n+1=_(·)n+1pre

and

(b) otherwise, we have plastic flow (or elasto-plastic evolution). By a traditional Newton-Raphson linearization we solve the following

_σ¯n+1pre-3μΔq-H(_qn+Δq)-_σy=0

and then we update the following variables at_tn+1,

ε_n+1p=ε_np+Δq32s_n+1σ¯n+1

_qn+1=_qn+Δq

ε_n+1e=ε_n+1-ε_n+1p

σ_n+1=λtrace(ε_n+1e)I+2με_n+1e

_Dn+1ep=2μ1-3μΔq_σ¯n+1preId+6^μ2Δq_σ¯n+1pre-13μ+H_M¯n+1⊗_M¯n+1+KI⊗I

where K is the bulk modulus,_M¯n+1≡23M_n+1=s_n+1pre∥s_n+1pre∥is the unit plastic flow vector and_Idis the fourth order deviatoric projection tensor defined as,

_Id≡_IS-13I⊗I

with_IS=12(_{δik δjl}+_{δil δjk})as the fourth order symmetric identity tensor.

Appendix A.4. Effective Incremental Energy

The effective incremental energy can be defined as follows,

_Wn(ε)=inf_qn+1gn(ε,_qn+1)=_gn(ε,_qn+1∗(ε))

where,

_qn+1∗(ε)∈argmin_gn(ε,_qn+1)

such that,

∂_gn∂_qn+1(ε,_qn+1∗(ε))=0

Conveniently know the effective potential depends solely inε=ε_n+1. Following this definition, we compute the first and second derivative of_Wn, where,

D_Wn(ε)=∂_gn∂ε(ε,_qn+1∗(ε))

=∂_gn∂ε(ε,_qn+1∗(ε))+∂_gn∂_qn+1(ε,_qn+1∗(ε))·∂_qn+1∗∂ε(ε)

but since∂_gn∂_qn+1(ε,_qn+1∗(ε))=0, we have that for a fully implicit scheme,

D_Wn(ε)=∂_gn∂ε(ε,_qn+1∗(ε))≡σ_n+1

whereσ_n+1are the stresses att=_tn+1.

Then we compute,

^D2 _Wn(ε)=^∂2 _gn∂ε∂ε(ε,_qn+1∗(ε))+^∂2 _gn∂ε∂_qn+1(ε,_qn+1∗(ε))·∂_qn+1∗∂ε(ε)

where we can redefine∂_qn+1∗∂ε(ε) deriving (A42) by∂∂ε, such that,

^∂2 _gn∂_qn+1∂ε(ε,_qn+1∗(ε))+^∂2 _gn∂_qn+1∂_qn+1(ε,_qn+1∗(ε))·∂_qn+1∗∂ε(ε)=0

⟹∂_qn+1∗∂ε(ε)=-^{∂2 _gn∂_qn+1∂_qn+1(ε,_qn+1∗(ε))-1}·^∂2 _gn∂_qn+1∂ε(ε,_qn+1∗(ε))

replacing in (A46), we have that for a fully implicit scheme,

^D2 _Wn(ε)=^∂2 _gn∂ε∂ε(ε,_qn+1∗(ε))-^∂2 _gn∂ε∂_qn+1(ε,_qn+1∗(ε))·^{∂2 _gn∂_qn+1∂_qn+1(ε,_qn+1∗(ε))-1}·^∂2 _gn∂_qn+1∂ε(ε,_qn+1∗(ε))

^D2 _Wn(ε)≡_Dn+1ep

where_Dn+1epis known as the consistent tangent modulus att=_tn+1.

Appendix A.5. Finite-Element Discretization of the I-Dir Formulation

Let:

u≈_uih:=∑A=1n_{NA uiA}=N(x)u

v≈_vjh:=∑A=1n_{NA vjA}=N(x)v

ε(∇u)≈ε(∇u^h):=∑A=1n_{BA uiA}=B(x)u

ε(∇v)≈ε(∇v^h):=∑A=1n_{BA vjA}=B(x)v

where_NAare the shape functions,Nis the matrix of global shape functions, andBis the global strain-displacement matrix.

N=_N100_N2⋯_Ni000_N100⋯0_Ni000_N10⋯00_Ni

B=_B1B2B3⋯_Bi

_Bi=∂_Ni∂x000∂_Ni∂y000∂_Ni∂z∂_Ni∂y∂_Ni∂x00∂_Ni∂z∂_Ni∂y∂_Ni∂z0∂_Ni∂x

Substituting the approximations (A51)-(A54) into the linear variational problem (40) we obtain the linear system of equations defined in (41), where the tangent matrix and residual vector are defined as

Kn:=α∫Ωh NT∇T(unh)⊗∇T(unh)+(T(unh)-R)∇∇T(unh)N+∫Ωh BT Dn+1ep(ε(unh))B,

Fn:=α∫Ωh NT(T(unh)-R)∇T(unh)+∫Ωh BT σn+1(ε(unh)),

which are constructed by numerically evaluating the element expressions and assembling their contributions into the global matrix and vector using standard finite-element techniques, see, e.g., [40].

References

1. Sotiras, A.; Davatzikos, C.; Paragios, N. Deformable medical image registration: A survey. IEEE Trans. Med Imaging 2013, 32, 1153-1190.

2. Modersitzki, J. Numerical Methods for Image Registration; Oxford University Press: Oxford, UK, 2003; p. 199.

3. Rueckert, D.; Sonoda, L.I.; Hayes, C.; Hill, D.L.; Leach, M.O.; Hawkes, D.J. Nonrigid registration using free-form deformations: Application to breast MR images. IEEE Trans. Med. Imaging 1999, 18, 712-721.

4. Ashburner, J. A fast diffeomorphic image registration algorithm. NeuroImage 2007, 38, 95-113.

5. Christensen, G.; Rabbitt, R.; Miller, M. Deformable template using large deformation kinematics. IEEE Trans. Image Process. 1996, 5, 1435-1447.

6. Reinhardt, J.M.; Ding, K.; Cao, K.; Christensen, G.E.; Hoffman, E.A.; Bodas, S.V. Registration-based estimates of local lung tissue expansion compared to xenon CT measures of specific ventilation. Med. Image Anal. 2008, 12, 752-763.

7. Oliveira, F.P.; Tavares, J.M.R. Medical image registration: A review. Comput. Methods Biomech. Biomed. Eng. 2014, 17, 73-93.

8. Foskey, M.; Davis, B.; Goyal, L.; Chang, S.; Chaney, E.; Strehl, N.; Tomei, S.; Rosenman, J.; Joshi, S. Large deformation 3D image registration in image-guided radiation therapy. Phys. Med. Biol. 2005, 5869-5892.

9. Gering, D.T.; Nabavi, A.; Kikinis, R.; Grimson, W.E.L.; Hata, N.; Everett, P.; Jolesz, F.; Wells, W.M. An integrated visualization system for surgical planning and guidance using image fusion and interventional imaging. In Proceedings of the International Conference on Image Computing and Computer-Assisted Intervention-MICCAI, Cambridge, UK, 19-22 September 1999; pp. 809-819.

10. Amelon, R.E.; Cao, K.; Ding, K.; Christensen, G.E.; Reinhardt, J.M.; Raghavan, M.L. Three-dimensional characterization of regional lung deformation. J. Biomech. 2011, 44, 2489-2495.

11. Hurtado, D.E.; Villarroel, N.; Retamal, J.; Bugedo, G.; Bruhn, A. Improving the accuracy of registration-based biomechanical analysis: A finite element approach to lung regional strain quantification. IEEE Trans. Med. Imaging 2016, 35, 580-588.

12. Jahani, N.; Choi, S.; Choi, J.; Iyer, K.; Hoffman, E.A.; Lin, C.l. Assessment of regional ventilation and deformation using 4D-CT imaging for healthy human lungs during tidal breathing. Appl. Physiol. 2015, 1064-1074.

13. Hurtado, D.E.; Villarroel, N.; Andrade, C.; Retamal, J.; Bugedo, G.; Bruhn, A. Spatial patterns and frequency distributions of regional deformation in the healthy human lung. Biomech. Model. Mechanobiol. 2017, 16, 1413-1423.

14. Choi, S.; Hoffman, E.A.; Wenzel, S.E.; Tawhai, M.H.; Yin, Y.; Castro, M.; Lin, C.l. Registration-based assessment of regional lung function via volumetric CT images of normal subjects vs. severe asthmatics. Appl. Physiol. 2013, 115, 730-742.

15. Bodduluri, S.; Bhatt, S.P.; Hoffman, E.A.; Newell, J.D.; Martinez, C.H.; Dransfield, M.T.; Han, M.K.; Reinhardt, J.M.; COPDGene Investigators. Biomechanical CT metrics are associated with patient outcomes in COPD. Thorax 2017, 72, 409-414.

16. Retamal, J.; Hurtado, D.; Villarroel, N.; Bruhn, A.; Bugedo, G.; Amato, M.B.P.; Costa, E.L.V.; Hedenstierna, G.; Larsson, A.; Borges, J.B. Does Regional Lung Strain Correlate With Regional Inflammation in Acute Respiratory Distress Syndrome During Nonprotective Ventilation? An Experimental Porcine Study. Crit. Care Med. 2018, 46, e591-e599.

17. Hurtado, D.E.; Erranz, B.; Lillo, F.; Sarabia-Vallejos, M.; Iturrieta, P.; Morales, F.; Blaha, K.; Medina, T.; Diaz, F.; Cruces, P. Progression of regional lung strain and heterogeneity in lung injury: Assessing the evolution under spontaneous breathing and mechanical ventilation. Ann. Intensive Care 2020, 10, 107.

18. Cruces, P.; Retamal, J.; Hurtado, D.E.; Erranz, B.; Iturrieta, P.; González, C.; Diaz, F. A physiological approach to understand the role of respiratory effort in the progression of lung injury in SARS-CoV-2 infection. Crit. Care 2020, 24, 494.

19. Rodarte, J.R.; Hubmayr, R.D.; Stamenovic, D.; Walters, B.J. Regional lung strain in dogs during deflation from total lung capacity. Appl. Physiol. 1985, 58, 164-172.

20. Yin, Y.; Hoffman, E.A.; Lin, C.L. Lung Lobar Slippage Assessed with the Aid of Image Registration. In Medical Image Computing and Computer-Assisted Intervention-MICCAI 2010; Jiang, T., Navab, N., Pluim, J.P.W., Viergever, M.A., Eds.; Springer: Berlin/Heidelberg, Germany, 2010; pp. 578-585.

21. Von Siebenthal, M.; Székely, G.; Gamper, U.; Boesiger, P.; Lomax, A.; Cattin, P. 4D MR imaging of respiratory organ motion and its variability. Phys. Med. Biol. 2007, 52, 1547-1564.

22. Hua, R.; Pozo, J.M.; Taylor, Z.A.; Frangi, A.F. Multiresolution eXtended Free-Form Deformations (XFFD) for non-rigid registration with discontinuous transforms. Med. Image Anal. 2017, 36, 113-122.

23. Amelon, R.E.; Cao, K.; Reinhardt, J.M.; Christensen, G.E.; Raghavan, M.L. A measure for characterizing sliding on lung boundaries. Ann. Biomed. Eng. 2014, 42, 642-650.

24. Hua, R.; Pozo, J.M.; Taylor, Z.A.; Frangi, A.F. Discontinuous Non-rigid Registration using Extended Free-Form Deformations. In Proceedings of the SPIE Medical Imaging 2015: Image Processing, Orlando, FL, USA, 21-26 February 2015; Volume 9413.

25. Schmidt-Richberg, A.; Werner, R.; Handels, H.; Ehrhardt, J. Estimation of slipping organ motion by registration with direction-dependent regularization. Med. Image Anal. 2012, 16, 150-159.

26. Pace, D.F.; Aylward, S.R.; Niethammer, M. A Locally Adaptive Regularization Based on Anisotropic Diffusion for Deformable Image Registration of Sliding Organs. IEEE Trans. Med. Imaging 2013, 32, 2114-2126.

27. Delmon, V.; Rit, S.; Pinho, R.; Sarrut, D. Registration of sliding objects using direction dependent B-splines decomposition. Phys. Med. Biol. 2013, 58, 1303-1314.

28. Thirion, J.P. Image matching as a diffusion process: An analogy with Maxwell's demons. Med. Image Anal. 1998, 2, 243-260.

29. Wu, Z.; Rietzel, E.; Boldea, V.; Sarrut, D.; Sharp, G.C. Evaluation of deformable registration of patient lung 4DCT with subanatomical region segmentations. Med. Phys. 2008, 35, 775-781.

30. Barnafi, N.A.; Gatica, G.N.; Hurtado, D.E. Primal and Mixed Finite Element Methods for Deformable Image Registration Problems. SIAM J. Imaging Sci. 2018, 11, 2529-2567.

31. Rueckert, D.; Schnabel, J.A. Medical Image Registration. In Biomedical Image Processing; Deserno, T.M., Ed.; Springer: Berlin/Heidelberg, Germany, 2011; pp. 131-154.

32. Schmidt-Richberg, A. Registration Methods for Pulmonary Image Analysis; Springer Vieweg: Wiesbaden, Germany, 2013.

33. Wells, W.M.; Viola, P.; Atsumi, H.; Nakajima, S.; Kikinis, R. Multi-modal volume registration by maximization of mutual information. Med. Image Anal. 1996, 1, 35-51.

34. Lu, W.; Chen, M.L.; Olivera, G.H.; Ruchala, K.J.; Mackie, T.R. Fast free-form deformable registration via calculus of variations. Phys. Med. Biol. 2004, 49, 3067-3087.

35. Lubliner, J. Plasticity Theory; Dover Publications: Mineola, NY, USA, 2013; p. 544.

36. De Souza Neto, E.A.; Perić, D.; Owen, D.R.J. Computational Methods for Plasticity: Theory and Applications; John Wiley & Sons: Chichester, UK, 2008.

37. Radovitzky, R.; Ortiz, M. Error estimation and adaptive meshing in strongly nonlinear dynamic problems. Comput. Methods Appl. Mech. Eng. 1999, 172, 203-240.

38. Ortiz, M.; Stainier, L. The variational formulation of viscoplastic constitutive updates. Comput. Methods Appl. Mech. Eng. 1999, 171, 419-444.

39. Hurtado, D.E.; Ortiz, M. Finite element analysis of geometrically necessary dislocations in crystal plasticity. Int. J. Numer. Methods Eng. 2013, 93, 66-79.

40. Hurtado, D.E.; Henao, D. Gradient flows and variational principles for cardiac electrophysiology: Toward efficient and robust numerical simulations of the electrical activity of the heart. Comput. Methods Appl. Mech. Eng. 2014, 273, 238-254.

41. Modat, M.; Mcclelland, J.; Ourselin, S. Lung registration using the NiftyReg package. In Proceedings of the 13th International Conference on Medical Image Computing and Computer Assisted Intervention, MICCAI2010 Workshop: Medical Image Analysis for the Clinic-A Grand Challenge 2010, Beijing, China, 20-24 September 2010; pp. 33-42.

42. Cao, K.; Christensen, G.E.; Ding, K.; Du, K.; Raghavan, M.L.; Amelon, R.E.; Baker, K.M.; Hoffman, E.A.; Reinhardt, J.M. Tracking regional tissue volume and function change in lung using image registration. Int. J. Biomed. Imaging 2012, 2012, 956248.

43. Genet, M.; Stoeck, C.T.; Von Deuster, C.; Lee, L.C.; Kozerke, S. Equilibrated warping: Finite element image registration with finite strain equilibrium gap regularization. Med. Image Anal. 2018, 50, 1-22.