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

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.


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

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 Ω ⊂ R n be a domain of interest (image support), R : Ω → R be the reference image and T : Ω → R be the target image. The DIR problem aims to establish an optimal transformation u : Ω → R n that best aligns the reference and target images. To this end, we consider the functional space V := H 1 (Ω, R n ) and define a similarity functional D : V → R that penalizes differences between the reference image R and the resampled target image T • (id + u). A popular choice for the similarity measure in mono-modal applications of DIR [31,32] is the sum of squared-differences 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 regularizer S : V → R that 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 where the elastic energy density takes the form with ∇ the gradient operator, div the divergence operator, and λ and µ the Lamé constants.
With these elements, we define the elastic registration functional as where α > 0 is a weighting parameter. Then, the optimal transformation u is the minimizer of the elastic registration functional, and the DIR problem is formulated as the following variational problem: Find u such that We note that the optimal transformation u 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 transformation u, we define the strain tensor operator and we note that the elastic energy density defined in (3) can be rewritten as 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 where I is the identity tensor.

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 where ε e corresponds 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 thaṫ where() indicates partial derivatives with respect to time. To reflect the path-dependent nature of inelastic solids, we consider the effective inelastic strain q ∈ M := L 2 (Ω, 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 where we assume that the stored plastic energy takes the form W p = 1 2 Hq 2 , with H being the hardening modulus. Then, the elastic constitutive relation reads σ(ε e ) = ∂W ∂ε e = 2µε e + λ trace (ε e )I, (12) and the relation between the hardening internal variable and its thermodynamic conjugate stress is 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 where s ij are the components of the deviatoric stress tensor s defined as Following an associative plasticity framework [35], the evolution of the inelastic strain tensor is governed by the flow rulė 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 Finally, an inelastic model must comply with the loading/unloading complimen- 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ε p and q is not independent. Letσ be the effective stress defined asσ Then, the relation betweenε p and q is assumed to follow the Prandtl-Reuss flow rule that readsε 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, (see Appendix A.1 for details).

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 returnmapping algorithm based on an implicit backward-Euler temporal discretization. To this end, the time variable is discretized in generic subintervals [t n , t n+1 ]. Then, a series of incremental problems are obtained, where the main variables of the inelasticity model are assumed to be known at the time t = t n , and need to be solved for t = t n+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 Solving for ε p n+1 from the non-linear Equation (22) delivers an incremental update for the inelastic strain tensor, which we express as which depends solely of ε n+1 and q n+1 . Based on this flow-rule update, we define the effective incremental energy density for t = t n as [38], where, where ∆t = t n+1 − t n , 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 = q n+1 − q n . The minimization problem involved in the definition (25) is equivalent to the stationary condition which, for the rate-independent case reads (see details in Appendix A.2), whereσ pre n+1 is the elastic predictor forσ n+1 . Substituting (13) into (27) we obtain where the sub-differential of ψ * is defined as: 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 field u n is assumed to be known, and we find the displacement field u n+1 by solving the problem where the inelastic DIR functional reads and the inelastic regularizing term takes the form To solve the minimization problem (30) we consider the stationary condition The residual in (33) takes the form where 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 where is the consistent tangent tensor, see Appendix A.4. Thus, the linearized version of the residual problem reads: Given an initial guess w ∈ V, find the increment ∆w such that 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 where {N 1 , . . . , N m } 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 guess u h , find the increment ∆u h such that Using standard arguments (e.g., see [40]) we can show that (40) is equivalent to solving the linear system of equations where ∆u is a vector with the nodal values of the increment ∆u h , K n is the tangent matrix and R n 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 t n+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.

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 of 80 × 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: In addition, the normalized target registration error (TRE) was also computed. The TRE is defined as where p i (x, y) and q i (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 and images of the elastic von Mises strain, which takes the form We note that for purely elastic methods (FFD, Elastic FEM), the elastic strain tensor ε e is 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.

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: BE 1 = 1 × 10 −9 , BE 2 = 2.5 × 10 −6 and BE 3 = 1 × 10 −4 , 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.

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.    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 that 0.1 result in RSS errors that do no change considerably, delivering the highest accuracy observed for all three methods. In contrast, yield stress values above 1.0 deliver a much higher RSS error, which approaches that of the Elastic FEM method, see Table 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.

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 ≤ 10 −2 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.

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: Going a step back and clearingε p from (A1) we have,ε which in view of the plastic flow rule, is analogous tȯε Then by integration of (A5) and again, assuming isotropic hardening, we have the following relation,ε Finally, we compute the driving forces for q as,
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 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].