Computational and Mathematical Methods in Medicine

Computational and Mathematical Methods in Medicine / 2016 / Article

Research Article | Open Access

Volume 2016 |Article ID 9504949 | https://doi.org/10.1155/2016/9504949

Silja Kiriyanthan, Ketut Fundana, Tahir Majeed, Philippe C. Cattin, "Discontinuity Preserving Image Registration through Motion Segmentation: A Primal-Dual Approach", Computational and Mathematical Methods in Medicine, vol. 2016, Article ID 9504949, 20 pages, 2016. https://doi.org/10.1155/2016/9504949

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

Academic Editor: Yuanjie Zheng
Received19 Feb 2016
Revised02 Jun 2016
Accepted14 Jun 2016
Published19 Sep 2016

Abstract

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.

1. 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 mono- and 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 direction-dependent 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, 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 so-called “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- 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--regularized functionals is, however, inherently difficult due to the nonsmoothness of the TV. This is an active field of research, such as [1720] 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].

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

2.1. Registration

Let be the domain of the pixel positions . We then define by the functions and our reference and template image. The aim of image registration is to find a transformation such that the relation holds. The function with , describes the displacement field and will be the intrinsic function we investigate. For convenience we will use the abbreviations , , and for , , and .

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 functionalwhere is a weighting parameter and and are the fidelity term and the smoothness term, respectively, which are defined asThe fidelity term incorporates the constraints for the grey value constancy and the gradient constancy and the corresponding weighting parameters are given by . The smoothness term results in the -norm and, respectively, the vectorial TV of .

The energy functional is motivated by the energy functional proposed by Brox et al. [8]. Instead of using the approximation function for the norm, here the pure 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.

2.2. 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 into and , we therefore choose a binary function where , with . By defining as a data term, we formulate our energy functional aswhere 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 function is used here. Furthermore the approximation function for the norm, which was used in [13], is omitted here. Instead the pure norm is used for the smoothness terms (4) and the fidelity terms are replaced by a sum of absolute values of the constraints (3).

As pointed out by Chan et al. in [14], (6) is strongly related to the Mumford-Shah functional [4] and can be written aswhere denotes the perimeter of the set .

Remark 1. One can show that a global minimizer of can be found by solving the convexified problem and finally setting for almost every with (see Proposition  1 in [16]). In fact, this holds for any data terms and that are measurable.

Finally, we obtain the aimed displacement field by setting .

In Figure 1 we illustrate our method by an example. The motion segmentation function splits the displacement field into the two displacement fields and .

3. Optimisation

3.1. Iterative Scheme

To facilitate the optimisation procedure we replace the fidelity term in (3) by its linearised versionwherewith fixed.

The minimization of the energy functional (6) with respect to , , and is then performed by the following iterative scheme:(1)For fixed and , solve(2)For fixed , solve(3)For fixed , solveAlthough 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 steps instead of . 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.

3.2. Basic Framework for the Primal-Dual Approach of Chambolle and Pock

First, we define by and two finite-dimensional 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 formwhere and 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 problemwith 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 and are easy to compute. For a convex function the resolvent of the operator at can be calculated in our case byIn this paper we will only make use of Algorithm  1 in [17] with the extrapolation parameter . 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 and .

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

3.3. Resolvent Operators for Problem (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 and and their inner products:Using the rectangle rule we can rewrite problem (10) in the discretized versionwhere the factor can be neglected, since it has no influence on the optimal solution. The norm , with , is defined by , where .

Comparing (20) to (13), we see that , , and . The solution of the resolvent operator with respect to can be derived as 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 variable the projection on the feasible convex set 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 we should have been using , since discretization is performed after the function is applied. In the following we will also make use of this sloppy notation for the functions , , and given in (9).

3.4. Resolvent Operators for Problem (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 on the domain where will be set to . Although theoretically reasonable, for the numerical calculations the implementation gets facilitated by having a smooth extension of to the domain . We therefore consider instead the problemComparing (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 and . The inner product of is the same as in (19) and the one for is defined byAfter the discretization of problem (23) we obtainwhere this time the norm , with , is defined by with .

Comparing (25) to (13), we see that the corresponding functions are , , andFrom the resolvent operator with respect to we obtainThe 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 solveFor we can conclude from (28) that . On the other hand, for we have to distinguish the caseswhich turn out to be in total. The chosen numbering for the different cases is shown in Table 1.


Case number

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

3.4.1. Geometric Interpretation of Problem (28) with

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 . To facilitate the notation in the following, we introduce the termsUsing the definitions in (9) and the notation in (30), we can rewrite , , and as Note that with a line is defined with its normal vector being parallel to the vector . Similarly, defines a line and 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 . The first term in (28), , is the squared distance of to , and the terms , , and define the unsigned distances to the lines , , and , respectively. Considering now all with a fixed distance to , we see that problem (28) with is minimized for closest to the three lines , , and . See Figure 2 for an illustration.

3.4.2. Derivation of the Explicit Solutions for Problem (28) with

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.

Case  1 and Similar Ones. Let us consider the first case in Table 1; that is,Solving problem (28) for this case leads to the equation and by using the notation in (31) the explicit solution can then be written asThe derivation of the explicit solutions for cases , , , , , , and in Table 1 is performed similarly. There are however other cases left which need another treatment, as, for example, case number .

Case  3 and Similar Ones. The condition for case in Table 1 is given by Here we have to distinguish the two situations and . (i) If , we have the following.

From we see that the solution 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 formwhere is the normal to . Geometrically this means that is first orthogonally projected to the line and afterwards moved along the line such that it minimizes the distance to and . The unknown is determined by replacing the term in the equation using (37), which leads toSolving problem (28) for the current case and replacing again the term using (37) lead after some simple calculations to the second parameter (ii) If , we have the following.

Since vanishes, there is also no line that we have to consider and we can directly derive the solution 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 , , , , , , , , , , and in Table 1 similar to case above.

Case  15 and Similar Ones. Let us now consider case number , 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 and . This time we have to distinguish the situations , , and . (i) If , we have the following.

From and we see that the solution 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 formSimilar to the last paragraph, replacing the term in the equation by making use of (42) leads to the same as in (38). For the parameter we getThe case where 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 we look for will be at this intersection and the corresponding is calculated by replacing the term in the equation using (42). For the case where in (43), which is equivalent to , we conclude that the lines and are parallel. This means that if the solution lies on the line it necessarily also lies on the line and the parameter is calculated by solving problem (28) for the current case and replacing the term using (42).(ii) If , we have the following.

Since vanishes, there is no need to consider the line and we focus instead on the line on which the solution has to lie too. We therefore assume the solution has the formReplacing this time the term in the equation with the help of (44) results inAfter solving problem (28) and replacing the term using (44) we get (iii) If , 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 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 , , , , and in Table 1 is done similarly to case explained above.

Case  21. There is now one case left, which was not discussed so far. The corresponding condition from Table 1 is This time we have three equal signs and we have to distinguish the situations , , , and . (i) If , we have the following.

Since , , and , we see that the solution has to lie this time on all three lines , , and . For we assume that the explicit solution is of the formLike in the paragraphs before, replacing the term in by making use of (49) leads to being equal to the one in (38). To determine in (49), we argue as follows: (a) If , which is equivalent to , we use (49) in and get (b) If and , that is, and , we use (49) in and get (c) If and , this means that and . Thus all the three lines , , and are parallel and by solving problem (28) and using (49) we get (ii) If , we have the following.

Similarly as in the discussion of (44) in case , we assume the solution is of the formReplacing the term in by making use of (53) leads to being equal to the one in (45). For in (53), we argue as follows: (a) If , which is equivalent to , we make use of (53) in and get (b) If , which means , we get by solving (28) and using (53) (iii) If , we have the following.

Here, and vanish and we focus therefore on the line on which the solution has to lie. We assume the explicit solution is of the formwhere is calculated by using (56) to replace in and is given by (iv) If , we have the following.

In this situation we do not have to consider the lines , , and at all and solving (28) results in

3.4.3. Reformulation of the Conditions for Problem (28) with

Since is unknown but in (9) is known, we have to rewrite the conditions for the cases, which are defined by (29), in terms of , , and . This can be done by incorporating the derived explicit solutions and finding proper partitioning of the space spanned by , , and (see Figure 3). A similar approach was used in the work of Zach et al. [25] to reformulate the conditions for their proposed thresholding scheme.

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 in the condition we get or rewrittenOne 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 , namely, by the values on the right-hand side in (60).

The reformulation of the conditions for cases , , , , , , and 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 finding proper partitioning of the space spanned by , , and . Our aim is to divide the space into parts, each one defining a region for which exactly one of the 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 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 , , , , , , , and as shown in Figure 6.
For face number , for example, the values and will be the boundary values to region (see Figure 5) and after having a look at (60) we can conclude that and .
For each face it always holds that For example, to see that we consider the reformulation of the condition for case number :which can be derived, as mentioned before, similarly to case number . For face number we see in the same manner as before by having a look at Figure 5 and (62) that . Since (see (31)), we can conclude that .
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 rectangles are shown, which we call the “clockwise” and “anticlockwise” division method.
Depending on the values for , , , , , , , 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 , , , , , , , 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 be used to achieve a division of the face into rectangular regions: (i)Situation S1 in Table 2 is illustrated in Figure 8 for face number and the corresponding region numbers. For this situation, neither the clockwise nor the anticlockwise division method will work, since the corner regions and 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 , have to be adapted by removing from them the region where they overlap.(ii)For the situations S2, S3, S5, and S6 in Table 2 and S9.1, S9.2, S9.3, S9.4, and S9.5 in Table 3 we can use the clockwise division method illustrated in Figure 7. We define the set of situations for which the clockwise division method can be applied by