Discontinuity Preserving Image Registration through Motion Segmentation: A Primal-Dual Approach

Image registration is a powerful tool in medical image analysis and facilitates the clinical routine in several aspects. There are many well established elastic registration methods, but none of them can so far preserve discontinuities in the displacement field. These discontinuities appear in particular at organ boundaries during the breathing induced organ motion. In this paper, we exploit the fact that motion segmentation could play a guiding role during discontinuity preserving registration. The motion segmentation is embedded in a continuous cut framework guaranteeing convexity for motion segmentation. Furthermore we show that a primal-dual method can be used to estimate a solution to this challenging variational problem. Experimental results are presented for MR images with apparent breathing induced sliding motion of the liver along the abdominal wall.


Introduction
Image registration became an indispensable tool for many medical applications including the Image-Guided Therapy systems. Today's image registration methods [1] are powerful at handling rigid as well as nonrigid motion on monoand multimodal images. With the introduction of imaging technologies capable of capturing 4D organ motion, such as 4D-MRI [2], 4D-CT, or 4D-Ultrasound, imagery of sliding organs became a significant focus of research. Most of the state-of-the-art image registration methods cannot yet properly deal with these data sets, as their regularization causes continuous motion fields over the organ boundaries.
In recent publications, the concept of the directiondependent regularization of the displacement field has been introduced; see, for example, [3]. The drawback of these methods is the necessity of providing a good manual segmentation of the boundaries. Approaches that do not rely on prior manual segmentations have been a topic of research for some decades in classical computer vision such as optical flow. Transfer of these research results into the medical field was however scarce. Already in 1989, Mumford and Shah [4] proposed in their pioneering work a functional for image segmentation that avoids spatial smoothing in certain locations of the image, thus preserving discontinuities. Some years later, Weickert and Schnörr [5] proposed an extension of nonquadratic variational regularization for discontinuities preserving optical flow that uses spatiotemporal regularizers instead of flow-driven spatial smoothness. The resulting convex function guaranteed global convergence that could be solved with standard gradient based optimisation schemes. Vese and Chan [6] then introduced a level set framework based approach to efficiently solve the Mumford and Shah minimization problem for segmentation. Another influential approach based on the Total Variation (TV) norm, known to preserve discontinuities, was proposed by Rudin et al. [7], the ROF model. The beneficial behaviour of the TV norm was also exploited in recent registration and optical flow methods, as, for example, by [8,9].
Motion segmentation has been a topic of research for quite some time. In 1993, Nesi [10] proposed a variational optical flow approach that incorporates a variable for the presence of discontinuities and leads to a piecewise smooth motion estimation. Cremers and Schnörr introduced in [11] a variational approach seamlessly integrating motion segmentation and shape regularization into one single energy functional. In [12], the authors formulated the motion segmentation problem in a Bayesian inference framework, 2 Computational and Mathematical Methods in Medicine whereas the motion discontinuity could be represented either as an explicit spline function or by an implicit level set formulation allowing for an arbitrary number of objects to segment. As motion discontinuities mainly appear at object boundaries, it seems natural to combine motion segmentation and registration. Given a good motion segmentation, a proper discontinuous motion field can be estimated and vice versa. Amiaz and Kiryati [13] tried to leverage this property by combining motion segmentation with the optical flow method of Brox et al. [8] in a level set framework. They showed that motion segmentation leads to better discontinuities in the motion field. However, their approach highly depends on the initialisation for the motion segmentation and displacement field. Furthermore, the level set approach is prone to local minima and a rather good initialisation has to be provided in advance.
In this paper we propose an elastic registration approach that handles discontinuities in the motion field by combining motion segmentation and registration in a variational scheme. For motion segmentation we make use of the socalled "continuous cuts" [14,15] scheme that guarantees a globally optimal solution for a fixed motion field. In contrast to our previous work [16], the proposed energy functional is this time TV-1 regularized for both the motion segmentation and the displacement field and the fidelity term is formulated as a sum of absolute values of the constraints. Minimizing the TV-1 -regularized functionals is, however, inherently difficult due to the nonsmoothness of the TV. This is an active field of research, such as [17][18][19][20] and the references therein which are mainly based on operator splitting methods in convex analysis that split energy functionals into the nonlinear and linear terms. We show how to solve the proposed complex energy functional with the primal-dual method of Chambolle and Pock [17], which is very efficient for a wide class of nonsmooth problems. The proposed method is then applied to register 2D MR image pairs with apparent breathing induced sliding motion of the liver along the abdominal wall. The preliminary version of this paper has been published in [21].

Method
In this section we describe the proposed registration method which integrates the displacement field estimation into the convex segmentation method of Chan et al. [14], in order to find smooth displacement fields whilst preserving the discontinuities.

Registration.
Let Ω ⊂ R 2 be the domain of the pixel positions x = ( 1 , 2 ). We then define by the functions : Ω → R and : Ω → R our reference and template image. The aim of image registration is to find a transformation Φ(x) fl x + w(x) such that the relation ∘ Φ ≈ holds. The function with , V : Ω → R, describes the displacement field and will be the intrinsic function we investigate. For convenience we will use the abbreviations w, , and V for w(x), (x), and V(x).
To solve a nonrigid registration problem with expected discontinuities in the displacement field, we want to follow a variational approach and we therefore first consider the energy functional where ∈ R + is a weighting parameter and and are the fidelity term and the smoothness term, respectively, which are defined as The fidelity term incorporates the constraints for the grey value constancy and the gradient constancy and the corresponding weighting parameters are given by 1 , 2 ∈ R + 0 . The smoothness term results in the 1 -norm and, respectively, the vectorial TV of w.
The energy functional is motivated by the energy functional proposed by Brox et al. [8]. Instead of using the approximation function Ψ for the 1 norm, here the pure 1 norm is used for the smoothness term and for the fidelity term the sum of the absolute grey value difference and the componentwise absolute gradient differences is taken.

Motion
Segmentation. Now we would like to integrate the formulation of the energy functional above into the convex segmentation model of Chan et al. [14]. To differentiate the displacement field w into w + and w − , we therefore choose a binary functioñ: where Σ ⊆ Ω ⊆ R 2 , with Σ fl {x ∈ Ω |̃(x) = 1}. By defining (w) fl (w)+ (w) as a data term, we formulate our energy functional as Computational and Mathematical Methods in Medicine where the last term of the above energy is regularization defined by the TV norm and ] is a constant parameter. The energy functional (6) is very much related to the energy functional proposed by Amiaz and Kiryati [13]. Instead of applying the Heaviside function to a level set function, the binary functioñis used here. Furthermore the approximation function Ψ for the 1 norm, which was used in [13], is omitted here. Instead the pure 1 norm is used for the smoothness terms (w ± ) (4) and the fidelity terms (w ± ) are replaced by a sum of absolute values of the constraints (3).
In Figure 1 we illustrate our method by an example. The motion segmentation functioñsplits the displacement field w into the two displacement fields w + and w − .

Iterative Scheme.
To facilitate the optimisation procedure we replace the fidelity term in (3) by its linearised version with w 0 fixed. The minimization of the energy functional (6) with respect to w + , w − , and̃is then performed by the following iterative scheme: (1) For fixed w + and w − , solve (2) For fixed̃, solve (3) For fixed̃, solve Although problem (10) is convex, it is important to note that the overall minimization of the energy functional (6) is a nonconvex optimisation problem. For the minimization of nonconvex energy functionals there exist convex relaxation methods, which are able to provide solutions close to a global minimum. Recently, Strekalovskiy et al. [22] proposed such a method for nonconvex vector-valued labelling problems. In this paper however we will not make use of these kinds of methods. Note that, compared to our previous work [16], this time we did not introduce an auxiliaryṼ in the energy functional. Therefore our iterative scheme consists of only 3 steps instead of 4. The reason for this change is that we intend to use a different numerical approach to solve our problem. More precisely, to solve subproblems (10), (11), and (12) in a fast and efficient way, we follow a primal-dual approach as described by Chambolle and Pock in [17]. We therefore recapitulate in the next section the basic notations and formulations.

Basic Framework for the Primal-Dual
Approach of Chambolle and Pock. First, we define by and two finitedimensional real vector spaces. Their inner products are denoted by ⟨⋅, ⋅⟩ and ⟨⋅, ⋅⟩ , respectively, and their induced norms are given by ‖ ⋅ ‖ = √⟨⋅, ⋅⟩ and ‖ ⋅ ‖ = √⟨⋅, ⋅⟩ , respectively. The general nonlinear primal problem we consider is of the form where : → [0, +∞) and : → [0, +∞) are proper, convex, and lower semicontinuous and the map : → is a continuous linear operator. The corresponding primal-dual formulation of (13) is the saddle-point problem with * : → R ∪ {+∞} being the convex conjugate of . We assume that the problems above have at least one solution (̂,̂) ∈ × and therefore the following holds: where * (̂) and (̂) are the subdifferentials of the convex functions * at̂and at̂. Furthermore we assume that and are "simple"; that is, the resolvent operators ( + * ) −1 and ( + ) −1 are easy to compute. For a convex function the resolvent of the operator at̃can be calculated in our case by In this paper we will only make use of Algorithm 1 in [17] with the extrapolation parameter = 1. Although interesting, the usage of the other proposed algorithms is left for the moment for later research.
To apply Algorithm 1 in [17] to the minimization problems (10), (11), and (12), we first need to rewrite them in their discretized version, then identify the functions and , and finally derive the resolvent operators ( + * ) −1 and ( + ) −1 . For the discrete setting we therefore define by the pixel positions in the image domain with ℎ being the spatial step size. For the calculations of the finite differences, the discrete divergence operator, the discretized inner products, and further details, we refer the reader to [17] and the references therein.
In the following sections we will derive the resolvent operators for the three given minimization problems (10), (11), and (12). (10). Let us consider the continuous problem (10). As mentioned already before, we need to rewrite the problem in its discretized version to be able to apply a primal-dual approach for the minimization. We use the spaces = R and = × and their inner products:

Resolvent Operators for Problem
, , , Using the rectangle rule we can rewrite problem (10) in the discretized version where the factor ℎ 2 can be neglected, since it has no influence on the optimal solution. The norm ‖∇ ‖ 1 , with ∈ , is defined by The solution of the resolvent operator with respect to * can be derived as Computational and Mathematical Methods in Medicine 5 and the one with respect to as Remark 2. The primal-dual formulation of problem (20) is a special case of the general problem considered in Pock et al. 's work [23] and the resulting primal-dual algorithm is then the same as described there. Namely, the dual variable is projected onto the convex set, a disc with radius ], and with a truncation of the primal variablẽthe projection on the feasible convex set [0, 1] is achieved.

Remark 3.
Observant readers probably noticed that we slightly misused the mathematical notation to make the formulas more pleasing to the eye. More precisely, instead of writing (w ± , ) we should have been using (w ± )| , , since discretization is performed after the function is applied. In the following we will also make use of this sloppy notation for the functions 1 , (1) 2 , and (2) 2 given in (9). (11). Now we consider the continuous problem (11). If we have a closer look at this minimization problem, we see that we only receive information about w + on the domain Σ wherẽwill be set to 1. Although theoretically reasonable, for the numerical calculations the implementation gets facilitated by having a smooth extension of w + to the domain Ω \ Σ. We therefore consider instead the problem

Resolvent Operators for Problem
Comparing (11) to (23) the only difference is that the factor is not applied to the smoothness term anymore. For the primal-dual approach we will use this time the spaces = R × R and = × . The inner product of is the same as in (19) and the one for is defined by , , , , , After the discretization of problem (23) we obtain where this time the norm ‖∇ ‖ 1 , with = ( 1 , 2 ) ∈ , is defined by Comparing (25) to (13), we see that the corresponding functions are = ∇, (∇w + ) = ‖∇w + ‖ 1 , and From the resolvent operator with respect to * we obtain The derivation of the resolvent operator with respect to is not that straightforward and more effort has to be put in to find a suitable solution. It is common that the resolvent operators of functions, which are sums of quadratic and absolute norms, lead to so-called thresholding schemes [24,25]. This will be also the case here. Having a closer look at the definition of (26) and (16), we see that we have to solve For̃, = 0 we can conclude from (28) that w + , = w + 0 , . On the other hand, for̃, = 1 we have to distinguish the cases which turn out to be 27 in total. The chosen numbering for the different cases is shown in Table 1.  (29)). (28) with̃, = 1.

Geometric Interpretation of Problem
Before explaining the derivation of the explicit solutions and the needed reformulation of the conditions for the different cases in (29) in more detail, we want to give a deeper insight of the geometric interpretation of problem (28) with̃, = 1.
To facilitate the notation in the following, we introduce the terms Using the definitions in (9) and the notation in (30), we can Note that with 1 (w + , ) = 0 a line , is defined with its normal vector being parallel to the vector , . Similarly, (1) 2 (w + , ) = 0 defines a line , and (2) 2 (w + , ) = 0 a line , . Now we can argue similar to Zach et al. in [25] for the geometric interpretation. The mathematical structure of their minimization problem for which they derive the thresholding scheme is very similar to the one we have in (28) with̃, = 1. The first term in (28), is the squared distance of w + , to w + 0 , , and the terms | 1 (w + , )|, | (1) 2 (w + , )|, and | (2) 2 (w + , )| define the unsigned distances to the lines , , , , and , , respectively. Considering now all w + , with a fixed distance to w + 0 , , we see that problem (28) with̃, = 1 is minimized for w + , closest to the three lines , , , , and , . See Figure 2 for an illustration. (28) with , = 1. We are ready now to derive for each of the cases an explicit solution by using (28). We show the derivation of the solutions only for four different cases. The explicit solutions for the remaining cases can be calculated in a similar fashion as the presented ones. Table 1; that is,

Case 1 and Similar Ones. Let us consider the first case in
Computational and Mathematical Methods in Medicine 7 Solving problem (28) for this case leads to the equation and by using the notation in (31) the explicit solution can then be written as The derivation of the explicit solutions for cases 2, 4, 5, 7, 8, 10, and 11 in Table 1 is performed similarly. There are however other cases left which need another treatment, as, for example, case number 3.

Case 3 and Similar
Ones. The condition for case 3 in Table 1 is given by Here we have to distinguish the two situations , ̸ = 0 and , = 0.
From 1 (w + , ) = 0 we see that the solution w + , has to lie on the line , and of course the distance of it to the other lines , and , should be minimal. Therefore we can assume that the solution is of the form where , is the normal to , . Geometrically this means that w + 0 , is first orthogonally projected to the line , and afterwards moved along the line such that it minimizes the distance to , and , . The unknown 1 is determined by replacing the term (w + , − w + 0 , ) in the equation 1 (w + , ) = 0 using (37), which leads to Solving problem (28) for the current case and replacing again the term (w + , − w + 0 , ) using (37) lead after some simple calculations to the second parameter (39) (ii) If , = 0, we have the following.
Since , vanishes, there is also no line , that we have to consider and we can directly derive the solution w + , by just solving (28), similar to the first case in Table 1, and we get the solution Finally, we can derive the explicit solution for cases 6,9,12,13,14,16,17,22,23,25, and 26 in Table 1 similar to case 3 above.

Case 15 and Similar Ones.
Let us now consider case number 15, which was not covered so far. The condition is this time given by and we see that we have now two equal signs, namely, for . This time we have to distinguish the situations , ̸ = 0, , = 0 ∧ , ̸ = 0, and , = 0 ∧ , = 0.
From 1 (w + , ) = 0 and (1) 2 (w + , ) = 0 we see that the solution w + , has to lie on the lines , and , and should have a minimal distance to the line , . Thus, we assume that the explicit solution is of the form Similar to the last paragraph, replacing the term (w + , − w + 0 , ) in the equation 1 (w + , ) = 0 by making use of (42) leads to the same 1 as in (38). For the parameter 00 we get The case where T , , ̸ = 0 in (43), which is equivalent to , ∦ , , indicates that the lines , and , are not parallel and that they therefore intersect each other in one point. The solution w + , we look for will be at this intersection and the corresponding 00 is calculated by replacing the term (w + , − w + 0 , ) in the equation (1) 2 (w + , ) = 0 using (42). For the case where T , , = 0 in (43), which is equivalent to , ‖ , , we conclude that the lines , and , are parallel. This means that if the solution w + , lies on the line , it necessarily also lies on the line , and the parameter 00 is calculated by solving problem (28) for the current case and replacing the term (w + , − w + 0 , ) using (42). Since , vanishes, there is no need to consider the line , and we focus instead on the line , on which the solution w + , has to lie too. We therefore assume the solution has the form Replacing this time the term (w + , − w + 0 , ) in the equation (1) 2 (w + , ) = 0 with the help of (44) results in After solving problem (28) and replacing the term (w + , − w + 0 , ) using (44) we get (iii) If , = 0 ∧ , = 0, we have the following.
This time both , and , vanish and the corresponding lines , and , do not play a role in the derivation of the solution w + , anymore. We can directly solve problem (28) without any prior assumption, similarly as done before in other cases, and we get The derivation of the explicit solutions for cases 18, 19, 20, 24, and 27 in Table 1 is done similarly to case 15 explained above.
In this situation we do not have to consider the lines , , , , and , at all and solving (28) results in
As one can see in Figure 3, we indicated a cube in the coordinate system that will facilitate for us the problem of finding proper partitioning of the considered space. The edge length of the cube is chosen in such a way that the critical values of the reformulated conditions are covered. The reader will understand the meaning of this sentence afterwards in a better way. For later need, we number the faces of the cube as shown in Figure 4.
The orientation of the face numbers in Figure 4 indicates from which direction we will look at the faces.
In the following we want to show how we can rewrite the conditions for the different cases. The same four groups of cases which were discussed in Section 3.4.2 are considered here too and for each of them the reformulation of the conditions can be achieved in a similar fashion.

Case 1 and Similar Ones.
We start again with the first case in Table 1, which belongs to the group of cases that can be handled in the easiest way compared to the other ones. The condition for this case is given in (33) and by using (35) to replace the term (w + , − w + 0 , ) in the condition we get or rewritten One of the critical values we mentioned before concerning the edge length of the cube shown in Figure 3 can be determined now for case number 1, namely, by the values on the righthand side in (60). The reformulation of the conditions for cases 2, 4, 5, 7, 8, 10, and 11 in Table 1 is achieved similarly to this case and the remaining seven critical values, which are necessary to define the edge length of the cube, are determined in the same way as above.
Remark 4. Before jumping to the next group of cases, we want to explain why the cube in Figure 3 will help us in 5 2 (w + 0 , ), and (2) 2 (w + 0 , ). Our aim is to divide the space into 27 parts, each one defining a region for which exactly one of the 27 cases may apply with their already derived explicit solutions. This partitioning should not be done just arbitrarily and has to be meaningful and to conform to the already set up formulations. Before, we mentioned once that we want the eight critical values from the paragraph before to be covered by the cube. The reason for this is that like this it is possible to reduce the partitioning problem of the whole space to a partitioning problem of the cube. The uniquely defined regions of the cube are the ones from the corresponding eight cases discussed in the paragraph above and they define the corners of it. In Figure 5 we show the location of the regions for different case numbers. For simplicity we drew all the regions with the same size, although they can and most likely will have different sizes. Note that the region for case number 21 is not visible in the figure, because it fully lies in the interior of the cube.
Let us focus now on one face of the cube. We know the values for the boundaries of the corner regions from the critical values we got from the previous paragraph. They are fixed and we introduce for them the names V , ℎ , V , ℎ , V , ℎ , V , and ℎ as shown in Figure 6.
For face number 6, for example, the values V and ℎ will be the boundary values to region 1 (see Figure 5) and after having a look at (60) we can conclude that V = − T , (− , − , − , ) and ℎ = − T , (− , − , − , ). For each face it always holds that For example, to see that V ≤ V we consider the reformulation of the condition for case number 2: which can be derived, as mentioned before, similarly to case number 1. For face number 6 we see in the same manner as before by having a look at Figure 5 and (62) that V = − T , (+ , − , − , ). Since T , , = 1 | , | 2 ≥ 0 (see (31)), we can conclude that V ≤ V .
Let us consider again an arbitrary face of the cube. Although we have boundaries for the corner regions, the ones for the regions in between are not fully defined. Since a rectangular division of the face seems to be adequate, we investigate the problem on how to achieve such one. In Figure 7 a face is depicted with its corner regions and two ways of possible divisions into 9 rectangles are shown, which we call the "clockwise" and "anticlockwise" division method.
Depending on the values for V , ℎ , V , ℎ , V , ℎ , V , and ℎ sometimes both the clockwise and the anticlockwise division method work, as, for example, in Figure 7. Depending on the situation however, there can be also values for which only the clockwise or the anticlockwise or even none of both methods will work, where in the latter situation another adequate division method has to be defined.
In Table 2 we list the different possible situations for the values V , ℎ , V , ℎ , V , ℎ , V , and ℎ . Furthermore, we split the last situation S9 into its parts for later use and introduce the corresponding labelling in Table 3.
In the following, we go through the different situations given in Tables 2 and 3 and decide which kind of method can Table 3: Splitting of situation S9 in Table 2 into its parts.
be used to achieve a division of the face into 9 rectangular regions: (i) Situation S1 in Table 2 is illustrated in Figure 8 for face number 6 and the corresponding region numbers. For this situation, neither the clockwise nor the anticlockwise division method will work, since the corner regions 1 and 5 overlap. An adapted division method, which is depicted in Figure 8, will be used instead for this situation. Note that in this situation certain corner regions, which were defined through the reformulation of the conditions for the group of cases similar to case 1, have to be adapted by removing from them the region where they overlap.
We define the set of situations for which the clockwise division method can be applied by (iii) Situation S9.7 in Table 3 needs again a special treatment to achieve a meaningful division, since neither the clockwise nor the anticlockwise division method can be applied because of the overlapping of the corner regions 2 and 4, if we consider face number 6, for example. The adapted division method, which 5. 6. 14.

13.
1. Figure 8: The division method used for situation S1 in Table 2. For the illustration we used face number 6 and the corresponding region numbers.

13.
1. Figure 9: The division method used for situation S9.7 in Table 3. For the illustration we used again face number 6 and the corresponding region numbers.
will be used instead, is shown in Figure 9 for face number 6 and the corresponding region numbers. Note again that also in this situation certain corner regions have to be adapted by removing from them the region where they overlap.
(iv) Finally, for the situations S4, S7, and S8 in Table 2 and S9.6, S9.8, and S9.9 in Table 3 we can use this time the anticlockwise division method illustrated in Figure 7.
The set of situations for which the anticlockwise division method can be applied is then given by We are now able to give an idea on how to reformulate the conditions for the remaining group of cases.

13
Case 3 and Similar Ones. The reformulation of the cases in this group with the same approach as for case number 1 will not lead to proper partitioning of the space shown in Figure 3. Instead we will use a more geometrically based approach which depends on the different situations the boundary values can take. First we can derive for sure, independently of which approach we use, that This can be seen by considering Figure 5, from which we can deduce that region 3 has to lie between regions 1 and 2. By using the critical values of these adjacent regions, which can be determined from (60) and (62), we arrive at the inequality above.
For the reformulation of the condition with respect to (1) 2 (w + 0 , ), we consider the different situations in which the boundary values of the corner regions can be for face number 6.
(i) If the boundary values are in situation S1 given in Table 2, the reformulation is given by (ii) If the boundary values are in one of the situations in the set S c defined in (63) or in situation S9.7 in Table 3, we get For the reformulation of the condition with respect to (2) 2 (w + 0 , ), we consider this time the different situations belonging to face number 4.
(i) If the boundary values are in situation S1, the reformulation is given by (ii) If the situation of the boundary values is in the set S c or corresponds to the situation S9.7, we get (iii) Finally, if the situation of the boundary values is in S ac , the reformulation is Finally, the reformulation of the conditions for cases 6, 9, 12, 13, 14, 16, 17, 22, 23, 25, and 26 in Table 1 is achieved similarly to this case.

Case 15 and Similar Ones.
For this group of cases we will again make use of the different situations of the boundary values, but this time only for face number 6.
For the reformulation of the condition with respect to 1 (w + 0 , ), we get the following list: (i) If the boundary values are in situation S1 (ii) If the situation of the boundary values is in the set S c (iii) If the boundary values are in situation S9.7 (iv) If the situation of the boundary values is in the set S ac For the reformulation of the condition with respect to (1) 2 (w + 0 , ), we get the following list: (i) If the boundary values are in situation S1 (iv) If the situation of the boundary values is in the set S ac Finally, we get for the reformulation of the condition with respect to (2) 2 (w + 0 , ) the following list: (i) If the boundary values are in situation S1 (2) The reformulation of the conditions for cases 18, 19, 20, 24, and 27 in Table 1 is done then similarly to case 15 explained above.
Case 21. The last case which is left is case number 21. The corresponding region contains the part of the cube which was not covered so far with the corresponding regions of the other cases. Additionally, since it is unfortunately also possible that there are overlappings for some of the regions which appear at different faces, as, for example, the possible overlapping of region 1 with region 11, these overlappings are also assigned to case number 21 to finally achieve proper partitioning.
One could get the impression that this handling of the overlappings is mathematically rather grubby, but to legitimate this we recall that an optimal solution of problem (28) will result in rather small values for 1 (w + , ), (1) 2 (w + , ), and (2) 2 (w + , ) and therefore should be ideally close to case 21 in Table 1. Extending therefore the region for case number 21 should not be harmful. Furthermore, for small values of 1 (w + , ), (1) 2 (w + , ), and (2) 2 (w + , ) a reformulation of the condition can lead more likely to overlappings of certain region parts. Therefore it seems to be adequate to assign these overlappings to region 21.
It is possible to define region 21 in mathematical terms, but since we explained the content of this region already above this will be rather uninspiring and for the implementation an explicit formulation of the region is also not needed, since it can easily be depicted with the help of the other regions. (12). This section is very similar to Section 3.4 and we mainly have to just replace w + by w − and̃by 1 −̃. To achieve a smooth extension of w − to the domain Σ, we consider now instead of problem (12) the following problem:

Resolvent Operators for Problem
The spaces and and their inner products are defined in the same way as in Section 3.4 and the discretization of (83) is performed in the same manner as before.
The resolvent operator with respect to * is identical to (27). Finally, compared to the section before, there are only very slight changes of the resolvent operator with respect to .

Implementation and Results
Although we have a convex minimization problem with respect to the motion segmentation functioñin (10) and linearised the fidelity term in (8), we should keep in mind that the overall minimization problem remains nonconvex and that we have to update w + 0 and w − 0 regularly. Therefore a coarse-to-fine strategy is applied to avoid the risk of getting stuck in a local minimum during the optimisation. The minimal size of the images at the coarsest level is set to min = 32 and the scaling factor for the pyramid to = 0.9. At the coarsest level we initialise the displacement fields w + and w − trivially with 0.
Because of Remark 1 in Section 2.2 one would expect that the choice of the initialisation for the functioñcan be arbitrary. Indeed, whatever initialisation we choose, a global minimizer for the motion segmentation can be found for fixed w + and w − . But in our case the values for w + and w − are updated regularly during optimisation.
The final iteration scheme consists of two loops. The outer loop iterates over the pyramid levels. In each level an inner loop updates the values for̃, w + , and w − , following steps (1)-(3) in Section 3.1. Thus, in each iteration of this inner loop one step of the primal-dual Algorithm 1 from [17] is performed with the appropriate resolvent operators from the previous sections to updatẽ, w + , and w − . To achieve a better convergence of̃we update w + and w − only every 50th iteration. The motion segmentation Σ is obtained by choosing = 0.5 and setting Σ = Σ( ) fl {x ∈ Ω |̃(x) ≥ } (see Remark 1 in Section 2.2). As soon as the finest level is reached, the inner loop is executed until a certain tolerance or the maximum number of iterations is reached. The final displacement field is then obtained by setting Bicubic interpolation is used to calculate the images (x+w ± ) during the iterations and to obtain the final registered image (x + w).
In Figure 10, we show the registration result of a synthetic example, where a textured circle is moving down diagonally. Although our energy functional is convex with respect to the motion segmentation functioñ, the overall minimization task remains a nonconvex problem. To decrease the influence of the nonconvexity during the minimization procedure, certain terms are linearised and a coarse-to-fine strategy is applied. Experiments for the synthetic example in Figure 10 show that due to this workaround different initialisation forc ould be used to achieve similar results, as shown in Figure 11. In order to show the effect of the smoothness parameter ]   Figure 11: The motion segmentation results̃(b) and the corresponding displacement field results w (c) obtained when using different initialisations for̃(a) for the example and parameters in Figure 10. For the results in Figure 10 itself the reference image was used for initialisation. for the TV norm, we show in Figure 12 that we can obtain smoother results as ] increases. In order to show the performance of the proposed method, we show qualitative and quantitative results. Since we are interested in the discontinuities of the displacement field, we compare the proposed method to methods that are able to preserve discontinuities in the displacement field, in this case, the demon algorithm with anisotropic diffusion filtering [26], the registration algorithm of Brox et al. [8], and our previous work [16].
We show the qualitative comparison in Figure 13. As clearly seen, the proposed method can achieve more crisp discontinuities in the displacement field compared to the demon algorithm with anisotropic diffusion filtering and the registration algorithm of Brox et al. Furthermore, the proposed method managed nicely to separate the motion of the abdominal wall and the one of the organs.
In Figure 14 a quantitative evaluation is shown for 22 different liver image pairs from the sequences of [3]. We chose the parameters for all the methods by optimising them with respect to these image pairs. The parameters of our method were set to 1 = 4, 2 = 1, = 0.2, and ] = 0.1. For the demon algorithm with anisotropic diffusion filtering we could use the suggested parameters, for Brox et al. 's method we used = 5, = 80, and = 0.9, and finally for our previous method we used the parameters suggested there, namely, = 0.4, = 0.05, = 0.05, = 4, and = 0.00001.
To quantitatively assess the performance of the four methods, we calculated the mean squared error (MSE) and the normalised mutual information (NMI), with the grey values scaled from 0 to 1. For all our examples the proposed method performed better than the demon algorithm with anisotropic diffusion and the registration algorithm of Brox et al. We used the Kolmogorov-Smirnov test to check for normality of the results using the R Software package (Version 2.10.1). We considered a significance level of 5% as significant. Thetest showed that the proposed method delivered significantly better results than the demon algorithm with anisotropic diffusion filtering and the method of Brox et al. with both < 0.05. There was however no significant difference in the performance of the proposed method compared to our previous method [16] for the MSE and NMI.
In a next step we compared the running times of the proposed method to the ones of our previous work [16] for the 22 liver image pairs. Both methods were implemented in MATLAB and the experiments were performed on a 64-bit Linux system with 1.2 GHz. We used again the parameters 1 = 4, 2 = 1, = 0.2, and ] = 0.1 for the proposed method and for our previous work the parameters = 0.4, = 0.05, = 0.05, = 4, and = 0.00001. In [16] the displacement fields w + and w − were updated every 10th iteration. To provide a fair comparison, we therefore used the same update frequency for the proposed new method. Furthermore, the maximum number of iterations in the finest level is set to 30000 for both methods. The timing results for the 22 liver image pairs are shown in Table 4 together with the mean and standard deviation. The running times of both methods are comparable and the proposed new method performed with around 100 s slightly faster than the old method. The old method makes use of some optimised in-built MATLAB routines, whereas for the proposed method there is still a high potential to optimise the code.

Conclusion
In this paper we presented a primal-dual method for discontinuity preserving nonrigid registration that makes use of the continuous cuts framework. The so-gained motion segmentation influences the motion estimation positively by sharpening the discontinuities in the displacement field. The minimization of the energy functional was implemented in a coarse-to-fine strategy and exploits the rapidity of the primal-dual algorithm studied in [17]. The experimental results demonstrated desirable performance of the proposed method in comparison with those of the demon algorithm with anisotropic diffusion filtering [26] and the registration algorithm of Brox et al. [8].  Ex. number