3D DIFFEOMORPHIC IMAGE REGISTRATION WITH CAUCHY–RIEMANN CONSTRAINT AND LOWER BOUNDED DEFORMATION DIVERGENCE

. In order to eliminate mesh folding in 3D image registration problem, we propose a 3D diffeomorphic image registration model with Cauchy–Riemann constraint and lower bounded deformation divergence. This model preserves the local shape and ensures no mesh folding. The existence of solution for the proposed model is proved. Furthermore, an alternating directional projection 3D image registration algorithm is presented to solve the proposed model. Moreover, numerical tests show that the proposed algorithm is competitive compared with the other four algorithms.


Introduction
Image registration is spatial position calibration of two given data sets.It is an important tool for medical image analysis, for example, atlas-based image segmentation and shape analysis.Moreover, image registration is also widely used in remote sensing, astronomy and some other fields in industry.An overview of image registration problem can be found in [1][2][3].Here we focus on mono-modality 3D nonrigid deformable image registration problem in variational framework.Mathematically, 3D image registration problem can be stated in the following way.Let Ω ⊂ R 3 be an open bounded domain on R 3 , and let ,  : Ω → R be two real functions. and  are viewed as two 3D images, for example,  and  are tomography data.In image registration,  and  are named the floating image and the target image, respectively.The registration of  and  is to find a 3D spatial transformation h : Ω → Ω such that the intensity of  ∘ h(•) looks like (•) as much as possible.For example, finding an h such that the sum of squared difference (SSD) ∫︀ Ω | ∘ h(x) − (x)| 2 dx reaches its minimum.Addressing this issue, many variational models [4][5][6][7][58][59][60] are proposed.However, mesh folding is not taken into consideration in most of these models.In physical view, mesh folding implies the volumetric vanishing which contradicts physical principle.As a consequence, eliminating mesh folding is an essential technique for 3D image registration.To eliminate mesh folding for 3D image registration, deformation h : Ω → Ω should be a local inverse function.In diffeomorphic image registration framework, one should ensure that the Jacobian determinant of h is not equal to 0 (see Appendix C in [8] for details).Moreover, Jacobian determinant of h denotes the volume stretching rate for the transformation.To make deformation physically meaningful and preserve the orientation, one should ensure that Jacobian determinant of h is greater than 0 for any x ∈ Ω.For the purpose of eliminating mesh folding, Beg et al. [9][10][11][12][13][14][15][16][17] developed a large deformation diffeomorphic metric mapping (LDDMM) framework.In the LDDMM framework, image registration is formulated as a variational problem by viewing the solution as a geodesic flow in space of diffeomorphism.In the LDDMM framework, Dupuis [12] proved that no mesh folding occurred by restricting velocity field to the Sobolev space  3,2  0 (Ω).Moreover, several numerical tests also showed the efficiency of LDDMM.One can refer to [14][15][16] for details.In addition, by giving the hypothesis that h(x) is approximately divided into an identity part x and a displacement part u(x) (i.e., h(x) = x + u(x)), the 3D image registration problem can be also formulated by the following framework: where (u) is regularization (see [4][5][6][7]) and  is an admissible set of solutions.Here and in what follows, Many canonical registration models arise from framework (1), for example, Demons model [18], Diffeo Demons model [19], Vectorical total variation (VTV) model [20,21] and Bounded deformation (BD) model [6].However, mesh folding is not taken into consideration in framework (1).Therefore, searching for a framework (1) based image registration model without mesh folding is a challenging task.To solve this problem, a 3D orientationpreserving variational model is proposed by Zhang and Chen [22] based on framework (1).There are also some other relevant works on 2D quasi-conformal or conformal mapping (see [23][24][25][26][27][28][29][30][31][32] for details).Particularly, in [31,33], authors of this article proposed a 2D diffeomorphic image registration model by constraining the solution into a 2D conformal set.As an extension, to eliminate mesh folding in 3D image registration, we propose a 3D diffeomorphic image registration model with Cauchy-Riemann constraint and lower bounded deformation divergence (see Sect. 2 for details).This model preserves the local shape and ensures no mesh folding.Compared with 2D diffeomorphic image registration model in [31,33] which is only constrained by two equalities, three inequalities are additionally added in the proposed model to ensure that Jacobian determinant of h is greater than 0 for any x ∈ Ω.In this view, the proposed model is much more technically challenging compared with the model in [31,33].Around the proposed model, a relaxed model is proposed for the convenience of numerical implementation.The existence of solution for the relaxed model is proved.Moreover, an alternating directional projection 3D image registration algorithm is proposed.Numerical tests show that the proposed algorithm is competitive compared with the other four algorithms.
The rest of this paper is organized as follows: In Section 2, a 3D diffeomorphic image registration model without mesh folding is proposed.Based on Section 2, a relaxed model is also proposed and the existence of its solution is proved in Section 3. In addition, an alternating directional projection algorithm is proposed to solve the proposed model in Section 3. In Section 4, several numerical tests are performed to show the efficiency of the proposed projection algorithm.In Section 5, we summarize our work and list some problems for future research.

Proposed 3D diffeomorphic image registration model
At the beginning of this section, let's recall the following canonical theorem in Appendix C of [8].
) with lim ‖x‖→+∞ ‖h(x)‖ = +∞.It follows from the corollary on page 200 in [34] that h : R 3 → R 3 is a  1 global diffeomorphism.That is, the local bijection in Theorem 2.1 can be replaced by global bijection if the condition h(x) = x on Ω is added.
By Theorem 2.1, we know that constraint det(∇h(x)) > 0 is necessary to eliminate mesh folding in 3D image registration models.To achieve det(∇h(x)) > 0 for any x ∈ Ω, one can find several recent works which impose this constraint in direct way [35,36].However, in 3D space, Directly controlling det(∇h(x)) > 0 seems to be a too complex inequality constraint for 3D image registration problem.It is also not convenient for numerical implementation.To simplify 3D image registration with direct inequality constraint det(∇h(x)) > 0, we present a novel 3D diffeomorphic image registration model in this section.Before giving the proposed model, we define Remark 2.2.There are three remarks on the set : Note that here and in what follows, ∠ denotes the included angle of edge  and .In elastic mechanics,   ,   ,   are named shearing strain which represents the angle changes caused by h.In addition, there holds By restricting u into the set ,   =   =   = 0.That is,  produces diffeomorphic mappings h which locally preserve the angle.Compared with the conformal set which preserves the angle and stretching rate at the same time,  is a much more relaxed set.(ii) Define   =   ( = 1, 2, 3), then   are named strain on   direction.In physical view, 1 +   denotes the stretching rate of  1  2  3  4 −  1  2  3  4 on   direction.This is the main reason why the inequality constraint (x)   > −1 is given.In addition, by giving the equality constraint 1(x) 1 = 2(x) 2 = 3(x) 3 , one can ensure that the local stretching rate is the same on three different directions, which preserves the local shape.Moreover, by restricting u into the set , one can notice that div(u) = ∑︀ 3 Based on these definitions, the proposed 3D diffeomorphic image registration model is formulated by: where  > 2.5,  > 0, and  1 (u) = ∫︀ Ω |∇  u| 2 dx.Note that throughout this paper, Ω = (0,  1 ) × (0,  2 ) × (0,  3 ), for x = ( 1 ,  2 ,  3 ) ∈ Ω,  = 1, 2, 3 and function  : Ω → R, define and the details of these two fractional-order derivatives can be founded in [9,38].
Remark 2.3.There are four remarks on model ( 5): (i) Here we use the fractional-order regularization ∫︀ Ω |∇  u| 2 dx because the differential order of the fractionalorder operator can change continuously from zero to some positive number, which makes it possible to adjust the differential order to give a more accurate and more general image registration model.[39], Thm.4.58).That is, deformation h : Ω → Ω is diffeomorphic.This ensures that the first-order derivatives in  are well defined.Therefore, model ( 5) can be extended into diffusion tensor imaging (DTI) (,  : Ω → SPD(3)) registration by replacing (u) with ∫︀ Ω ‖ ◇ h(x) − (x)‖ 2 dx.Note that SPD(3) is a set which contains all 3 × 3 Symmetric Positive Definite real matrixes, and  ◇ h(x) is defined by One can refer to [11,16] for details.
That is, no mesh folding occurs in proposed model (5).This is the main motivation for us to add constraint u ∈  in the proposed model ( 5).(iv) By letting 1 3 = 2 3 = 0 in , then  becomes a 3D conformal set   .In this view, equation ( 5) is a much more relaxed model compared with 3D conformal image registration model.

Numerical solution of proposed model
At the beginning of this section, we introduce some definitions.
Here and in what follows, we define and where Θ > 0 is a large number.
To transform proposed model ( 5) into a variational problem with less constraints, we make use of the following approximation procedure: where In fact, the constraint set in ( 13) should be Here we introduce an infinite small variable  > 0 to construct a closed set ℬ  which is necessary in the proof of Theorem 3.5.Moreover, the solution of ( 13) approximately satisfies 29) and (30) for details).Concerning the existence of solution for the relaxed model ( 13), we have the following results: By the fact (u  ) ≤ (0) = (0) < +∞; there exists a constant  1 such that On the other hand; by the Lemma 3.3 in [9], we know that for any u  ∈ [  0 (Ω)] 3 , there holds for some constant  = (, Ω) > 0. Therefore, for some constant M .By ( 17) and ( 18), we know the semi-norm defined by 39], Thm.4.58) that there exists a subsequence of {u  } which is still labeled by  and u ∈ [ 1 (Ω)] 3 such that Based on these conclusions, now we claim that First, by the lower semi continuity of norm (see D4 on page 639 of [8]) and (18), we obtain that Furthermore, by (19), we obtain that, because of the fact that  (x) is continuous on Ω∖∆  and  (x) is bounded on Ω.Note that here h , where ∆ T is also a zero measure set.For simplicity, here we still write ∆ T as ∆  .
That is, Moreover, it follows from (19) that
Next, we give some results on numerical implementation of problem (13).
Theorem 3.3.Let u ∈ ℬ  be the solution of (13), then u satisfies the following variational inequality where Set ( ) (u +  (w − u)), by the fact (0) ≤ ( )(0 ≤  ≤ 1), we know that Therefore, Note that for any  ∈   0 (Ω) implies,   (x) Sobolev space in [8,9,39]).That is, any function on   0 (Ω) satisfies the homogeneous boundary conditions.By giving homogeneous boundary conditions, the following two important properties of fractional-order derivatives are used when deriving ( 35): (i) Riemann-Liouville derivatives, Grunwald-Letnikov derivatives and Caputo derivatives are equivalent [38]; (ii) The following integration by parts formula holds Now, we discuss the problem of how to find numerical solution of variational inequality (32).Before we give our results, let's recall an important lemma [41] which is necessary in our proof.Lemma 3.4.Suppose  is a closed, convex set,  :  → R 3 is monotone, then u is the solution of variational inequality if and only if u solves the fixed point equation for an arbitrary positive constant .
For technical demand in later proof, we define ≥ −1 +  for any x ∈ Ω}, where  is a very large number.Note that the solution of ( 13) is bounded on [  0 (Ω] 3 (i.e., (u ), ℬ   contains all the solutions of (13) when  is large enough.Hence, throughout this paper, we do not distinguish ℬ   and ℬ  except for the following theorem.Moreover, we assume that  (•) and ∇ (•) are Lipschitz continuous on Ω with Lipschitz constant .
where  > 0 is a small number and  ℬ   (v) is projection of v on convex set ℬ   .
Proof.Firstly, we claim that Since  2 (u) and  3 (u) are linear operators on u, by using the integration by parts formula, one can easily check that On the other hand, Note that here we use the assumption that  (•) and ∇ (•) are Lipschitz continuous on Ω with Lipschitz constant .
It follows from ( 42) and ( 43) that By the Lemma 3.3 in [9], there exists a Note that here we use the assumption  ≥ 3CML to ensure   − 3  ≥ 0. This concludes the claim (40).On the other hand, ℬ   is a convex and closed set.By Lemma 3.4, we conclude that u ∈ ℬ   is a solution of (32) if and only if u ∈ ℬ   satisfies Algorithm 1. Projection algorithm for 3D Diffeomorphic image registration.
Motivated by Theorem 3.5 and Chapter 5 in [42], variational inequality (32) can be solved by iteration formula Based on this formula, we give the following 3D image registration algorithm (Algorithm 1).
Remark 3.6.The discretization of ℱ(u  ) can be found in Appendix A. The convergence of Algorithm 1 is out of the scope of this paper.Perhaps, we will give some results on this topic in forthcoming work.
Based on (49), the projection  ℬ (v) for any given v ∈ R 3 can be formulated as Because there is no coupled term of  ,,  ( = 1, 2, 3) in (50), equation ( 50) is equivalent to the following three independent subproblems: These are three convex optimization problems with inequality constraints.Note that here   is a block diagonal matrix.Motivated by Han [9,43], we hope to divide 3D problem (50) into several parallel 1D problems.For this purpose, equation ( 51) is divided into the following 3 parallel 1D subproblems: with constraints  1  , + 1, ≥ 0, where here and in what follows,  , is a  1 ×1 vector,   is identity matrix of order ,   , = −2( with constraints  2  , +  2, ≥ 0, where  , is a with constraints  3  , +  3, ≥ 0, where  , is a Remark 3.8.By ignoring the constant term in objective function, equations ( 52)-( 54) are equivalent to and Based on previous discussion, the algorithm to compute  ℬ (v) can be summarized as the following Algorithm 2.
Apply PPA to the problem (58) by letting Both ( 61) and ( 62) are unconstrained optimization problems, which implies the gradient of objective functions equal to 0 at optimal points.Based on this conclusion, the following solution can be easily obtained: where Based on previous discussion, PPA for the quadratic programming problem (58) can be summarized as the following Algorithm 3.

Algorithm 3. PPA for quadratic programming problem (58).
Input: , , , ,  ,  0 ,  0 Set  = 0, norm=1, accuracy=0.001 and  = 50.While norm> accuracy and  ≤  do: EndWhile Return: X =   0 for some 0 ≤ .Remark 3.9.One can also use MATLAB code "x=quadrog (H,f,-A,-b)" to solve the quadratic programming problem (58).The numerical techniques for quadrog function in MATLAB can be found in Section 3 of [48].However, quadrog function converges much slower than Algorithm 3.For example, considering the following quadratic programming problem: where Both MATLAB quadrog function and Algorithm 3 achieve the optimal solution X = (−0.51,−1.98, −2.49)  .However, it takes 0.321356 s for MATLAB quadrog function to converge to the optimal solution while Algorithm 3 costs only 0.018727 s.In fact, in numerical tests which will be introduced in Section 4, we have noticed that Algorithm 1 converges much slower when Algorithm 3 is replaced by quadrog function in MATLAB.

Numerical test
In this section, two numerical tests are performed to show the efficiency of the proposed Algorithm 1.The aim of Test 1 is to show the sensitivity of Algorithm 1 with respect to the parameters , Θ,  and .In addition, for the purpose of showing the efficiency of Algorithm 1, in Test 2, four numerical tests are performed using a 3D synthetic image pair (see Test 2 for details) and three 3D medical image pairs, where the 3D medical images in Test 2 can be downloaded from [49][50][51].All the numerical tests are performed under Windows 7 and MATLAB R2012b with Intel core i7-6700 CPU @3.40 GHz and 8 GB memory.In these four tests, the proposed algorithm is compared with Diffeomorphic Log Demons Image Registration (DLDIR for short) algorithm [19,52], MIRT3D algorithm [53], 3D conformal model (One can refer to Appendix B or [54] for details) and MATLAB image registration function "imregister", where the open MATLAB code of DLDIR algorithm and MIRT3D algorithm can be downloaded from [52] and [53], respectively.For quantitative comparison, we compare the following four indexes: -Relative sum of squared difference (Re-SSD for short) which is defined by Re-SSD(, , u) = SSD( (x + u), ) where SSD(, ) = 1 2 ∑︀ ,, ( ,, −  ,, ) 2 .-Mesh folding number (MFN for short) which is defined by where det and for any set , ♯() denotes the number of elements in .
By giving different , Θ, , , we use Algorithm 1 to perform the registration on data I. Final registration results are recorded for each given , Θ,  and .One can see Figure 2 for details.
By Figure 2a, the change of  nearly has no obvious effect on final Re-SSD.This is mainly caused by the fact that (u) plays a dominant role in cost functional (u), because Θ is a very large number.For this reason, the role of ∫︀ Ω |∇  u| 2 dx is very limited compared with (u).By Figure 2b, we know that Re-SSD is also decreased when Θ becomes smaller.However, though we can get a much smaller Re-SSD when Θ is smaller, mesh folding occurs.Hence, Θ is greater than 5000 in our numerical tests.By Figure 2c, we know that smaller  leads to smaller Re-SSD.In addition, Re-SSD changes only a little when  ∈ (0, 0.5).Especially, there is no obvious change on Re-SSD as  ∈ (0, 0.01).In fact,  1 ⊃  2 ( 1 <  2 ).This implies that (u) achieves a much smaller minimum on  1 .Therefore, our numerical results coincide with qualitative analysis.At last, we emphasize that  should not be too large (In our test,  = 0.01).As it is shown on Figure 2c, final Re-SSD = 46.2% as  = 1.This is a very bad registration result.Concerning the sensitivity of , one can notice from Figure 2d that the proposed algorithm is not sensitive to the fractional order .This is mainly caused by the fact that Θ is a large number which weakens the role of fractional regularization ∫︀ Ω |∇  u| 2 dx.
Remark 4.2.In Theorem 3.5, a sufficient technical condition  ≥ 3CML (a not small number) is given to ensure that ℱ is a monotone operator.By Figure 2a, one can notice that the proposed Algorithm 1 is not  sensitive to the parameter .Therefore, in practice, one can also set a small  which may not satisfy the technical condition  ≥ 3CML but may help to get a more accurate result.
Test 2. In this test, we perform registration on a synthetic 3D image pair (data II) and three 3D medical pairs: brain image (data III), liver (data IV) image, lung image (data V).Data II is of the size 80 × 80 × 80.The floating image  (•) and target image (•) of data II are defined by Data III is a 3D brain image pair of size 181 × 217 × 165.The original data can be downloaded from [49].Data IV is a 3D liver image pairs of size 512 × 512 × 122.The original data can be downloaded from [50].Data V is a 3D lung image pairs of size 128 × 128 × 64.The original data can be downloaded from [51].To verify the efficiency of Algorithm 1, we make a quantitative comparison among five different image registration algorithms: DLDIR algorithm, MIRT3D algorithm, MATLAB, 3D conformal model [54] and proposed Algorithm 1.For this purpose, these five algorithms are used to perform the registration on data II-V, respectively.The registration results are listed on Table 1.Note that here the key parameters in the Algorithm 1 are set as  = 2.85, Θ = 7500,  = 5 × 10 −6 and the parameters used in DLDIR, MIRT3D and MATLAB are default values of the code.In addition, to show the 3D data clearly, for each 3D data pair and its registration result, the middlemost

Figure 1 .
Figure 1.Role of  for producing the diffeomorphic deformation.

Remark 4 . 1 .Test 1 .
Theoretically, MFN(u) = 0 implies no mesh folding.In addition, the smaller Re-SSD and the larger SSIM and MI, the better registration result.The aim of this test is to show the sensitivity of Algorithm 1 with respect to the parameters , Θ,  and .For this purpose, numerical tests are performed on a 3D synthetic image pair (data I) with size 80 × 80 × 80.In data I, the floating image  (•) and target image (•) are defined by  (x) = 255 Ω1 (x), (x) = 255 Ω2 (x),

Table 1 .
Quantitative comparison between registration results of five different 3D registration algorithms.