On Finite Difference Jacobian Computation in Deformable Image Registration

Producing spatial transformations that are diffeomorphic is a key goal in deformable image registration. As a diffeomorphic transformation should have positive Jacobian determinant |J| everywhere, the number of voxels with |J|<0 has been used to test for diffeomorphism and also to measure the irregularity of the transformation. For digital transformations, |J| is commonly approximated using a central difference, but this strategy can yield positive |J|'s for transformations that are clearly not diffeomorphic -- even at the voxel resolution level. To show this, we first investigate the geometric meaning of different finite difference approximations of |J|. We show that to determine if a deformation is diffeomorphic for digital images, the use of any individual finite difference approximation of |J| is insufficient. We further demonstrate that for a 2D transformation, four unique finite difference approximations of |J|'s must be positive to ensure that the entire domain is invertible and free of folding at the pixel level. For a 3D transformation, ten unique finite differences approximations of |J|'s are required to be positive. Our proposed digital diffeomorphism criteria solves several errors inherent in the central difference approximation of |J| and accurately detects non-diffeomorphic digital transformations. The source code of this work is available at https://github.com/yihao6/digital_diffeomorphism.


Introduction
The goal of deformable image registration is to establish a nonlinear spatial transformation that aligns two images.For many tasks, it is reasonable to assume that the anatomy in the two images share the same topology.Therefore, for deformable registration algorithms, the ability to produce a topology preserving transformation is preferred.Following the work of Christensen et al. (Christensen, 1994), many registration algorithms either penalize (Avants et al., 2008;Chen et al., 2017;Mok and Chung, 2020) or constrain (Beg et al., 2005;Chen et al., 2015;Dalca et al., 2018) their output transformations to be diffeomorphic and preserve topology.In the continuous domain, a diffeomorphic transformation is a smooth and invertible mapping with a smooth inverse that is guaranteed to maintain the topology of the anatomy being transformed.A diffeomorphic transformation should have positive Jacobian determinant |J| everywhere 1 , and testing if a transformation is diffeomorphic involves local computation of |J|.When a transformation is not diffeomorphic, it is common to use the number of voxels with negative |J| (Balakrishnan et al., 2019;Chen et al., 2022Chen et al., , 2021;;Liu et al., 2022) or the standard deviation of the logarithmic transformed |J| (Hering et al., 2022;Kabus et al., 2009;Johnson and Christensen, 2002) to measure the irregularity of the transformation.
Given a digital transformation that is defined on a regular grid, a widely accepted practice for computing the Jacobian is to use finite difference approximations of spatial derivatives.There are three standard methods for computing finite differences along each axis.We denote the forward, backward, and central differences that operate along the x axis (and similarly for the y and z axes) as D +x , D −x , and D 0x , respectively.To approximate |J| in 2D or 3D transformations, either the same or a different type of finite difference can be used for different axes.We denote the central difference approximation of |J| in 2D and 3D as D 0x D 0y |J| and D 0x D 0y D 0z |J|, respectively.
Despite its current popularity, the central difference approximation of the Jacobian does not always work as expected in the evaluation of the diffeomorphism property of spatial transformations.For example, Fig. 1(a) shows a transformation around center point p that is not diffeomorphic despite the fact that D 0x D 0y |J| = 1.In fact, the transformation at p has no effect on the computation of D 0x D 0y |J|(p), even if p moves outside the field of view.We call this the checkerboard problem because the central difference based Jacobian computations on the transformations of the "black" and "white" pixels (of a checkerboard) are independent of each other.A possible solution to this problem is to use D  In both figures, the center point p is transformed to p t .The transformation around the center point is visualized as displacements (shown as dotted arrows pointing toward solid dots) and the displacement for the center point is highlighted in red examples illustrate the problem with naive application of finite differences in this application.We seek a better approach that still involves finite differences, but avoids these types of contradictory situations.
In this work, we first investigate the geometric meaning of finite difference approximations of |J|.We show that when using forward or backward differences, the sign of |J|(p) determines if the underlying transformation T is invertible and orientation-preserving in a triangle (2D) or tetrahedron (3D) adjacent to p. Reversing the orientation indicates folding in space.We formally define digital transformations that are globally invertible and free of folding as digital diffeomorphisms.In order to determine if a transformation is a digital diffeomorphism, at each point it is necessary to consider four finite difference approximated |J|'s for 2D transformations and ten |J|'s for 3D transformations.We also demonstrate that because of the checkerboard problem and other errors that are inherent in the central difference based |J|, the number of non-diffeomorphic voxels it reports is always less than or equal to the actual number.Finally, we propose to use non-diffeomorphic area (2D) and volume (3D) as more meaningful measurements of irregularity in computed transformations.Consider the standard Euclidean space R 2 that follows the right-hand rule.Let T be a digital transformation for R 2 that is defined for every grid point p.When using backward differences on both x and y axes for approximating |J| at point p, we have the formulation where T x (p) and T y (p) are the x and y components of T (p).We denote the 4-connected neighbors of p as p −x , p +x , p −y , and p +y , as shown in Fig. 2(a).Their transformed locations are denoted with subscripts t, as shown in Fig. 2(b).For example, p +x t := T (p +x ).We denote the triangular region defined by the vectors pp −x and pp −y as △pp −x p −y , and we assume that the 2D transformation T is linearly interpolated on △pp −x p −y .
Proof Since T is linearly interpolated on △pp −x p −y , T is linear for △pp −x p −y and can be written as a 2 × 2 matrix with p t p −x t and p t p −y t as its columns.Thus, T is invertible if and only if p t p −x t and p t p −y t are linearly independent (not colinear).D −x D −y |J|(p) can be written as a triple product: where n is the unit vector perpendicular to vectors p t p −x t and p t p −y t and is positively oriented following the right-hand rule.Therefore, T is invertible if and only if Proof By the right-hand rule, △pp −x p −y is positively oriented.(⇒) When T is free of folding, the orientation of △pp −x p −y is preserved by T .Because of linear interpolation, △pp −x p −y is transformed to △p t p −x t p −y t , which is also positively oriented.Equation 2 shows that D −x D −y |J|(p) equals twice the signed area of △p t p −x t p −y t .Therefore, D −x D −y |J|(p) > 0.
(⇐) When D −x D −y |J|(p) > 0, from Eq. 2 △p t p −x t p −y t is positively oriented.Because of linear interpolation, △pp −x p −y is transformed to △p t p −x t p −y t and both of them are positively oriented.Therefore, T is free of folding for △pp −x p −y .
Definition 2 A 2D transformation T is digitally diffeomorphic for the region △pp −x p −y if T is invertible and free of folding for △pp −x p −y .
Proposition 3 A 2D transformation T is digitally diffeomorphic for the region △pp −x p −y if and only if Proof Proposition 3 is a direct consequence of Propositions 1 and 2.
In conclusion, D −x D −y |J| for the point p informs us about the digitally diffeomorphic property of a triangle adjacent to p, under the assumption that the transformation is linearly interpolated.

Digital Diffeomorphism in Two Dimensions
Similar to D −x D −y |J|(p), we can replicate Proposition 3 to establish that any |J|(p) approximated using any combination of forward and backward differences is testing if the transformation is digitally diffeomorphic for a triangular region around p (see Fig. Since there are two ways of dividing the square-size area between grid points, there are many more piecewise linear transformations that correspond to the same digital transformation, e.g., Fig. 3(c).To anticipate all possible choices of triangulation-each of which leads to its own invertible transformation on the plane-we should therefore consider all finite difference approximations in determining whether a transformation is digitally diffeomorphic.This leads naturally to the following definition: The proposed digital diffeomorphism definition guarantees the transformation to be free of folding and invertible regardless of the piecewise linear transformation (on triangles that divide the squares between pixel centers) that is used.

Central Difference Based Jacobian Determinant
In this section, we analyze the central difference approximation of |J| given the previous analysis.
We first ask where p can be positioned to yield a digital diffeomorphism when its neighbors have fixed transformations.We start with the following definitions: Definition 4 For grid point p with fixed p −x t , p −y t , p +x t , and Definition 5 Let a and b be points in R 2 .The halfplane H(a, b) is defined as p ∈ R 2 such that △abp is positively oriented.
Proof Equation 1 can be written as:   Proof D 0x D 0y |J|(p) can be written as: (⇐) On the other hand, if D 0x D 0y |J| > 0, the polygon p −x t p −y t p +x t p +y t is positively oriented (Braden, 1986) (i.e.interior to the left).Thus, if there exists a p t that is visible to all vertices of the polygon, then p t ∈ R(p) by Proposition 4. Following (Chvátal, 1975), such p t always exists for simple polygons with five or fewer vertices.Therefore, R(p) is non-empty.
Although a positive D 0x D 0y |J|(p) ensures that R(p) is non-empty, p can still be digitally nondiffeomorphic when p t ̸ ∈ R(p), which we see in the checkerboard problem.In addition to the checkerboard problem, D 0x D 0y |J| also fails to provide a meaningful interpretation when the polygon p −x t p −y t p +x t p +y t exhibits self-intersection (see Fig. 4(d)).In conclusion, D 0x D 0y |J|(p) ≤ 0 always indicates that T is digitally non-diffeomorphic but D 0x D 0y |J|(p) > 0 does not mean it is digitally diffeomorphic because of the checkerboard or selfintersection problems.Therefore, use of central differences alone to estimate |J| is an inadequate characterization of digital diffeomorphism.

Digital Diffeomorphism in Three Dimensions
Consider the standard Euclidean space R 3 that follows the right-hand rule.Let T be a digital transformation of R 3 that is defined for every grid point p.We denote the 6-connected neighbors of p as p ±x , p ±y , p ±z (see Fig. 5(a)) and their transformed locations are denoted with the subscript t.
We denote the tetrahedron defined by the vectors pp −x , pp −y , and pp −z as pp −x p −y p −z , and we assume that the 3D transformation T is linearly interpolated on pp −x p −y p −z .
Proposition 6 A 3D transformation T is invertible for pp −x p −y p −z if and only if Proof Since T is linearly interpolated on pp −x p −y p −z , T is linear for pp −x p −y p −z and can be written as a 3 × 3 matrix with p t p −x t , p t p −y t , and p t p −z t as its columns.Thus, T is invertible if and only if p t p −x t , p t p −y t , and p t p −z t are linearly independent (not colinear).D −x D −y D −z |J|(p) can be written as a triple product: Proof Proposition 8 is a direct consequence of Propositions 6 and 7.
Similarly, we can prove that a |J|(p) approximated using any combination of forward and backward differences is testing if the transformation is digitally diffeomorphic in a tetrahedron adjacent to p.One of these tetrahedra (for D −x D +y D −z |J|) is visualized in Fig. 5(a).However, unlike the 2D case where all the triangular regions completely cover the entire 2D space, the union of all these adjacent tetrahedra does not fill the entire 3D space.Therefore, even if all forward and backward difference approximations of |J| are positive, the transformation can still cause folding in the spaces not covered by these adjacent tetrahedra.
To solve this issue, we introduce another tetrahedron pp −x−y p −x−z p −y−z , as shown in Fig. 5 There are other ways to divide a cube into tetrahedra (Carr et al., 2006).These other schemes are not considered here because 1) their computation involve finite difference approximations at different points than those we consider, or 2) their computation requires interpolating the transformation at non-grid point.For similar reason as in 2D, our definition of digital diffeomorphism involves both schemes shown in Figs.5(b) and (c).Specifically, each of the two schemes corresponds to a particular piecewise linear transformation.But for a given digital transformation there are many plausible piecewise linear transformations (see Fig. 3(c)).By considering both schemes, we avoid the ambiguity of choosing different schemes for every cube-sized volume. The

Non-diffeomorphic Space Measurement
For non-diffeomorphic transformations in 3D, the number and percentage of non-diffeomorphic voxels is often used to measure its irregularity.Specifically, a voxel is considered non-diffeomorphic if its center location p has a central difference based Jacobian determinant that is not positive (i.e. in 3D: D 0x D 0y D 0z |J|(p) ≤ 0).However, as we have demonstrated in Sections 1 and 2, the central difference based Jacobian determinant underestimates the non-diffeomorphic space, in general.Given the limitations of the central difference-based Jacobian determinant, we seek an alternative way of evaluating the diffeomorphic property of digital transformations.We aim to find a quantitative measure that involves the computation of finite difference-based Jacobians and is consistent with our definition of digital diffeomorphism.As a result, we propose non-diffeomorphic volume (NDV) to quantify the size of the non-diffeomorphic space in 3D caused by a digital transformation.While it is possible to report the minimum or maximum achievable non-diffeomorphic volume by tetrahedralizing each cube between voxels based on the specific transformation, doing so would inevitably associate the given digital transformation with one particular piecewise linear transformation.Therefore, averaging the two tetrahedralization schemes provides a more comprehensive and unbiased measure for a digital transformation.Further discussion on the impact of spatially varying tetrahedralization schemes can be found in Section 4. The proposed NDV is connected to the ideas of simplex counting (SC) (Yushkevich et al., 2010) and surface propagation (SP) (Pai et al., 2016), which are used to assess the degree of volume change on a per-voxel basis.However, those methods typically use only one scheme to discretize the space.In contrast, NDV considers the two tetrahedralization schemes shown in Figs.5(b)  and (c).This choice is motivated by the definition of digital diffeomorphism but it also offers a rotation-invariant property.Specifically, rotating the transformation by any multiple of 90 degrees does not affect its result.This cannot be achieved using only one scheme.Further discussion on the relationship and distinction between SC, SP, and the proposed NDV can be found in Section 4.
For completeness we note the 2D version of non-diffeomorphic space, which we term nondiffeomorphic area (NDA), follows from Def. 3, where . . , 4], are the Jacobian determinants approximated using the four possible combinations of forward and backward differences.
A demonstration of NDA is provided in Fig. 6.When using the central difference based |J|, all three cases in Fig. 6 would be considered diffeomorphic (D 0x D 0y D 0z |J|(p) > 0) because of the checkerboard problem.The forward difference based |J| is able to identify that Figs.6(b) and (c) exhibit folding, but only NDA can provide the observation that the non-diffeomorphic space caused by the transformation shown in Fig. 6(c) is larger than the non-diffeomorphic space in Fig. 6(b).
Table 1 Our proposed non-diffeomorphic volume and several other measures on the IXI dataset and the seven comparison algorithms.For 'D 0x D 0y D 0z |J| ≤ 0' and 'Any |J i | ≤ 0' we report the mean number of voxels (voxel #) over the 115 test subjects and the corresponding standard deviation (±), as well as the percentage (%) with respect to the brain mask.We also report our proposed measure of non-diffeomorphic space-i.e.non-diffeomorphic volume (NDV) in 3D-and the corresponding standard deviations and percentages.Methods are listed in the order in which they were published.

Experiments
We compared the commonly used central difference based Jacobian determinant (D 0x D 0y D 0z |J|) and our proposed non-diffeomorphic volume (NDV) using several deformable registration algorithms on two publicly available datasets: IXI: A total of 576 T1-weighted brain magnetic resonance (MR) images from the publicly available IXI dataset were used.403 scans were used in training for the task of atlas-to-subject registration (Kim et al., 2021) and 58 scans were used for validation.The transformations generated from registering an atlas brain MR images to 115 test scans were evaluated.
Learn2Reg OASIS: We also used the brain T1weighted MR images from the 2021 Learn2Reg challenge (Hering et al., 2022;LaMontagne et al., 2019).Scans were preprocessed using FreeSurfer (Fischl, 2012;Hoopes et al., 2021).All algorithms were trained using the training set of 414 scans and the transformations for the 19 validation pairs were evaluated.
The central difference Jacobian determinant approximation was implemented directly from the 2021 Learn2Reg challenge evaluation script.The implementation details and hyper-parameters for each of the algorithms were adopted from (Chen et al., 2022) and (Liu et al., 2022).
A visualization of a result from the IXI dataset is shown in Fig. 7. Since it is difficult to visualize displacements across slices (anterior-to posterior direction in this case), only the displacements within the coronal plane were considered.Thus, we computed the non-diffeomorphic area instead of non-diffeomorphic volume.For each pixel in Fig. 7(d), a higher red intensity indicates larger non-diffeomorphic area around that pixel.The nondiffeomorphic pixels computed from D 0x D 0y |J| are highlighted in Fig. 7(c) for comparison.
The results on the IXI dataset are summarized in Table 1 and for the Learn2Reg OASIS in Table 2.Only the voxels within the brain were considered and the percentages were calculated relative to the brain volume of the fixed image.In both tables, we report the number of non-diffeomorphic voxels (# of voxel) and its percentage (%) based on the D 0x D 0y D 0z |J|.We also report the number (and percentage) of voxels that have at least one |J i | ≤ 0, denoted as 'Any |J i | ≤ 0' in the tables.We observe that in most of the cases, there are actually more than twice the number of voxels having |J i | ≤ 0 for some finite difference than found using only the central difference.The differences between 'D 0x D 0y D 0z |J| ≤ 0' and 'Any |J i | ≤ 0' highlight that errors in using the central difference approximation are very common in practice.
The proposed average non-diffeomorphic volume (NDV) and its percentage are shown in the last three columns of Tables 1 and 2. For algorithms that impose strong regularization on the transformations (e.g., MIDIR (Qiu et al., 2021), deedsBCV (Heinrich et al., 2015), NiftyReg (Modat et al., 2010), andSyN (Avants et al., 2008)), a small-even zero-NDV is observed.However, for the deep learning methods that directly output deformation fields (Kim et al., 2021;Balakrishnan et al., 2019;Chen et al., 2022;Liu et al., 2022), we usually have higher NDVs.It is important to note that the results shown in Tables 1  and 2 do not reflect the accuracy of the algorithms, just the proportion of their deformation that is non-diffeomorphic.

Discussion
The Jacobian determinant of a transformation is a widely used measure in deformable image registration, but the details of its computation are often overlooked.In this paper, we focused on the finite difference based approximation of |J|.Contrary to what one might expect, the commonly used central difference based |J| does not reflect if the transformation is diffeomorphic or not.Our investigation shows that each of the finite difference approximated |J| corresponds to the signed area of a triangle in 2D or the signed volume of a tetrahedron in 3D when the digital transformations are assumed to be piecewise linear.Following this, we propose the definition of a digital diffeomorphism that allows diffeomorphisms-a concept in continuous domain-to be applied to digital transformations.It solves several problems that are inherent in the central difference based |J|.We further propose to use non-diffeomorphic volume to measure the irregularity of 3D transformations and non-diffeomorphic area for 2D transformations.As demonstrated in Fig. 6 and Fig. 7, our proposed approach measures the severity of the irregularity whereas the commonly used central difference based |J| is only a binary indicator of folding (also with errors).The transformation shown in Fig. 6(b) is obviously more favorable than the one shown in Fig. 6(c) in terms of regularity.As such, it is important for us to be able to draw distinctions between these two scenarios.
The non-diffeomorphic measures presented in Eqs.7 and 8 are averages of two choices of tetrahedralization of the volume.It is possible, however, to tetrahedralize each cube between voxels based on the specific registration outcome, for example, to yield the minimum or maximum achievable nondiffeomorphic volume.These measures must be computed for the entire image domain since the brain mask is defined at voxel level.For the IXI dataset, the average non-diffeomorphic volume for the best case is 33227 voxel 3 for Voxelmorph and 31488 voxel 3 for Transmorph while the average non-diffeomorphic Volume for the worse case is 41903 voxel 3 for Voxelmorph and 40964 voxel 3 for Transmorph.In comparison, our proposed NDV is 37565 voxel 3 for Voxelmorph and 36226 voxel 3 for Transmorph for the entire image domain.This result shows that, at least for these two algorithms, the non-diffeomorphic volume cannot be made substantially different from the average value that we specified in Eq. 7 by alternative selection of tetrahedralization.
The idea of discretizing the space using triangles or tetrahedrons has been presented before for regularizing deformable transformations and deformation-based volume change estimation (Pai et al., 2016;Holland et al., 2011;Yushkevich et al., 2010).It was previously believed that computing the volume change in discretized space using tetrahedrons is more accurate than using Jacobian determinants because the latter involves finite difference approximation (Yushkevich et al., 2010).Our analysis shows that discretizing the space using triangles or tetrahedrons is in fact the consequence of finite difference approximation of Jacobian determinants.Therefore, the two approaches are equivalent when the corresponding finite differences are chosen for a given partition of space.Haber et al. (Haber andModersitzki, 2004, 2007) and Burger et al. (Burger et al., 2013) used the volume of triangles in 2D or tetrahedrons in 3D as a regularization term for their registration algorithm.Initially, Haber et al. proposed a hard equality constraint (Haber and Modersitzki, 2004) that enforces the preservation of the discretized volume (or area in 2D) of every deformed box.Recognizing that this approach could not detect Table 2 Results for the Learn2Reg OASIS dataset.For 'D 0x D 0y D 0z |J| ≤ 0' and 'Any |J i | ≤ 0' we report the mean number of voxels (voxel #) over the 19 validation subject pairs and the corresponding standard deviation (±), as well as the percentage (%) with respect to the brain mask.We also report our proposed measure of non-diffeomorphic space-i.e.non-diffeomorphic volume (NDV) in 3D-and the corresponding standard deviations and percentages.Methods are listed in the order in which they were published."twists" (i.e., folding), they later proposed an inequality constraint (Haber and Modersitzki, 2007) that calculates the volumes of the tetrahedrons to prevent such twists by imposing positive volumes of the tetrahedrons.However, their methods only account for a single combination of Jacobian determinants (as shown in Fig. 3(a)), which is insufficient to guarantee a digital diffeomorphism.Moreover, these previous works were motivated by the fact that calculating the volume of a "twisted" polygon or octahedron is difficult.Our analysis on the central difference approximated |J| explains why using a polygon or octahedron is inaccurate.With respect to deep network based registration methods, many ensure that the output of their network is diffeomorphic, by adopting a scalingand-squaring layer as the output layer (Dalca et al., 2019;Chen et al., 2022;Hoopes et al., 2021).However, several works have reported that the scaling-and-squaring approach produces voxels with negative central difference based Jacobian determinant.The analysis presented in this paper explains why the scaling-and-squaring approach cannot guarantee a folding-free digital transformation.Specifically, when composing two digital transformations, one of the transformation needs to be sampled at non-grid locations, which usually involves bilinear or trilinear interpolation.These interpolation methods, however, are inconsistent with the piecewise linear transformation that is implicitly assumed by the finite difference based Jacobian determinant computation.As a result, the sampling process can introduce folding and result in locations with negative Jacobian determinant.
For recent deep learning based registration methods, our digital diffeomorphism criteria has the potential to be used as a loss function to improve the smoothness and promote digitally diffeomorphic transformations (Mok and Chung, 2020).This is a promising direction for future research.

Statements and Declarations
Data availability The datasets analysed during the current study are available in the OA-SIS repository, https://www.oasis-brains.org/ and the IXI repository, https://brain-development.org/ ixi-dataset/ Ethical approval declarations This work involved human subjects or animals in its research.Approval of all ethical and experimental procedures and protocols was granted by the relevant local Institutional Review Boards, and performed in line with the Declaration of Helsinki.Conflict of interest The authors have no competing interests to declare that are relevant to the content of this article −x D −y |J| or D +x D +y |J| instead.However, Fig. 1(b) shows an example in which D −x D −y |J| and D +x D +y |J| have opposite signs, leading to contradictory conclusions.These 1 We do not consider the case of all negative |J|, which would produce a reflection of the entire image.

Figure 1
Figure 1 (a) is an illustration of the checkerboard problem for D 0x D 0y |J|.(b) is a transformation that illustrates the inconsistency issue between D −x D −y |J| and D +x D +y |J|.In both figures, the center point p is transformed to p t .The transformation around the center point is visualized as displacements (shown as dotted arrows pointing toward solid dots) and the displacement for the center point is highlighted in red

D
Figure 2 (a) shows the notations for p and its 4-connected neighbors.(b) shows p t , the transformed version of p as well as the transformed position of the 4-connected neighbors.(c) shows the triangular regions and their corresponding forward and backward difference based |J|'s 2(c)), assuming that the region is linearly interpolated.Since these |J|(p)'s cover different regions, their signs are independent of each other.For example, T can have a positive D −x D −y |J|(p) and a negative D +x D +y |J|(p) at the same time (see Fig.1(b)).For each of these approximations, the implied triangular regions for all p's taken together cover only half of the space (e.g., D −x D −y |J| only considers the red tiles in Fig.3(a)).Consequently, even if a transformation T has D −x D −y |J|(p) > 0 for all p's, half the space is not considered and can potentially exhibit folding or be non-invertible.This is also the case for any other forward and backward difference computations of |J|.Although combining the triangular regions of D −x D −y |J| and D +x D +y |J| can cover the entire space, positive D −x D −y |J| and D +x D +y |J| only guarantee that the transformation is digitally diffeomorphic when the transformation is piecewise linearly interpolated as shown in Fig.3(a).A different choice, for example that shown in Fig.3(b), which corresponds to D −x D +y |J| and D +x D −y |J| can give contradictory conclusions.

Figure 4 Figure 3
Figure4shows examples of R(p).In Figs.4(c) and (d), no matter where p t is located, at least one forward or backward difference based |J|(p) is non-empty, for every p t ∈ R(p) all its forward and backward difference based |J|(p)'s are positive by Definition 4, and therefore from Eq. 3, D 0x D 0y |J|(p) > 0.
4) Therefore, T is invertible if and only if D −x D −y D −z |J|(p) ̸ = 0. Definition 6 A 3D transformation T is said to cause folding for pp −x p −y p −z if the orientation of pp −x p −y p −z is reversed by T .Proposition 7 A 3D transformation T is free of folding for pp −x p −y p −z if and only if T has D −x D −y D −z |J|(p) > 0. Proof (⇒) When T is free of folding, the orientation of pp −x p −y p −z (negatively oriented by the right-hand rule) is preserved by T .Because of linear interpolation, pp −x p −y p −z is transformed to p t p −x t p −y t p −z t , which is also negatively oriented.Equation 4 shows that D −x D −y D −z |J|(p) equals to six times the negative signed volume of p t p −x t p −y t p −z t .Therefore, D −x D −y D −z |J|(p) > 0. (⇐) D −x D −y D −z |J|(p) > 0 indicates that p t p −x t p −y t p −z t is negatively oriented by Eq. 4. Because of linear interpolation pp −x p −y p −z is mapped to p t p −x t p −y t p −z t and both of them are negatively oriented.Therefore, T is free of folding for pp −x p −y p −z .Definition 7 A 3D transformation T is digitally diffeomorphic for the region pp −x p −y p −z if T is invertible and free of folding for pp −x p −y p −z .Proposition 8 A 3D transformation T is digitally diffeomorphic for the region pp −x p −y p −z if and only if D −x D −y D −z |J|(p) > 0.

,
(b).When combined with four existing tetrahedra from finite difference based |J|'s, they completely cover the entire volume.We define |J ⋆ 1 |(p) as |J ⋆ 1 |(p) = (p t p −x−y t × p t p −x−z t ) • p t p −y−z t formed locations of p −x−y , p −x−z , p −y−z .The signed volume of the extra tetrahedron after applying T equals 1 6 |J ⋆ 1 |(p).Similar to 2D, there are two ways of dividing the cube-size volume inbetween grid points into five tetrahedra, as shown in Figs.5(b) and (c).The signed volume of the extra tetrahedron in Fig. 5(c) can be computed as 1 6 |J ⋆ 2 |(p), where

Figure 5
Figure 5 (a) shows the notations for p and its 6-connected neighbors in 3D and the tetrahedron considered by D −x D +y D −z |J|.(b) and (c) are illustrations of the two schemes to divide the cube volume in-between grid points in 3D

Figure 6
Figure6A demonstration of measuring non-diffeomorphic space in 2D.The transformations are visualized as displacement fields.In each of the three subfigures, only the center points are transformed to their corresponding p t 's and the other grid points remain fixed.The triangular region corresponds to the forward difference based |J| is shown in green if |J| > 0, otherwise it is shown in red

Figure 7
Figure 7 A visualization of the proposed non-diffeomorphic area (only the displacement inside the coronal plane was considered in this example).(a) The warped image.(b) The grid line representation of the transformation (generated using Voxelmorph (Balakrishnan et al., 2019)).(c) The warped image with non-diffeomorphic pixels (marked in red) as measured by the central difference Jacobian determinant (D 0x D 0y |J|) highlighted in red.(d) The warped image with a map overlay indicating the non-diffeomorphic area with brighter shades of red indicating larger non-diffeomorphic area Analogous statements can be made for the other forward and backward difference based |J|'s.When all the four |J|'s are positive, p t must be inside the four half-planes H(p −x D 0x D 0y D 0z |J| ≤ 0 Any |J i | ≤ 0 Proposed