Algorithm for the reconstruction of dynamic objects in CT-scanning using optical flow

TL;DR

This paper introduces a novel method that employs optical flow to estimate and correct for motion in CT scans, improving image quality when objects move during data acquisition.

Contribution

The paper presents a new approach using optical flow to correct motion artifacts in CT imaging, which was not previously demonstrated.

Findings

Successfully estimates motion in simulated CT data

Corrects motion artifacts in reconstructed images

Demonstrates effectiveness of optical flow in CT correction

Abstract

Computed Tomography is a powerful imaging technique that allows non-destructive visualization of the interior of physical objects in different scientific areas. In traditional reconstruction techniques the object of interest is mostly considered to be static, which gives artefacts if the object is moving during the data acquisition. In this paper we present a method that, given only scan results of multiple successive scans, can estimate the motion and correct the CT-images for this motion assuming that the motion field is smooth over the complete domain using optical flow. The proposed method is validated on simulated scan data. The main contribution is that we show we can use the optical flow technique from imaging to correct CT-scan images for motion.

Figures

Tables

Table 1. Table 1 : The RMSE 𝔄 subscript RMSE 𝔄 \text{RMSE}_{\mathfrak{A}} ( 21 ) when estimating the motion using the algorithm 1 with depth d 𝑑 d equal to 0 0 (first and third column) and depth equal to 3 3 3 (second and last column). The noise is normally distributed with mean zero and standard deviation 2 2 2 . We can derive from this table that doing the coarse-to-fine algorithm 1 improves the performance for big motions but that for smaller motion it is better to use depth d = 0 𝑑 0 d=0 . We see that when adding noise to the data, it is preferable to choose the depth d 𝑑 d equal to 3 3 3 .

	depth $d=0$	depth $d=3$	depth $d=0$	depth $d=3$
	no noise added	no noise added	with noise	with noise
Motion 1: Shift	$0.3994$	$0.6664$	$0.7873$	$0.6640$
Motion 2: Rotation	$3.1863$	$1.2257$	$3.4893$	$1.2315$
Motion 3	$2.6076$	$0.4679$	$3.5171$	$0.4658$

Table 2. Table 2 : The error with respect to the exact solution f Δ ​ t / 2 subscript 𝑓 Δ 𝑡 2 f_{\Delta t/2} . First column is the error of the reconstruction if object f 𝑓 f is stationary in the state at time Δ ​ t / 2 Δ 𝑡 2 \Delta t/2 . The second column is the error when we correct for the exact motion. The third column is the error when we correct for the estimated motion. In the last column we put the error when we do not correct for the motion. We see that the error in the first two columns is similar which means our method works. We can see that when using the estimated motion, we still have a far more precise solution than when we do not correct at all.

	${\Vert f_{\Delta t/2}^{\text{ex}}-f_{\Delta t/2}\Vert }_2$	${\Vert f_{\Delta t/2}^{\text{corr}}-f_{\Delta t/2}\Vert }_2$	${\Vert {\hat{f}}_{\Delta t/2}^{\text{corr}}-f_{\Delta t/2}\Vert }_2$	${\Vert f_{\Delta t/2}^{\text{rec}}-f_{\Delta t/2}\Vert }_2$
Motion 1: Shift	$1.3589$	$2.7775$	$4.0105$	$5.4246$
Motion 2: Rotation	$2.6729$	$2.0863$	$3.2873$	$12.8958$
Motion 3	$2.5758$	$2.0819$	$2.9457$	$13.0636$

Table 3. Table 3 : The error with respect to the exact solution f Δ ​ t / 2 subscript 𝑓 Δ 𝑡 2 f_{\Delta t/2} when noise is added to the images. The first column is the error of the reconstruction if object f 𝑓 f is stationary in the state at time Δ ​ t / 2 Δ 𝑡 2 \Delta t/2 . The second column is the error when we correct for the exact motion. The third column is the error when we correct for the estimated motion. The last column shows the error in absence of motion correction. We see that the error in the first two columns is similar which means our method works. We can see that when using the estimated motion, we get an error of similar magnitude

	${\Vert f_{\Delta t/2}^{\text{ex}}-f_{\Delta t/2}\Vert }_2$	${\Vert f_{\Delta t/2}^{\text{corr}}-f_{\Delta t/2}\Vert }_2$	${\Vert {\hat{f}}_{\Delta t/2}^{\text{corr}}-f_{\Delta t/2}\Vert }_2$	${\Vert f_{\Delta t/2}^{\text{rec}}-f_{\Delta t/2}\Vert }_2$
Motion 1: Shift	$3.3296$	$4.4059$	$5.0106$	$6.1845$
Motion 2: Rotation	$3.5557$	$3.4496$	$4.2850$	$13.386$
Motion 3	$3.4983$	$3.4416$	$3.6837$	$13.506$

Equations

$$\mathfrak{R}f(\alpha ,u)={\int }_{L(\alpha ,u)}f(x,y)\mathrm{d}x\mathrm{d}y$$

$$\mathfrak{R}f:[0,\pi [\to \mathbb{R}:(\alpha ,u)\mapsto \mathfrak{R}f(\alpha ,u)$$

$$T:[0,m\pi [\to [0,m\Delta t[:\alpha \mapsto \frac{\alpha }{\pi }\Delta t.$$

$${\mathfrak{R}}_{T(\alpha )}f(\alpha ,u)={\int }_{L(\alpha ,u)}f(x,y;T(\alpha ))\mathrm{d}x\mathrm{d}y.$$

$$t_i:=T(i\pi ),i=0,\dots ,m,$$

$${\mathfrak{R}}_t^{\delta t}f:[\frac{t\pi }{\Delta t}-\pi /2,\frac{t\pi }{\Delta t}+\pi /2[\times \mathbb{R}\to \mathbb{R}:(\alpha ,u)\mapsto {\mathfrak{R}}_{T(\alpha )+\delta t}f(\alpha ,u)$$

$$f^{\text{rec}}=\mathfrak{BR}g(x,y)$$

$$\displaystyle=\int^{\pi}_{0}\mathfrak{F}_{1}^{-1}\Big{(}|\rho|\mathfrak{F}_{1}\big{(}\mathfrak{R}f(\alpha,u)\big{)}\big{(}\rho\big{)}\Big{)}\Big{(}x\cos(\alpha)+y\sin(\alpha)\Big{)}\mathop{\mathrm{\mathstrut d}}\nolimits\!\alpha$$

$$f_t^{\text{rec}}(x,y)$$

$$\frac{\partial f(x,y;t)}{\partial t}+v(x,y)\cdot \nabla f(x,y;t)=0\text{for all}x,y\in \Omega ,$$

$$\forall x,y,t\exists v_x,v_y:f(x,y,t)=f(x+\Delta t{v}_x,y+\Delta t{v}_y,t+\Delta t)$$

$$f(x+\Delta t{v}_x,y+\Delta t{v}_y,t+\Delta t)\approx f(x,y,t)+\frac{\partial f(x,y;t)}{\partial t}+\Delta t{v}_x\frac{\partial f(x,y;t)}{\partial x}+\Delta t{v}_y\frac{\partial f(x,y;t)}{\partial y}$$

$$E(v(x,y))=\int {\int }_{\Omega }{\left(\Delta t\frac{\partial f(x,y;t)}{\partial t}+v(x,y)\cdot \nabla f(x,y;t)\right)}^2+\lambda {\Vert \nabla v(x,y)\Vert }_2^2\mathrm{d}x\mathrm{d}y$$

$$\{\begin{array}{r}\frac{\partial f}{\partial x}\left(\frac{\partial f}{\partial x}{v}_x+\frac{\partial f}{\partial y}{v}_y+\Delta t\frac{\partial f}{\partial t}\right)-\lambda \left(\frac{{\partial }^2{v}_x}{\partial x^2}+\frac{{\partial }^2{v}_x}{\partial y^2}\right)=0\\ \frac{\partial f}{\partial y}\left(\frac{\partial f}{\partial x}{v}_x+\frac{\partial f}{\partial y}{v}_y+\Delta t\frac{\partial f}{\partial t}\right)-\lambda \left(\frac{{\partial }^2{v}_y}{\partial x^2}+\frac{{\partial }^2{v}_y}{\partial y^2}\right)=0.\end{array}$$

$${\mathbf{A}}_{HS}\left[\begin{array}{c}\text{vec}({\mathbf{v}}_x)\\ \text{vec}({\mathbf{v}}_y)\end{array}\right]={\mathbf{b}}_{HS}$$

$${\mathbf{A}}_{HS}=\left[\begin{array}{cc}\text{diag}\left({\left({\mathbf{D}}_x{\mathbf{f}}_t\right)}_{\odot }^2\right)-\lambda \left({\mathbf{D}}_x^2+{\mathbf{D}}_y^2\right)& \text{diag}\left({\mathbf{D}}_x{\mathbf{f}}_t\odot {\mathbf{D}}_y{\mathbf{f}}_t\right)\\ \text{diag}\left({\mathbf{D}}_x{\mathbf{f}}_t\odot {\mathbf{D}}_y{\mathbf{f}}_t\right)& \text{diag}\left({\left({\mathbf{D}}_y{\mathbf{f}}_t\right)}_{\odot }^2\right)-\lambda \left({\mathbf{D}}_x^2+{\mathbf{D}}_y^2\right)\end{array}\right]\in {\mathbb{R}}^{2n^2\times 2n^2}\text{and}$$

$${\mathbf{b}}_{HS}=\left[\begin{array}{c}-\Delta t{\mathbf{D}}_x{\mathbf{f}}_t\odot {\mathbf{D}}_t\mathbf{f}\\ -\Delta t{\mathbf{D}}_y{\mathbf{f}}_t\odot {\mathbf{D}}_t\mathbf{f}\end{array}\right]\in {\mathbb{R}}^{2n^2\times 1}.$$

$$\frac{\partial \mathfrak{R}f}{\partial u}(\alpha ,u)=\cos (\alpha )\mathfrak{R}\frac{\partial f}{\partial x}(\alpha ,u)+\sin (\alpha )\mathfrak{R}\frac{\partial f}{\partial y}(\alpha ,u).$$

$$\left[\begin{array}{c}u\\ w\end{array}\right]=\left[\begin{array}{cc}\cos (\alpha )& \sin (\alpha )\\ -\sin (\alpha )& \cos (\alpha )\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]\text{ and }\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{cc}\cos (\alpha )& -\sin (\alpha )\\ \sin (\alpha )& \cos (\alpha )\end{array}\right]\left[\begin{array}{c}u\\ w\end{array}\right].$$

$$\frac{\partial \mathfrak{R}f}{\partial u}(\alpha ,u)$$

$$\displaystyle=\dfrac{\partial}{\partial u}\int_{-\infty}^{\infty}f\big{(}\cos(\alpha)u-\sin(\alpha)w,\sin(\alpha)u+\cos(\alpha)w\big{)}\mathop{\mathrm{\mathstrut d}}\nolimits\!w$$

$$\displaystyle=\int_{-\infty}^{\infty}\dfrac{\partial}{\partial u}f\big{(}\cos(\alpha)u-\sin(\alpha)w,\sin(\alpha)u+\cos(\alpha)w\big{)}\mathop{\mathrm{\mathstrut d}}\nolimits\!w$$

$$\displaystyle=\int_{-\infty}^{\infty}\cos(\alpha)\dfrac{\partial f}{\partial x}\big{(}\cos(\alpha)u-\sin(\alpha)w,\sin(\alpha)u+\cos(\alpha)w\big{)}+$$

$$\displaystyle\qquad\qquad\sin(\alpha)\dfrac{\partial f}{\partial y}\big{(}\cos(\alpha)u-\sin(\alpha)w,\sin(\alpha)u+\cos(\alpha)w\big{)}\mathop{\mathrm{\mathstrut d}}\nolimits\!w$$

$$=\cos (\alpha ){\int }_{L(\alpha ,u)}\frac{\partial f}{\partial x}\mathrm{d}x\mathrm{d}y+\sin (\alpha ){\int }_{L(\alpha ,u)}\frac{\partial f}{\partial y}\mathrm{d}x\mathrm{d}y.$$

$$\frac{\partial f_t^{\text{rec}}}{\partial x}=\mathfrak{B}{\mathfrak{R}}_t\frac{\partial f}{\partial x}+\mathfrak{B}\sin (\alpha )\left(\cos (\alpha ){\mathfrak{R}}_t\frac{\partial f}{\partial y}-\sin (\alpha ){\mathfrak{R}}_t\frac{\partial f}{\partial x}\right)$$

$$\frac{\partial f_t^{\text{rec}}}{\partial y}=\mathfrak{B}{\mathfrak{R}}_t\frac{\partial f}{\partial y}-\mathfrak{B}\cos (\alpha )\left(\cos (\alpha ){\mathfrak{R}}_t\frac{\partial f}{\partial y}-\sin (\alpha ){\mathfrak{R}}_t\frac{\partial f}{\partial x}\right).$$

$$\frac{\partial f_t^{\text{rec}}}{\partial x}(x,y)$$

$$\displaystyle=\int_{\frac{t\pi}{\Delta t}-\frac{\pi}{2}}^{\frac{t\pi}{\Delta t}+\frac{\pi}{2}}\int_{-\infty}^{\infty}\cos(\alpha)|\rho|2\pi i\rho\mathfrak{F}_{1}\big{(}\mathfrak{R}_{t}f(\alpha,u)\big{)}\big{(}\rho\big{)}\exp\left(2\pi i\rho\big{(}x\cos(\alpha)+y\sin(\alpha)\right)\big{)}\mathop{\mathrm{\mathstrut d}}\nolimits\!\rho\mathop{\mathrm{\mathstrut d}}\nolimits\!\alpha$$

$$\displaystyle=\int_{\frac{t\pi}{\Delta t}-\frac{\pi}{2}}^{\frac{t\pi}{\Delta t}+\frac{\pi}{2}}\int_{-\infty}^{\infty}\cos(\alpha)|\rho|\mathfrak{F}_{1}\Big{(}\dfrac{\partial\mathfrak{R}_{t}f(\alpha,u)}{\partial u}\Big{)}\big{(}\rho\big{)}\exp\left(2\pi i\rho\big{(}x\cos(\alpha)+y\sin(\alpha)\right)\big{)}\mathop{\mathrm{\mathstrut d}}\nolimits\!\rho\mathop{\mathrm{\mathstrut d}}\nolimits\!\alpha$$

$$\displaystyle=\mathfrak{B}\bigg{(}\cos(\alpha)\dfrac{\partial\mathfrak{R}_{t}f(\alpha,u)}{\partial u}\bigg{)}(x,y)$$

$$\displaystyle=\mathfrak{B}\bigg{(}\mathfrak{R}_{t}\frac{\partial f}{\partial x}(\alpha,u)\cos^{2}(\alpha)\bigg{)}(x,y)+\mathfrak{B}\bigg{(}\mathfrak{R}_{t}\frac{\partial f}{\partial y}(\alpha,u)\sin(\alpha)\cos(\alpha)\bigg{)}(x,y)$$

$$\displaystyle=\mathfrak{B}\bigg{(}\mathfrak{R}_{t}\frac{\partial f}{\partial x}(\alpha,u)\bigg{)}(x,y)-\mathfrak{B}\bigg{(}\sin^{2}(\alpha)\mathfrak{R}\frac{\partial f}{\partial x}(\alpha,u)\bigg{)}(x,y)+$$

$$\displaystyle\quad\quad\mathfrak{B}\bigg{(}\sin(\alpha)\cos(\alpha)\mathfrak{R}_{t}\frac{\partial f}{\partial y}(\alpha,u)\bigg{)}(x,y)$$

$$\frac{\partial f_t^{\text{rec}}}{\partial y}$$

$$=\mathfrak{B}{\mathfrak{R}}_t\frac{\partial f}{\partial y}-\mathfrak{B}\cos (\alpha )\left(\cos (\alpha ){\mathfrak{R}}_t\frac{\partial f}{\partial y}-\sin (\alpha ){\mathfrak{R}}_t\frac{\partial f}{\partial x}\right).$$

$$\cos (\alpha )\mathfrak{R}\frac{\partial f}{\partial y}(\alpha ,u)-\sin (\alpha )\mathfrak{R}\frac{\partial f}{\partial x}(\alpha ,u)=0.$$

Full text

††journal: Computational and Applied Mathematics

1 Introduction

CT

Computed Tomography (CT) is a powerful imaging technique that allows non-destructive visualization of the interior of physical objects in medical applications, bio-mechanical research, material science, geology, etc. In current applications, a certain imaging resource and detector, e.g. an X-ray source and detector, are used to acquire two-dimensional projection images of an object, each measured from different directions. From these projections, a three-dimensional virtual reconstruction can then be computed. We refer to Webb (1990) for a review on the origin of computed tomography.

Reconstruction techniques

In practice, the most commonly used analytical method for CT reconstruction is filtered back projection Pan et al. (2009). A major drawback of this method is its inflexibility to different experimental set-ups and its inability to include reconstruction constraints, which can be used to exploit possible prior information to improve the reconstruction of the object. Iterative Algebraic Reconstruction Techniques (ARTs) are a powerful alternative to the aforementioned analytical method by describing the reconstruction problem as a system of linear equations. Algebraic reconstruction methods include SIRT Gregor and Benson (2008), SART Andersen and Kak (1984) and DART Batenburg and Sijbers (2011) and the general class of Krylov solvers such as CG, BiCGStab, GMRES, CGLS and LSQR (of which the latter two are applicable to non-square systems), an overview of which can be found in Simoncini and Szyld (2007).

Time-dependent CT

In both the analytical and algebraic methods, the object of interest is traditionally considered to be static. However, when the object is moving or changing form during the data acquisition process, the flow (direction of movement) also needs to be reconstructed. This requires the calculation of a full 3D reconstruction of both the object and the flow field from the projection measurements over a period of time, yielding a 4D tomographic reconstruction problem, i.e. 3D in space plus 1D for time.

During the previous years, significant progress is made to account for motion during the acquisition process. A first approach is to model the motion in the reconstruction model Mooser et al. (2013); Van Eyndhoven et al. (2012); Li et al. (2005); Van Eyndhoven et al. (2014). A second method sorts the data in subsets such that each subset contains data acquired from a static object. This technique is used, for example when there is periodic motion such as breathing Lu et al. (2006) or heartbeats. Within each subset, where the object does not change, a reconstruction is performed. Finally, if the motion is known upfront, it is possible to compensate for the motion of the object in the CT-scanning Hahn (2014) by using the method of the approximate inverse Louis (1996).

Optical flow

Throughout this paper, we extract the motion using optical flow from the scan data. This is a widely used technique in imaging (see Bardow et al. (2016); Baker et al. (2010) and it exploits the differences between the images to identify patterns of motion. Let $f(x,y,t)$ and $f(x,y,t+\Delta t)$ be two pictures taken with a small time difference $\Delta t$ (we consider 2D images for convenience). Assuming that for every $(x,y)$ holds that $f(x,y,t)=f(x+\Delta x,y+\Delta y,t+\Delta t)$ for some $\Delta x$ and $\Delta y$ (this is the brightness constancy constraint), we can deduce the linear system $\frac{\partial f}{\partial x}v_{x}+\frac{\partial f}{\partial y}v_{y}=-\frac{\partial f}{\partial t}$ using Taylor series. The unknowns $v_{x}$ and $v_{y}$ are the $x$ and $y$ components of the optical flow or the velocity components of the object. There are many ways to use optical flow to estimate the motion from a sequence of images. The most widely used methods are the differential methods that can be classified into local methods like the Lucas-Kanade technique Lucas et al. (1981), global methods like the Horn-Schunck approach Horn and Schunck (1981) and some extensions Basnayake et al. (2013). Local methods calculate the flow per point whereas global methods solve one equation for the entire domain. Recently, methods were developed that are a combination of a local and a global method Bruhn et al. (2005) and also Newton-Krylov methods with regularization have been applied to this problem Mang and Biros (2015).

Outline

The paper is structured as follows. Section 2 introduces the general notations and basic concepts used throughout this article. In section 3 we look at how we can retrieve the motion given just the scan results from a theoretical as well as from a practical point of view. In section 4 we present our method to correct CT-scan images for the motion. Numerical results on both the motion estimation and the correcting of the images are shown in section 5. Finally, we draw some conclusions and we review some research possibilities in section 6. The main contribution is that techniques in imaging to detect motion can also be used to detect the motion in CT-scan images. The numerical results show that it gives also in practice the desired results.

We note discretisations of analytic variables with bold small letters and matrices with bold big letters.

2 Notation and key concepts

We describe everything for 2D CT-scanning, but all definitions can easily be expanded to 3D.

2.1 Modelling scan data

Let $f$ be a time-dependent object, i.e. a function $f:\mathbb{R}^{2}\times[0,T_{f}]\rightarrow\mathbb{R}$ where $[0,T_{f}]$ is the time interval. For notational purposes, let $f_{t}:\mathbb{R}^{2}\rightarrow\mathbb{R}$ be the object $f$ at time $t$. We assume, for all $t$, that $f_{t}$ is an element of $L^{2}(\mathbb{R}^{2})$ and that it is zero outside a square domain $\Omega\subset\mathbb{R}^{2}$ around the center. We first define the used scan model for a stationary object (independent of $t$). In the next paragraph we extend this to time-dependent objects. The scan data per X-ray are modelled by the so-called Radon transformation. The Radon transform for a specific angle $\alpha\in[0,2\pi]$ and shift $u\in\mathbb{R}$ is defined as the integral

[TABLE]

where $L(\alpha,u)=\{(x,y)\in\mathbb{R}^{2}|x\cos(\alpha)+y\sin(\alpha)=u\}$. This integration area is the line perpendicular to the direction $(\cos(\alpha),\sin(\alpha))$ at a shift $u$ of the origin. An illustration of this definition can be found in figure 1. A full scan consists of projection data over angles in an interval of size $\pi$ and shifts $u$ such that the X-ray bunch covers the hole domain $\Omega$. The set of all projection data is called a sinogram and is thus defined as

[TABLE]

for a stationary object. An example is given in figure 2.

If the object moves or changes form we need to generalise the aforementioned definitions. We assume that all projection data for an angle $\alpha$ (so for all shifts $u$) is recorded instantaneously but that each angle has a different recording time $t$. Each complete scan has a total scanning time $\Delta t$ and we assume we scan consecutively $m$ times, so the total scan time is $m\Delta t$. The relation between the angle and the time we acquire data for this angle, is given by

[TABLE]

This means the Radon transform (1) is generalised as

[TABLE]

For notational purposes, we define

[TABLE]

as the start of the $i+1\,$-th scan and/or the end of the $i\,$-th scan. We make the sinograms continually so we define the sinogram at time $t\in[\dfrac{\Delta t}{2},m\Delta t-\dfrac{\Delta t}{2}[$ as

[TABLE]

where the constant $\delta t$ is the time shift. Mostly $\delta t$ is equal to zero, then we just write $\mathfrak{R}_{t}f$ instead of $\mathfrak{R}_{t}^{0}f$.

For the reconstruction of the data, we make use of the filter backprojection theorem see for example Natterer (2001). The backprojection of a sinogram $\mathfrak{R}f$ of a time-independent function $f\in L^{2}(\mathbb{R}^{2})$ is defined by

[TABLE]

where $\mathfrak{F}_{1}$ is the 1-dimensional Fourier transform. Note that the domain of integration for $\alpha$ just needs to be an interval of length $\pi$ and that the interval itself depends on the angles for which we have projection data.

For a time-dependent object $f$, we define the reconstructed object $f^{\text{rec}}_{t}$ at time $t$ as

[TABLE]

This is in fact the filtered backprojection theorem applied on the time-dependent sinogram $\mathfrak{R}_{t}f$. Note that the operator $\mathfrak{R}_{t}$ has a non-trivial null space $\aleph\mathfrak{R}_{t}$ which means that $\mathfrak{B}\mathfrak{R}_{t}(.)$ has a non-empty null space too. This means that in theory, it is possible that a reconstruction is not unique. In practice, the motion is small enough to assume that this does not causes any problem. An illustration of these definitions can be found in figure 3.

2.2 Discretization of Radon transformations

A discrete sinogram for stationary objects (2) can be constructed by dividing the volume of the unknown object in pixels. With each pixel center $(x_{i},y_{j})$, we associate an unknown value $f_{ij}$ and we assume that every pixel is a square. We take in each scan projection data for 180 angles uniformly distributed over $[0,\pi[$. The line integral of a single Radon transform (see (1)) can now be approximated by a weighted sum over these pixels. These weights are the length of the line segment of the projection direction through the pixel (as it needs to approximate an integral) and are entered in a large sparse matrix to form a linear algebra problem $\mathbf{A}\mathbf{x}=\mathbf{b}$. In this matrix $\mathbf{A}$ each row contains the weights of one X-ray and the vectors $\mathbf{x}$ and $\mathbf{b}$ are the unknown pixel values and the discrete measurements respectively. The solution of this linear system is in fact the discretized version of the reconstructed image $f^{\text{rec}}$. In this context, the matrix $\mathbf{A}$ is called the projection operator. There are multiple ways to discretize Radon transformations and to construct this matrix $\mathbf{A}$. We use in this paper the projection scheme described by Joseph in Joseph (1982). In practice we generate this matrix using the Astra Toolbox, see Palenstijn et al. (2011). In section 4 we explain how we can adapt this matrix $\mathbf{A}$ for the motion.

2.3 Evolution

We assume that each pixel moves in a 2D flow field $v(x,y)$ that is yet unknown and that needs to be extracted from the measurements. In this paper we assume that the evolution of the object $f$ itself can be modelled by the optical flow PDE

[TABLE]

where the flow field $v(x,y)$ represents the motion of $f$ on a time interval of size $\Delta t$. The operator $\nabla$ represents the gradient. If we model the flow of $f$ by this equation, we accept that the brightness constancy constraint

[TABLE]

holds for a given time period $\Delta t$. In our case this $\Delta t$ is the scan time of one complete scan. For a specific $x$ and $y$, (5) now follows from the Taylor expansion

[TABLE]

and (6), if we define $v(x,y)=(v_{x},v_{y})$.

This flow field is unknown and needs to be extracted from the measurements. We assume that this flow field is, for short time horizons, independent of the time $t$. This assumption can be made because a CT-scan is applied fast and changes in the flow happens smoothly at a slower time scale (so changes are small in short time periods).

We describe the method to estimate the motion for objects at different time points, in section 3 we describe how we use this method if only scan results are available.

To estimate the flow field $v(x,y)$, we use the Horn-Schunck method Horn and Schunck (1981). This method is a global method (one equation for the entire domain) that prefers flow fields that are smooth.

Here we minimize the following expression with respect to the flow field $v(x,y)$,

[TABLE]

with $\lambda>0$ a regularisation parameter. The first part of the integral is the optical flow expression (see (5)) and the other part is a regularisation term. This means we are searching for a solution $v(x,y)$ such that it satisfies the optical flow equation and that the solution itself is smooth. The extent to which the solution needs to be smooth, is represented by the parameter $\lambda$. Denote with $v_{x}$ and $v_{y}$ the restriction of $v(x,y)$ to the first respectively second variable. If we minimize (7) analytically by the Euler-Lagrange equations, we get

[TABLE]

A proof for (8) can be found in Tarnec et al. (2014). We want to discretize (8) in order to make it useful for doing calculations. We use second order central formulas for all derivatives and we set Neumann boundary conditions on our domain. This means we need three images to make an estimation of the applied motion. We denote by $\mathbf{D}_{i}$ and $\mathbf{D}^{2}_{i}$ the matrix we use to estimate the first respectively the second derivative in the $i$th direction. Further we define $\mathbf{z}_{1}\odot\mathbf{z}_{2}$ as the pointwise multiplication and $\mathbf{z}^{2}_{\odot}$ as the pointwise multiplication with itself. Similar to the discretisation of an object, we associate with the centre the horizontal and vertical component $\mathbf{v}_{x}(x_{i},y_{j})$ and $\mathbf{v}_{y}(x_{i},y_{j})$ of the discretised flow field $\mathbf{v}(x_{i},y_{j})$. If we denote with $\mathbf{D}_{t}\mathbf{f}$ the approximated time derivative with time step $\Delta t$, it is verifiable that the discretisation of (8) leads to a system

[TABLE]

with

[TABLE]

[TABLE]

We denote with $\mathbf{v}$ the exact motion and with $\hat{\mathbf{v}}$ the estimated motion. In this work, we set the parameter $\lambda$ equal to 1. This choice gives us in practice good reconstructions. In future work this can be set automatically.

3 Motion estimation from scan results

3.1 Theoretical derivation

In section 2.3 we have seen how we can estimate the motion in three successive images. In practice, we do not have images but only scan results, so we need to estimate the motion starting from the reconstructed images $f_{t}^{\text{rec}}$ (see (4)). Before we can prove the optical flow equation for the reconstructions (see theorem 5) we first need some lemmas.

Lemma 1**.**

For $f\in L^{2}(\mathbb{R}^{2})$, it holds, for all angles $\alpha$ and $u$, that

[TABLE]

Proof.

We make use of the following transformation and its inverse in the next calculations.

[TABLE]

We obtain that

[TABLE]

∎

We can use this result in the next lemma which gives us a link between the derivatives of the reconstructions $f_{t}^{\text{rec}}$ and the derivatives of the object $f_{t}$.

Lemma 2**.**

For the reconstruction of the object $f$ at time $t\in[\Delta t/2,m\Delta t-\Delta t/2]$, it holds that

[TABLE]

and

[TABLE]

Proof.

For every $(x,y)\in\Omega$ holds that

[TABLE]

where we have used previous lemma in (11). Similarly we can prove that

[TABLE]

∎

In our main theorem 5 we encounter the term $\cos(\alpha)\mathfrak{R}_{t}\dfrac{\partial f}{\partial y}(\alpha,u)-\sin(\alpha)\mathfrak{R}_{t}\dfrac{\partial f}{\partial x}(\alpha,u)$. It is a consequence of following lemma that this term is equal to zero.

Lemma 3**.**

For every angle $\alpha$ it holds, for all $u$, that

[TABLE]

Proof.

We use the same transformation (10) as in lemma 1. For a particular angle $\alpha$ and $u\in\mathbb{R}$, we obtain

[TABLE]

where (12) follows from the fact that $f$ is zero outside a closed area $\Omega$. ∎

The following lemma give us the necessary relation between the object and its reconstruction.

Lemma 4**.**

It holds for every $t$ and $z\in\mathbb{Z}$ that

[TABLE]

Proof.

For every $z\in\mathbb{Z}$ holds that

[TABLE]

∎

We prove an alternative on the optical flow differential equation (9). Instead of $\frac{\partial f}{\partial t}$ we take the second-order Taylor estimation $\dfrac{f_{t+\Delta t}-f_{t-\Delta t}}{2\Delta t}$.

Theorem 5**.**

(Optical flow equation for the reconstructions)

If for all $t\in[\Delta t,(m-1)\Delta t]$

[TABLE]

then it holds for all $t\in[3\frac{\Delta t}{2},(m-\frac{3}{2})\Delta t]$

[TABLE] 2. 2.

If for all $t\in[3\frac{\Delta t}{2},(m-\frac{3}{2})\Delta t]$

[TABLE]

then it holds for $t\in[\Delta t,(m-1)\Delta t]$ that

[TABLE]

for a certain $g\in\aleph\mathfrak{B}\mathfrak{R}_{t}$ (= null space of operator $\mathfrak{B}\mathfrak{R}_{t}$) .

Proof.

If we apply the operator $\mathfrak{B}\mathfrak{R}_{t}$ on (13) we obtain

[TABLE]

where (14), (15) and (16) follows from respectively lemma 4, 2 and 3. 2. 2.

This follows from

[TABLE]

∎

It is possible to achieve a similar result by approximating the time derivative by for example first order approximation formulas as long as the time mesh width is $\Delta t$ (duration of one scan) but tests showed us that the reconstruction had a significant lower quality in that case. The function $g\in\aleph\mathfrak{B}\mathbb{R}_{t}$ in previous theorem is considered to be negligible in practice.

3.2 Implementation

Because the optical flow equation is still valid for the reconstructions (see previous theorem), we can use the method of Horn-Schunck (see section 2.3) to estimate the flow from the reconstructions $f^{\text{rec}}_{t_{i}-\Delta t/2},i=1,\ldots,m$. This means we calculate a reconstruction per scan. In fact, we are not limited to the boundaries of a scan, we can use data which is coming partly from the $i\,$-th scan and partly from the $i+1\,$-th scan. We can reconstruct for example new images $g_{i},i=1,\ldots,m-1$ for which the data for angles between $[\pi/2,\pi[$ is coming from the $i+1\,$-th scan and the data for angles between $[0,\pi/2[$ is coming from the $i\,$-th scan. However, tests showed us that the use of more images does not add significant value and furthermore it is computationally more expensive. From section 2.3, we can retrieve systems

[TABLE]

between three successive reconstructions. Because we assume the motion is constant over time, we know the motion between successive reconstructions is equal. This means we can see the unknown motion $\mathbf{v}$ as the solution of the big system

[TABLE]

By using a single system, we make use of all information available in the scans. It is known in the literature that the Horn-Schunck method can estimate well small motions, but that large displacements can not be retrieved. We can deal with this by applying a coarse-to-fine scheme where we calculate the motion on a lower resolution and use this as an initial solution when calculating the motion on a higher resolution, see for example Anandan (1989). This results in following algorithm where the upper index represents the resolution of the used image. We assume the reconstructed images are objects on a $2^{p}\times 2^{p}$ grid for a certain $p\in\mathbb{N}$. The variable $d$ determines the factor we reduce the resolution in the first estimation of the motion. For the examples in section 5 we take always $d$ equal to $3$. This means we first estimate the motion on a $\frac{2^{p}}{2^{3}}\times\frac{2^{p}}{2^{3}}$ grid.

4 Correcting images for motion

In this section we present a method to correct the reconstruction of a CT-scan for motion when we know the deformation $v(x,y)$ (or an estimate of it) that is performed during one scan. For this we only make use of one scan. The strategy is to move the data recorded at different time steps to a single reference time point and execute the reconstruction there. In fact, we look at how the X-rays are deformed by the motion over time. A consequence is that we also need to have an estimate of the motion on shorter time intervals. We choose to interpolate linearly since no other information is available about the motion on shorter time intervals. So the quality of the reconstruction depends on the extent the linear approximation corresponds with the real deformation. In the next calculations we choose to take the middle of the scan time interval (= $\Delta t/2$) as reference time point, but alternative choices are also possible. It will be clear that the proposed method can easily be adapted such that it reconstructs at another time point. To move the data to the middle of scan, we do following calculations for all angles. For the projection under an angle $\alpha$ performed at time $T(\alpha)$ now holds, for all $u$, that

[TABLE]

So instead of integrating over the path $L(\alpha,u)$, we integrate over the path

[TABLE]

and we obtain that the integral in (18) is equal to

[TABLE]

The objective is now to approximate the integral in (19) by a matrix-vector product $\mathbf{A}_{\text{moved}}\mathbf{f}$ where $\mathbf{f}$ is the discretised and vectorised image at time $\Delta t/2$. Just as the weights on each line in the projection operator $\mathbf{A}$ (see section 2.2) are determined by the length of the segment of the projection line $L(\alpha,u)$ through the pixels, the weights on each row $\mathbf{A}_{\text{moved}}$ are the length of the adapted projection line $L_{\text{moved}}(\alpha,u)$ through the pixels. In practice we determine the values in the matrix $\mathbf{A}_{\text{moved}}$ with following algorithm for a projection angle $\alpha$ and $u\in\mathbb{R}$. Let $\Omega_{\delta}$ be the domain $\Omega$ with an extra border such that we are ensured that all adapted projection lines start and finish outside $\Omega$.

After performing this algorithm we get a new linear system

[TABLE]

with $\mathbf{x}$ the vectorized version of $f_{\Delta t/2}$. We use $f_{\Delta t/2}^{\text{corr}}$ to refer to the solution of this system using the LSQR algorithm.

The accuracy of the algorithm can be further improved by taking more intermediate points between successive intersection points. When doing this, take in mind that the computational time is much higher.

5 Numerical results

In this section, the algorithms on motion estimation and the correcting of CT-scan images are tested and validated on the examples in figure 5. The object in all our simulated examples is the logo of University of Antwerp, see figure 2 (a) on a $256\times 256$ pixel grid. We choose a shift (see figure 5 (a)-(c)), a rotation (see figure 5 (d)-(f)) and a motion which is not an affine transformation (see figure 5 (g)-(i)) as the applied motions. We simulated for every motion field 10 successive scans while the object was moving where per scan we acquire projection data for 180 angles uniformly distributed over $[0,\pi[$.

In the first section we discuss the results for the motion estimation based on algorithm 1 described in section 3.2. To validate the proposed algorithm 2 in section 4 we correct the images for the exact motion and the estimated motion.

5.1 Estimation of motion

In this section, we check the quality of the motion estimate. Because the error of the motion estimate is less relevant in the regions where the information comes only from the regularisation, we calculate the error only for points in the set

[TABLE]

This is the set of points for which the derivative is in absolute value, for at least one moment in the time interval, greater than a certain small default value $\beta>0$. This corresponds with the points for which we have information about the motion. This set is represented in figure 6 for each example. We define the error $\text{RMSE}_{\mathfrak{A}}$ as

[TABLE]

with $n$ the number of elements in the set $\mathfrak{A}$ and $\hat{v}(x,y)$ is the estimated motion from algorithm 1.

In figure 7 the estimated motion is presented together with a plot of the error with the exact motion. In table 1 we calculate the error $\text{RMSE}_{\mathfrak{A}}$ for all three examples for the exact data and for data where we add some normally distributed noise with mean [math] and standard deviation equal to $2$. We choose to calculate it using depth $d=0$ (so without changing the resolution) and depth $d=3$ (we start algorithm 1 by estimating the motion on a resolution which is $8(=2^{3})$ times smaller).

5.2 Correcting CT-images

In this section we validate the algorithm 2 in section 4 by applying it on the examples from figure 5. For this we assume we have access to the exact motion $\mathbf{v}$ or the approximation $\hat{\mathbf{v}}$ made in the previous section 5.1. We assume we have the exact scan data, in the next section we add noise to our data. As mentioned before, we make our reconstructions at the middle of the first scan, i.e. at time $\Delta t/2$. For this we only need the sinogram of the first scan. In table 2 we compare the error with the exact object $f_{\Delta t/2}$ of the following reconstructions: a) the error we make when we perform a reconstruction $f^{\text{ex}}_{\Delta t/2}$ of object $f_{\Delta t/2}$ assuming the object is stationary in this state during the data acquisition, b) reconstructed figure $f_{\Delta t/2}^{\text{corr}}$ when correcting the image for the exact motion $\mathbf{v}$, c) reconstructed figure $\hat{f}_{\Delta t/2}^{\text{corr}}$ when correcting the image for the estimated motion $\hat{\mathbf{v}}$ and d) a reconstruction $f^{\text{rec}}_{\Delta t/2}$ where we do not correct for the motion. It is logical that it is not possible to do better than the reconstruction when there is no motion present so if the reconstruction when correcting for motion has a similar error then we can conclude our method works. We can conclude from this table that the error the algorithm make (columns 2) is similar to the error we make when there is no motion (column 1), so the algorithm to correct for motion works. If we compare the last column with all the other columns we can see that correcting the system for the motion has a significant effect on the quality of the reconstruction. This can also be seen in figure 8.

5.3 Correcting CT-images with noise

We have seen that our algorithm works given exact scan data. The question arises how it performs with added noise on the data. We apply the same noise as in section 5.1. We first check the quality of the reconstruction when we correct for the exact motion using the noisy data and then we test how the method performs when the motion is calculated using the noisy data and the reconstruction is corrected for this estimated motion.

The results of the reconstruction can be found in the next table. We use the same abbreviations as in table 2. We see that adding normal distributed error does not affect much the quality of the reconstruction. The error with the exact motion (see second column in table 3) is also in this case comparable with the error we make when the object is stationary (see first column in table 3). The same conclusions can be drawn from the images of the reconstructions (see figure 9).

We see that adding normal distributed noise has no big impact on the quality of the restrictions. This can be seen in the table as well as in the above figure.

6 Conclusion

We have shown in this paper that the motion can be determined starting from CT-scan data and that we can use this to correct CT-scan images. In fact, we show that techniques from imaging such as optical flow an be used for the dynamic CT-problem. Furthermore, the correction for the motion can be done very efficiently because all x-rays can be calculated in parallel. Although the results look promising, there is still work to refine the proposed method. First of all we have only tested the presented algorithms with simulated data so it would be very interesting to look how it performs with real data. Secondly, we calculate the motion in every pixel although in certain areas we can see that there is no motion present. Certainly when we look at extending our methods to the 3D case, it is not possible anymore to calculate the motion in every single pixel. We already did some experiments how to split up sinogram data into parts where there is motion and into parts where there is no motion present. In the areas where we know that the object stay stationary, we obviously do not need to calculate the motion here. Moreover, we only need to reconstruct this data once. The question arises how to adapt the methods if we only need to perform everything on a small area. Also the automatic determination of the regularisation parameter for the estimation of the motion is something that needs to be done. If this would work, this could significantly improve our methods.