Deformation Estimation with Automatic Sliding Boundary Computation

Abstract

We present a novel method for image registration via a piecewise diffeomorphic deformation which accommodates sliding motion, such as that encountered at organ boundaries. Our method jointly computes the deformation as well as a coherent sliding boundary, represented by a segmentation of the domain into regions of smooth motion. Discontinuities are allowed only at the boundaries of these regions, while invertibility of the total deformation is enforced by disallowing separation or overlap between regions. Optimization alternates between discrete segmentation estimation and continuous deformation estimation. We demonstrate our method on chest 4DCT data showing sliding motion of the lungs against the thoracic cage during breathing.

1 Introduction

Dense image registration driven by a pointwise dissimilarity term is an ill-posed problem, requiring regularization to produce meaningful results. For medical images, smoothness and invertibility of the resulting deformations have become the standard requirements used to define regularization penalties. This represents assumptions regarding the physical system being modeled – that objects should not appear or disappear, and that nearby points should remain nearby. However, certain anatomical motions violate these assumptions. In particular, many organs are not directly attached to adjacent anatomical structures, but are able to slide freely along a boundary. Most notable in CT imaging is the sliding of the lower lungs against the thoracic cage.

We propose a method for jointly computing a segmentation representing regions of smooth motion, and a set of deformations modeling the motion in each region. Driven by a novel formulation of the sliding constraint as a discrete pairwise penalty, optimization alternates between a graph-cut based estimation of the motion segmentation and a continuous optimization of the constituent deformations. The resulting composite deformation is globally invertible and piecewise-smooth. The motion segmentation guarantees that discontinuities due to sliding occur only along coherent boundaries in the image, and the joint estimation of mask and deformations allows anatomical structures to be automatically grouped based on the similarity of their motion.

Previous work has attempted to address the errors caused by using globally smooth deformations to represent lung motion. Wu et al. [11] segment regions and separately register them, while adding a penalty to avoid ‘gaps’ (or overlaps) in the resulting deformation. Recent work in this area has either required a precomputed segmentation of the sliding boundary (e.g. [2, 6]) or used a discontinuity preserving regularization on the deformation (such as Total Variation [10] or Bilateral Filtering [5]), which may introduce spurious discontinuities and does not guarantee the invertibility of the deformation. Vandemeulebroucke et al. [9] give an automated method for creating a ‘motion mask’ of the anatomy defining the sliding boundaries, but base this on anatomical features, and not the observed motion. Schmidt-Richberg et al. [7] propose a diffusion-based regularization which allows sliding along organ boundaries. Boundaries are either precomputed or estimated based on an ad-hoc measure of deformation discontinuity at precomputed image edge locations. Work from computer vision has also considered motion-based segmentation (e.g. [8]), but the modalities dictate that deformations model occlusion instead of invertibility/sliding. The most similar to our work is [2], which also formulates the sliding conditions in terms of a region segmentation, but assumes a predefined segmentation on one of the images.

2 Methods

We first formulate or registration problem and constraints for a piecewise diffeomorphic deformation $\bar{\mathrm {\phi }}$. Given two $D$-dimensional scalar-valued images $\mathrm {f}_{0}:\varOmega _{0}\rightarrow \mathbb {R}$ and $\mathrm {f}_{t}:\varOmega _{t}\rightarrow \mathbb {R}$ (assigned artificial time indices $0$ and $t$) where $\varOmega _{0}, \varOmega _{t} \subset \mathbb {R}^{D}$, we attempt to find a deformation $\bar{\mathrm {\phi }}: \varOmega _{t} \rightarrow \varOmega _{0}$ such that $\mathrm {f}_{0}\circ \bar{\mathrm {\phi }}\approx \mathrm {f}_{t}$, where the dissimilarity is measured by a functional of the form

$$\begin{aligned} \mathrm {E}_{\mathrm {data}}(\mathrm {f}_{0}\circ \bar{\mathrm {\phi }}, \mathrm {f}_{t}):=\int _{\varOmega _{t}}\mathrm {D}(\mathrm {f}_{0}\left( \bar{\mathrm {\phi }}(\varvec{x})\right) , \mathrm {f}_{t}(\varvec{x}))\,\mathrm {d}\varvec{x}, \end{aligned}$$

(1)

where $\mathrm {D}:\mathbb {R}\!\times \!\mathbb {R}\rightarrow \mathbb {R}$ is a pointwise dissimilarity term.

Our method formulates piecewise-diffeomorphic deformation $\bar{\mathrm {\phi }}$ via $K$ diffeomorphic deformations, $\varvec{\mathrm {\phi }}:=\{\mathrm {\phi }^{k}\}_{k=1\ldots K}$, $\mathrm {\phi }^{k}: \varOmega _{t} \rightarrow \varOmega _{0}$ (where we assume identical codomains $\varOmega _{0}$ for simplicity), and a segmentation $\mathrm {\ell }_{0}: \varOmega _{0}\rightarrow \{1\ldots K\}$ defining the region over which each deformation has effect.

We define indicator functions of each label as

$$\begin{aligned} \mathrm {\chi }^{k}_{0}(\varvec{y}) := \left\{ \begin{array}{rl} 1 &{}\text { if } \mathrm {\ell }_{0}(\varvec{y}) = k\\ 0 &{}\text { otherwise} \end{array}\right. \end{aligned}$$

(2)

and the indicator function deformed by the corresponding deformation as

$$\begin{aligned} \mathrm {\chi }^{k}_{t} := \mathrm {\chi }^{k}_{0}\circ \mathrm {\phi }^{k}. \end{aligned}$$

We assume that each $\mathrm {\phi }^{k}$ is diffeomorphic, and therefore invertible. In order to maintain invertibility of the composite deformation we require that it does not result in ‘tearing’ or ‘overlap’ between regions. This can be expressed as a pointwise constraint on the sum of indicator functions

$$\begin{aligned} \sum _{k=1}^{K}\mathrm {\chi }^{k}_{t}(\varvec{x}) = 1 \quad \forall \,\varvec{x}\in \varOmega _{t}. \end{aligned}$$

(3)

Intuitively, $\mathrm {\chi }^{k}_{t}=1$ asserts that deformation $\mathrm {\phi }^{k}$ is in effect at point $\varvec{x}$. The constraint stated in terms of the labeling is then at each $\varvec{x}$, $\mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x})\right) =k$ should be true for exactly one $k$. If the statement is not true for any label k, then we have tearing in the deformation at time $t$, and if it is true for two or more labels we have overlap at this point.

We can now define the segmentation at time $t$

$$\begin{aligned} \mathrm {\ell }_{t}(\varvec{x}) := k\ s.t.\,\,\mathrm {\chi }^{k}_{t}(\varvec{x}) = 1 \end{aligned}$$

(4)

which, for $\varvec{x}\in \varOmega _{t}$, is guaranteed to exist and be unique by (3). We can then define a single piecewise-smooth deformation as

$$\begin{aligned} \bar{\mathrm {\phi }}(\varvec{x}) := \mathrm {\phi }^{\mathrm {\ell }_{t}(\varvec{x})}(\varvec{x}). \end{aligned}$$

(5)

It is easy to verify that $\bar{\mathrm {\phi }}(\varvec{x})$ is invertible. It is surjective since each point $\varvec{y}\in \varOmega _{0}$ can be mapped to a point $\varvec{x}\in \varOmega _{t}$ via $\varvec{x}= (\mathrm {\phi }^{\mathrm {\ell }_{0}(\varvec{y})})^{-1}(\varvec{y})$ via the invertibility of the constituent deformations, and injective because $\bar{\mathrm {\phi }}(\varvec{x})$ s.t. $\varvec{x}\in \mathrm {\chi }^{k}_{t}$ is a unique mapping among $\hat{\varvec{x}}\in \mathrm {\chi }^{k}_{t}$ via the invertibility of the constituent deformations, and unique among $\hat{\varvec{x}}\notin \mathrm {\chi }^{k}_{t}$ as $\hat{\varvec{x}}\notin \mathrm {\chi }^{k}_{t}$ implies $\mathrm {\chi }^{k}_{0}(\bar{\mathrm {\phi }}(\hat{\varvec{x}})) \ne k$, contradicting $\varvec{x}\in \mathrm {\chi }^{k}_{t}$.

2.1 Optimization Criteria

We have now given constraints on the labeling $\mathrm {\ell }_{0}$ and deformations $\varvec{\mathrm {\phi }}$ guaranteeing a piecewise-diffeomorphic deformation, and we choose among these valid deformations by optimizing over a composite objective function that balance the tradeoff between data matching and the regularity of the deformation and labeling

$$\begin{aligned} \mathrm {E}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}) := \mathrm {E}_{\mathrm {data}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}) + \lambda _{\mathrm {per}}\mathrm {E}_{\mathrm {per}}(\mathrm {\ell }_{0}) + \lambda _{\mathrm {pdf}}\mathrm {E}_{\mathrm {pdf}}(\mathrm {\ell }_{0}) + \lambda _{\mathrm {reg}}\sum _{k=1}^{K}\mathrm {E}_{\mathrm {reg}}(\mathrm {\phi }^{k}), \end{aligned}$$

(6)

where we reformulate the data term (1) as a function of $\varvec{\mathrm {\phi }}$ and $\mathrm {\ell }_{0}$

$$\begin{aligned} \mathrm {E}_{\mathrm {data}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}) := \sum _{k=1}^{K}\int _{\varOmega _{t}} \mathrm {\chi }^{k}_{t}(\varvec{x})\,\mathrm {D}(\mathrm {f}_{0}\left( \mathrm {\phi }^{k}(\varvec{x})\right) \!\!, \mathrm {f}_{t}(\varvec{x}))\,\mathrm {d}\varvec{x}, \end{aligned}$$

(7)

and the regularization is taken over each constituent deformation (the form of the deformation regularization $\mathrm {E}_{\mathrm {reg}}$ is dependent on the registration method used). The term $\mathrm {E}_{\mathrm {per}}$ is a ‘perimeter’ penalty on the size of the segmentation boundaries, enforcing the notion that discontinuity should be allowed at a limited set of points. In addition, the penalty $\mathrm {E}_{\mathrm {pdf}}$ is a region-based intensity segmentation penalty, encouraging the labeling $\mathrm {\ell }_{0}$ to segment regions of similar intensity. This helps produce ‘anatomically reasonable’ results (where segmentations follow organ boundaries) in regions such as the top of the lungs, where no sliding is observed and therefore motion segmentation is ambiguous.

The problem we wish to solve can then be stated as

$$\begin{aligned} \underset{\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}}{{\text {argmin}}\;} \mathrm {E}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}) \quad s.t.\,\quad \sum _{k=1}^{K}\mathrm {\chi }^{k}_{t}(\varvec{x}) = 1 \quad \forall \,\varvec{x}\in \varOmega _{t}. \end{aligned}$$

(8)

This is clearly a challenging nonconvex optimization problem. In order to achieve an approximate solution, the hard constraint (3) is relaxed to a penalty,

$$\begin{aligned} \mathrm {E}_{\mathrm {sum}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}):=\left\| \textstyle \sum _{k=1}^{K}\mathrm {\chi }^{k}_{t}- \mathbf {1}\right\| _{1}\!\!, \end{aligned}$$

(9)

We write the relaxed version of (8) as

$$\begin{aligned} \underset{\varvec{\mathrm {\phi }},\mathrm {\ell }_{0}}{{\text {argmin}}\;} \mathrm {E}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}) + \lambda _{\mathrm {sum}}\mathrm {E}_{\mathrm {sum}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}), \end{aligned}$$

(10)

and optimization alternates between finding an optimal labeling $\mathrm {\ell }_{0}$ for a fixed set of deformations $\varvec{\mathrm {\phi }}$, and optimizing over $\varvec{\mathrm {\phi }}$ while keeping $\mathrm {\ell }_{0}$ fixed. This alternating optimizations do not increase the objective (10) at each step, guaranteeing convergence.

2.2 Segmentation Estimation

In the segmentation estimation step we hold the deformations $\varvec{\mathrm {\phi }}$ fixed and optimize (10) over $\mathrm {\ell }_{0}$

$$\begin{aligned} \underset{\mathrm {\ell }_{0}}{{\text {argmin}}\;} \mathrm {E}_{\mathrm {data}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}) + \lambda _{\mathrm {per}}\mathrm {E}_{\mathrm {per}}(\mathrm {\ell }_{0}) + \lambda _{\mathrm {pdf}}\mathrm {E}_{\mathrm {pdf}}(\mathrm {\ell }_{0}) + \lambda _{\mathrm {sum}}\mathrm {E}_{\mathrm {sum}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}). \end{aligned}$$

(11)

Note that for a large enough value of $\lambda _{\mathrm {sum}}$, this becomes a hard constraint; since a single-label segmentation $\mathbf {k}$ s.t. $\mathbf {k}(\varvec{x}) = k\ \forall \,\varvec{x}$ is guaranteed to satisfy the constraints and also have zero boundary cost ($\mathrm {E}_{\mathrm {per}}(\mathbf {k}) = 0$), setting $\lambda _{\mathrm {sum}}= \underset{\mathbf {k}}{{\text {min}}}\;\left( \mathrm {E}_{\mathrm {data}}(\varvec{\mathrm {\phi }}, \mathbf {k})\right) + \epsilon$ guarantees that the optimal solution satisfies the constraint.

As the segmentation $\mathrm {\ell }_{0}$ takes discrete values, we propose a discrete optimization formulation for its estimation. We discretize our problem by defining $\mathrm {f}_{0}$ and $\mathrm {\ell }_{0}$ only on a discrete set of points arranged on a regular grid, $\{\varvec{y}_{n}\}_{n=1}^{N}\subset \varOmega _{0}$, and similarly define $\mathrm {f}_{t}$, $\mathrm {\ell }_{0}$, $\mathrm {\phi }^{k}$, etc. on $\{\varvec{x}_{n}\}_{n=1}^{N}\subset \varOmega _{t}$. For simplicity we will assume equal numbers of points indexed by $n$. For interpolation of images linear interpolation is used, and for label images nearest-neighbor interpolation is used. The objective is written as a discrete set of unary and pairwise functions of node labelings which can be solved exactly for binary labelings (K=2), and approximately otherwise, via graph-cut algorithms [3].

The data matching term $\mathrm {E}_{\mathrm {data}}$ (7) can be written discretely as a unary pointwise penalty on the labeling $\mathrm {\ell }_{0}$, as can $\mathrm {E}_{\mathrm {pdf}}$. We write $\mathrm {E}_{\mathrm {pdf}}$ as a negative log likelihood penalty on $\mathrm {\ell }_{0}$

$$\begin{aligned} \mathrm {E}_{\mathrm {pdf}}(\mathrm {\ell }_{0}) := -\log (\mathrm {p}(\mathrm {f}_{t}|\mathrm {\ell }_{0})) = -\sum _{n=1}^{N}\mathrm {p}(\mathrm {f}_{t}(\varvec{y}_{n})|\mathrm {\ell }_{0}(\varvec{y}_{n})), \end{aligned}$$

(12)

where $\mathrm {p}(\mathrm {f}_{t}(\varvec{y}_{n})|\mathrm {\ell }_{0}(\varvec{y}_{n})=k)$ is a kernel density estimate of the distribution of $\{\mathrm {f}_{t}(\varvec{y}_{n})\ s.t.\,\,\left( \mathrm {\ell }_{0}\right) ^{-}(\varvec{y}_{n})=k\}$, and $\left( \mathrm {\ell }_{0}\right) ^{-}$ is the segmentation from the previous iteration. As this penalty is intended only to resolve ambiguous situations, and not lead to an intensity based segmentation, its influence is kept relatively weak.

We now formulate $\mathrm {E}_{\mathrm {per}}$ and $\mathrm {E}_{\mathrm {sum}}$ discretely as sums of pairwise terms. For ease of notation when working with label values, we define the functions $\mathrm {q}: \{1\ldots K\}\rightarrow \{0,1\}$ (and similarly $\lnot \mathrm {q}$) for indicating equality (or inequality) of labelings ($\mathrm {q}$ for equals):

$$\begin{aligned} \begin{array}{cc} \mathrm {q}(i, j) := \left\{ \begin{array}{rl} 1 &{}\text { if } i = j \\ 0 &{}\text { otherwise} \end{array}\right. \quad &{}\quad \lnot \mathrm {q}(i, j) := \left\{ \begin{array}{rl} 1 &{}\text { if } i \ne j \\ 0 &{}\text { otherwise} \end{array}\right. \end{array} \end{aligned}$$

(13)

Our ‘perimeter’ penalty $\mathrm {E}_{\mathrm {per}}$ is written as a simple Potts model smoothness term on $\mathrm {\ell }_{0}$,

$$\begin{aligned} \mathrm {E}_{\mathrm {per}}(\mathrm {\ell }_{0}) := \sum _{n=1}^{N}\sum _{\varvec{y}_{m}\in \mathcal {N}(\varvec{y}_{n})} \lnot \mathrm {q}\left( \mathrm {\ell }_{0}(\varvec{y}_{n}),\mathrm {\ell }_{0}(\varvec{y}_{m})\right) \!, \end{aligned}$$

(14)

where $\mathcal {N}(\varvec{y}_{n})$ is a set containing the neighbors of $\varvec{y}_{n}$, in our case a 4- or 6-connected neighborhood in 2D or 3D, respectively.

We can directly write the discretized version of $\mathrm {E}_{\mathrm {sum}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0})$ (9) as the sum of $N$ terms (one associated with each $\varvec{x}_{n}$), term $n$ taking as inputs the labelings $\{\mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x}_{n})\right) \}_{k=1}^{K}$

$$\begin{aligned} \mathrm {E}_{\mathrm {sum}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}):= \sum _{n=1}^{N}\left| \left( \sum _{k=1}^{K} \mathrm {q}\left( \mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x}_{n})\right) \!, k\right) \right) - 1\right| \!. \end{aligned}$$

(15)

Due to nearest-neighbor interpolation of the labeling, each $\mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x}_{n})\right)$ is equivalent to $\mathrm {\ell }_{0}(\varvec{y}_{n})$ for some $n\in \{1\ldots N\}$. However, the $K$-input terms associated with each point are difficult and costly to optimize over. Instead, we introduce auxiliary variables representing the labeling $\mathrm {\ell }_{t}$, and optimize over both $\mathrm {\ell }_{0}$ and $\mathrm {\ell }_{t}$. Under the sum constraint (9), $\mathrm {\ell }_{t}$ is uniquely determined via $\mathrm {\ell }_{0}$ and $\varvec{\mathrm {\phi }}$, so we have not introduced any additional unknowns. However it allows us to rewrite each $K$-input term in (15) as the sum of $K$ pairwise terms

$$\begin{aligned} \mathrm {E}_{\mathrm {sum}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}, \mathrm {\ell }_{t}):= \sum _{n=1}^{N}\sum _{k=1}^{K} \mathrm {c}^{k}\left( \mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x}_{n})\right) \!,\mathrm {\ell }_{t}(\varvec{x}_{n})\right) \!, \end{aligned}$$

(16)

where $\mathrm {c}^{k}(i, j) := \lnot \mathrm {q}(\mathrm {q}(i, k), \mathrm {q}(j, k)).$ The function $\mathrm {c}^{k}\left( \mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x}_{n})\right) \!,\mathrm {\ell }_{t}(\varvec{x}_{n})\right)$ indicates whether labels $\mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x}_{n})\right)$ and $\mathrm {\ell }_{t}(\varvec{x}_{n})$ agree that deformation $k$ is (or is not) in effect at $\mathrm {\phi }^{k}(\varvec{x}_{n})$. It takes value zero if $\mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x}_{n})\right) = \mathrm {\ell }_{t}(\varvec{x}_{n}) = k$ (the labelings agree that $\bar{\mathrm {\phi }}(\varvec{x}_{n}) = \mathrm {\phi }^{k}(\varvec{x}_{n})$) or if $\mathrm {\ell }_{0}\left( \mathrm {\phi }^{k}(\varvec{x}_{n})\right) \ne k$ and $\mathrm {\ell }_{t}(\varvec{x}_{n}) \ne k$ (the labelings agree that $\bar{\mathrm {\phi }}(\varvec{x}_{n}) \ne \mathrm {\phi }^{k}(\varvec{x}_{n})$), but takes value 1 if the labelings are in disagreement.

We now have the objective function (11) written as unary ($\mathrm {E}_{\mathrm {data}}$ and $\mathrm {E}_{\mathrm {pdf}}$) and pairwise ($\mathrm {E}_{\mathrm {per}}$ and $\mathrm {E}_{\mathrm {sum}}$) terms on the label values $\mathrm {\ell }_{0}(\varvec{y}_{n})$ and $\mathrm {\ell }_{t}(\varvec{x}_{n})$, which can be optimized by efficient discrete methods [3].

In practice, the effects of discretization mean that enforcing hard constraints give unsatisfactory results, even if the true deformations are known. Instead, $\lambda _{\mathrm {sum}}$ is chosen to allow some constraint violations. In locations of constraint violation the labeling of $\mathrm {\ell }_{t}$ is ambiguous, but the given optimization chooses among ambiguous labels by minimizing the data energy.

2.3 Deformation Estimation

In the deformations estimation step, we hold the labeling $\mathrm {\ell }_{0}$ fixed and optimize (10) over the deformations $\varvec{\mathrm {\phi }}$

$$\begin{aligned} \underset{\varvec{\mathrm {\phi }}}{{\text {argmin}}\;} \mathrm {E}_{\mathrm {data}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}) + \lambda _{\mathrm {reg}}\sum _{k=1}^{K}\mathrm {E}_{\mathrm {reg}}(\mathrm {\phi }^{k}) + \lambda _{\mathrm {sum}}\mathrm {E}_{\mathrm {sum}}(\varvec{\mathrm {\phi }}, \mathrm {\ell }_{0}), \end{aligned}$$

(17)

We represent each $\mathrm {\phi }^{k}$ via a b-spline control grid, and impose additional regularization via $\mathrm {E}_{\mathrm {reg}}(\mathrm {\phi }^{k}):=\Vert \nabla \mathrm {\phi }^{k}\Vert ^2$. While invertibility is not guaranteed by this regularization, it can be easily verified. In practice the regularization is significant (as allowed by the piecewise formulation) and the deformations are sufficiently well-behaved that we do not observe noninvertibility during optimization. However, the optimization is agnostic towards the deformation model, and in particular a flow-based diffeomorphic model could be substituted without additional changes.

The $L^1$ norm in $\mathrm {E}_{\mathrm {sum}}$ (9) is replaced by a differentiable approximation, and optimization is performed via gradient descent.

2.4 Segmentation Initialization

As we are optimizing a highly nonconvex function, a reasonable initialization is needed. As the goal is to segment deformations along regions of motion discontinuity, we automatically create an initial segmentation which separates regions showing different motion. We first generate a deformation on downsampled versions of the input data enforcing only the anisotropic TV penalty

$$\begin{aligned} \underset{\mathrm {\nu }}{{\text {argmin}}\;} \int _{\varOmega _{t}} \mathrm {D}(\mathrm {f}_{0}(\varvec{x}+\mathrm {\nu }(\varvec{x})), \mathrm {f}_{t}) + \sum _{d=1}^{D} |\nabla _{d}\mathrm {\nu }(\varvec{x})|\,\,\mathrm {d}\varvec{x}, \end{aligned}$$

(18)

where $\nabla _{d}$ is the gradient along dimension $d$.

This penalty results in a discontinuous motion field containing piecewise-constant regions. The elements of $\mathrm {\nu }$ are then grouped into $K$ clusters using the k-means algorithm, and the cluster labels are then mapped back to the corresponding spatial locations creating an initial segmentation.

3 Results

We test our method on the publicly available DIR-lab dataset [1], which consists of ten 4DCT datasets as well as landmark correspondences. Resolution (voxel size) is $\approx 1 \!\times \!1 \!\times \!2.5$ mm. Datasets 6–10 contain a large volume of empty space around the patient, and were cropped to a region of interest around the patient body. We use a binary motion segmentation ($K=2$) to represent the breathing motion. A multiscale optimization was used, both on the image resolution and the b-spline grid resolution (two image resolutions, with a total of three b-spline grid refinements for cases 1–5 and four for cases 6–10). Other parameters for both the initialization and optimization steps were empirically chosen, but fixed across all datasets. In order to better register fine features and overcome intensity differences due to density changes in the lung during breathing, we chose a Normalized Gradient Field (NGF) dissimilarity [4] as the pointwise image dissimilarity measure. Figure 1 shows the computed motion segmentation for each dataset. The segmented region includes the lungs as well as regions of visible motion continuing into the organs of the abdominal cavity, while excluding structures such as bone which show little motion. Table 1 gives the landmark errors on these datasets, showing similar results to methods which rely on a precomputed segmentation, or do not guarantee the invertibility of the deformation (average correspondence error of 1.07 mm across all ten datasets, compared to 1.00 mm for method NGF(b) of [4] and 1.01 mm for method pTV of [10]). Figure 2 compares the effect of the piecewise-smooth deformation to that of a single smooth deformation run with the same parameters.

Fig. 1.

Motion segmentation results for all datasets. Top row shows a coronal slice near the spine with motion mask region highlighted in red. Bottom row shows transparent 3D rendering of this region. Bottom number indicates the dataset. (Color figure online)

Table 1. landmark errors (in millimeters) for the $n=300$ landmarks defined at the extreme breathing phases for each dataset

Fig. 2.

Comparison of our piecewise-smooth deformation ((a) and (c)), and a single globally smooth deformation ((b) and (d)), on the DIR-lab dataset 1. The figures on the left show the jacobian determinant of the deformation on a log-scale colormap. Note the clear sliding boundary (dark/red colors where the jacobian is undefined along the tearing boundary) in (a), and the nonphysical deformation of the spine region in (b). On the right are difference images between exhale and deformed inhale for the two methods. Again notice the motion of the spine with a single deformation. (Color figure online)