Directional Averages for Motion Segmentation in Discontinuity Preserving Image Registration

Abstract

The registration of abdominal images is central in the analysis of motion patterns and physiological investigations of abdominal organs. Challenges which arise in this context are discontinuous changes in correspondence across sliding organ boundaries. Standard regularity criteria like smoothness, are not valid in such regions. In this paper, we introduce a novel regularity criterion which incorporates local motion segmentation in order to preserve discontinuous changes in the spatial mapping. Based on local directional statistics of the transformation parameters it is decided which part of a local neighborhood influences a parameter during registration. Thus, the mutual influence of neighboring parameters which are located on opposing sides of sliding organ boundaries is relaxed. The motion segmentation is performed within the regularizer as well as in the image similarity measure and is thus implicitly updated throughout the optimization. In the experiments on the 4DCT POPI dataset we achieve competitive registration performance compared to state-of-the-art methods.

1 Introduction

Abdominal image registration is an active field of research with many applications such as the analysis of respiratory dynamics or the physiology of abdominal organs as for example the lung. The particular challenges which arise in this scenario are discontinuous changes in the correspondence that occur between organs sliding along each other. Standard smoothness regularity is thus inappropriate at the boundaries between the sliding organs.

Fig. 1.

The dashed circles mark the neighborhood of $c_i$ while the radial shading visualizes the amount of influence a displacement has to the average. In the standard average (left), all displacements in the neighborhood, colored in red, are taken into account while in the directional average (right), only displacements which are approximately aligned to $c_i$ are considered.

In this paper, we present a parametric image registration method comprising a novel regularity criterion. The idea is that the spatial mapping is ought to be locally homogeneous while this criterion should be relaxed at sliding organ boundaries. With local homogeneity we mean that a displacement at a particular image location should be aligned to the average displacement in its neighborhood. To cope with sliding interfaces, we replace the average displacement by a directional average where only displacements which are aligned up to a certain degree are considered (see Fig. 1). The proposed directional average performs a motion segmentation for each parameter separately. The homogeneity requirement for a certain parameter therefore only extends to an aligned neighborhood.

There are major trends which target the sliding organ problem. In [10, 14], image masks of parts in the image sliding along each other are used to register them independently. However, the expected sliding interfaces have to be known in advance. Various approaches [1, 3, 6, 7, 12] adapt their regularity criterion to local image features, as for example image gradients, in order to reduce the influence of the regularity across sliding organ boundaries. The assumption here is that sliding interfaces occur where high intensity changes in the images are present. This applies well to the interfaces between the lung and the thoracic cavity for example. However, it does not hold for generic sliding interfaces. In [5, 9], this problem is addressed by motion segmentation, which is performed exclusively on the spatial transformation. The generic recognized sliding interfaces are thus refined during the optimization. The motion segmentation and the image registration are separately formalized and intertwined by alternating optimization. To achieve a stable convergence remains, however, challenging. In [4, 13, 16], stationary sparse regularity is applied without explicitly considering sliding interfaces. The regularity only amounts to correct the image similarity term at sliding organ boundaries.

The contribution of this paper is twofold. First, we formulate a novel regularity criterion which is based on motion segmentation which is exclusively performed on the transform parameters and thus is implicitly integrated in the registration objective. Hence, it obviates alternated optimization schemes. Second, we modify the image similarity such that the dependency of gradients on opposing sides of sliding interfaces is relaxed. In the experiments on a 4DCT dataset, we compare our method to state-of-the-art registration methods and show competitive registration results.

2 Background

We recap the kernel-framework for image registration on which we base our method. It was elaborated in [3, 4] and we borrow the notation used therein. Let a reference and target image map the d-dimensional input domain to intensity values and a spatial mapping transform the reference coordinate system. Image registration is performed by optimizing

$$\begin{aligned} \mathop {\mathrm{arg\,min}}\limits _u \int _\mathcal {X} \mathcal {L}\left( I_R\left( x+f\left( x\right) \right) ,I_T\left( x\right) \right) dx + \eta \mathcal {R}[f], \end{aligned}$$

(1)

where $\mathcal {L}$ is a loss-function and $\mathcal {R}$ is a regularization term with $\eta$ as trade-off parameter. As transformation model a reproducing kernel Hilbert space (RKHS) is defined

(2)

where is a reproducing kernel and $\Vert \cdot \Vert _\mathcal {H}$ is the RKHS norm. For more details about kernel methods we refer to [2]. In [3], the existence of a finite dimensional solution to Eq. 1 with N pair-wise distinct sampled domain points $x_i$ was shown applying a regularization term operating solely on the finite many transform parameters $c:=\{c_i\}_{i=1}^N$

$$\begin{aligned} \mathcal {R}[f]:=g(p(c)), \end{aligned}$$

(3)

where is a strictly increasing function and is weakly semi-continuous and bounded from below. As k, the combined Wendland kernel of [3] is used. Note, that in a homogeneous region the direction of a transform parameter is similar to an actual displacement if a radial basis function is used as kernel.

3 Method

Within organ regions, the transformation should be locally homogeneous while discontinuities across sliding interfaces are required. Thus, neighboring parameters should point in similar direction to fulfill the local homogeneity requirement. However, at a certain misalignment this similarity requirement should be relaxed. The assumption is that such misalignments appear for neighboring parameters which are located on opposing sides of sliding interfaces.

Fig. 2.

The geometrical interpretation of s and $\omega$ (a). The angle between $c_i$ and neighboring $c_j$ is given by $\omega$, while the sigmoid function s has the value 1 as long as the angle $\omega <\theta$. In (b), s is plotted with increasing slopes $\sigma =1,2,4,6,\dots ,20$.

As a criterion for the misalignment, we propose that the angle

$$\begin{aligned} \omega (c_i,c_j) = \cos ^{-1}\left( \frac{c_i^Tc_j}{\Vert c_i\Vert _\epsilon \Vert c_j\Vert _\epsilon }\right) \end{aligned}$$

(4)

between the average direction of different sliding regions has to exceed a certain threshold $\theta$, where the $\epsilon$-norm¹ $\Vert \cdot \Vert _\epsilon =\sqrt{\Vert \cdot \Vert ^2+\epsilon }$. We specify a sigmoid function which is one if the angle $\omega$ is small and zero if $\omega$ exceeds a threshold $\theta$

$$\begin{aligned} s(\omega ) = 1-\frac{1}{2}\left( 1+\tanh \left( \sigma \left( \omega -\theta \right) \right) \right) , \end{aligned}$$

(5)

where $\sigma$ controls the slope of s (see Fig. 2).

3.1 Directional Average Regularizer

We define the directional average regularizer using Eq. 3 as

$$\begin{aligned} \mathcal {R}_\text {da}(f) := \sum _{i=1}^N \left\| c_i - \mu (c_i)\right\| _\epsilon . \end{aligned}$$

(6)

The directional average $\mu$ for a certain parameter $c_i$ becomes

$$\begin{aligned} \mu (c_i) = \frac{1}{Z(c_i)} \sum \limits _{j=1, j\ne i}^N c_j k(x_i,x_j) s(\omega _{ij}),\quad Z(c_i) = \epsilon +\sum \limits _{j=1, j\ne i}^N k(x_i,x_j) s(\omega _{ij}). \end{aligned}$$

(7)

For convenience, we write $\omega _{ij}$ for $\omega (c_i,c_j)$. Actually, $\mu (c_i)$ is a weighted average over the remaining parameters $c_j$, where the kernel k serves as weighting function. As we use compact kernels, only neighboring parameters $c_j$ are considered. The sigmoid function s additionally decides if a certain $c_j$ contributes to the average dependent on the angle between $c_i$ and $c_j$.

3.2 Directional Average Similarity Metric

With the directional average regularizer of Eq. 6 we specify what kind of parameter configurations should be preserved or penalized respectively. The problem is, that it is only an additive term in the overall registration objective in Eq. 1 and the similarity metric does not take the directional average into account. This is critical because image similarity gradients across sliding interfaces influence a parameter update in the optimization. This can be verified when looking at the derivative of

$$\begin{aligned} \frac{\partial \mathcal {L}}{\partial c_i}=\frac{\partial \mathcal {L}}{\partial I_R(x+f(x))}\nabla I_R(x+f(x))k(x_i,x) \end{aligned}$$

(8)

where to the right side we differentiate $\mathcal {L}$ with respect to its first argument. The gradient for a parameter $c_i$ is influenced by all neighboring points x regardless of their alignment to $c_i$. We unravel now the strict distinction between similarity and regularization term and propose to modify the similarity metric as follows

$$\begin{aligned} \overrightarrow{\mathcal {D}}[I_R,I_T,f] := \int _\mathcal {X} \mathcal {L}\big (I_R(x+f(x)), I_T(x)\big ) \mathcal {A}\big (f(x),c\big )dx. \end{aligned}$$

(9)

We simplify notation with $f_x = f(x)$ and specify the factor $\mathcal {A}$ as

$$\begin{aligned} \mathcal {A}(f_x,c) = 1-\frac{1}{2}\Bigg \Vert \frac{f_x}{\Vert f_x\Vert _\epsilon } - \frac{\mu _A(f_x,c)}{\Vert \mu _A(f_x,c) \Vert _\epsilon } \Bigg \Vert _\epsilon , \end{aligned}$$

(10)

where the average function $\mu _A$ is defined as

$$\begin{aligned} \mu _A(f_x,c) = \frac{\sum _j c_j k(x_j,x)s(\omega (f_x,c_j))}{\epsilon + \sum _j k(x_j,x)s(\omega (f_x,c_j))}. \end{aligned}$$

(11)

The term $\mathcal {A}$ and the regularizer $\mathcal {R}_\text {da}$ mainly differ in two ways. First, $\mathcal {A}$ compares a full displacement $f_x$ with surrounding parameter vectors $c_j$; there is no explicit exclusion of a parameter $c_i$ in the directional average. Second, we avoid to compare apples and pears by considering only normalized versions of $f_x$ and $\mu _A$. Thus, $\mathcal {A}$ evaluates to one if $f_x$ is aligned to neighboring parameter vectors. In this case, the full image similarity is considered. Otherwise, $\mathcal {A}$ becomes zero and the contribution of the image similarity is discarded. Since k smoothes across sliding organ boundaries the direction of an actual displacement $f_x$ does not necessarily correspond to those of its surrounding transform parameters. In these particular regions, $\mathcal {A}$ masks out the contribution of such a point x to the similarity metric.

4 Results

We evaluate our new Directional Averages Motion Segmentation (DAMS) method on a publicly available abdominal 4DCT dataset. We register on three scale levels using an average stochastic gradient descent optimizer [8] where we upsample the transform parameters using nearest neighbor interpolation. A prominent sliding organ boundary in the dataset is located between the lung and the thoracic cavity well expressible by the structure tensor of the image. Therefore, in the first scale level, we apply the anisotropic kernel introduced in [3], where the basis functions are stretched with respect to the structure tensor of the reference image. In the remaining levels, we use the combined Wendland kernel [3]. We apply the robust Cauchy loss $\mathcal {L}(x,x'):=\frac{\beta ^2}{2}\log \left( 1+\left( x-x'\right) ^2/\beta ^2\right)$ where $\beta =1$ controls the scale. In the directional average computation we have set the angle $\theta =\frac{\pi }{3}=60^\circ$. We will make our implementation and configuration files publicly available².

Table 1. Expected TRE and std dev [mm]. Last column: average mean Ø.

Fig. 3.

The transformation magnitude of a coronal slice through case 7 is visualized (all three scale levels) where we have overlayed the reference image to emphasize the spine. A clean outline of the lung at the lower vertebrae is visible.

Fig. 4.

Left: transform parameters $\{c_i^0\}_{i=1}^{N^0}$ after the first scale level as yellow arrows (amplified by 10). Right: sub-sampled $f^3$ after the third scale level as yellow arrows. Background: a coronal slice of the transformed target image.

4.1 POPI Model

We test our method on the 4DCT POPI dataset [15] containing 10 3D-states of a respiratory cycle of a thorax. We empirically set the registration parameters for image number 7 and fixed them for all other time steps. The initial step sizes for the optimizer have been adjusted to each case separately. Image number 1 was the target image. We compare our method to the FFD [11], the pTV [16], the SKM [4] and the bKM [3] method. For the FFD method, we took the target registration errors (TRE) from the POPI homepage³, for the other methods the authors of [16] resp. [3, 4] have kindly provided their TRE values. We calculated the expected TRE [mm] of the 40 first provided ground truth landmarks which are listed in Table 1. Our method performs on par with the tested methods in terms of TRE.

In Fig. 3, we show a sample coronal slice through case 7 where we highlight the transformation magnitude. The transformation becomes finer in higher scale levels and greatly outlines the lung with a clear cutting transition at the lower vertebrae. In Fig. 4, the parameters of the first scale level and final displacements of the third level are plotted. One can clearly identify abrupt directional changes of neighboring parameters at the lower left rib and the lower vertebra.

5 Conclusion

We presented a novel regularity criterion which is targeted to discontinuity preserving image registration. The main contribution is the motion segmentation which is performed exclusively on the spatial transformation and which we integrated into the registration objective. It is based on a directional average of the transform parameters and can be directly optimized using gradient descent. In the experiments with a 4DCT dataset we have achieved competitive registration performance. It would be interesting to investigate further possibilities to adjust the similarity metric based on local motion segmentation.