Contrast enhanced computerized tomography measurement of vascular blood flow

In this work, we study the measurement of blood velocity with contrast enhanced computed tomography. The transport equation is used as a constraint to obtain stable solutions. The inverse problem is formulated as an optimal control problem. The density of the contrast agent is reconstructed together with the flow field. The inversion scheme is tested on a simple phantom. The reconstruction of the velocity is improved but the convergence of the method is slow.


Introduction
The blood velocity in arteries is an important clinical parameter. It can be evaluated with contrast agent enhanced tomography. Yet, in contrast enhanced computed tomography imaging this velocity measurement is not a well established method. Some blood velocity measurement with CT is known as the time-of-flight TOF technique and is based on the difference in arrival time of the contrast agent at two different locations. This measurement of the intravascular blood velocity is based on the temporal changes in the sinogram obtained in sequential scanning mode. The method is limited to the estimation of the mean of the speed of the flow component along the flow field propagation axis [1].
The reconstruction of the flow velocity and of the density of the contrast agent can be formulated as a dynamic inverse problem in which the measuring process and the data are time dependent [2,4,10]. Such problems arise in the field of medical imaging for moving or beating organs. In order to obtain stable solutions, time-dependent regularization techniques are required [5,6,7,8]. Some regularization methods based on optimal control problem can be found in [2].
We present here a method to evaluate the blood velocity from spiral CT scans, consisting of repeated spiral acquisitions with alternating scan directions. The method is based on the integration of physical constraints on the flow for the iterative reconstruction. The scalar transport equation and the divergence free condition are used for the reconstruction of two types of functions: the density of the contrast agent f (x, t) and the flow field w(x, t) for x ∈ R 3 and for x ∈ [0, T ]. The inverse problem is thus adressed in a 3D+t framework. The main contribution of this work is to formulate the blood velocity and density reconstruction problem as an optimal control problem. We consider here a simple geometry where the flow field to recover is uniform to validate the approach. This paper is organized as follows. In Section 1, the dynamic tomography problem considered is presented and it is formulated as an optimal control problem. To solve this problem, we derive first-order optimality conditions. We obtain a system of PDE with initial and boundary conditions. We display numerical results with simple synthetic images that demonstrate the efficiency of the optimal control formulation to improve the reconstruction of the flow field and of the density of the contrast agent.

Optimal control formulation for dynamic tomography
The flow phantom and the CT acquisition geometry used in this study are displayed in Figure 1. This phantom simulates a straight arterial vessel segment oriented along the z axis. We consider here a 4-dimensional (4D) spiral acquisition mode, consisting of repeated spiral acquisitions with projection directions perpendicular to the z axis of rotation.
where p θ(t) are the dynamic measured projections, R θ(t) is the linear Radon operator for the projection angles θ(t) that is available at time t. We assume the contrast agent follows a scalar transport equation with a velocity field, It is injected on a disk D in the plane z = 0. The following transport equation is satisfied: We have added a divergence-free constraint of the flow field w, div(w) = 0 to make the flow smooth and varying not too much inside the moving object. This transport equation and the divergence free condition will be used as a constraints to obtain stable solutions. The cost functional without any constraint can be written: where p δ θ(t) represent the noisy projections and ∇w 2 the Frobenius norm of the Jacobian matrix of the flow field w. The first term is a data term that measures the discrepancy between the measured projections and the projected images. The term α 2 T 0 ∇w 2 2 dt is a smoothing term and α a regularization parameter. We consider two dual variables p : [0, T ] × Ω → R and q : [0, T ] × Ω → R for the constraints, and the Lagrangian L of the constrained problem is given by; where β and γ are regularization parameters.

Optimality system
To minimize the functional, we have to derive the first order optimality condtion with the Euler-Lagrange equations, from which the optimal state f , the optimal speed w, and the adjoint variable p and q can be determined.
2.1.1. Optimality system of p The adjoint equation specifies the first-order necessay condition with respect to the state variable f . The first variation of the Lagrangian L at f ∈ L 2 (Q) in the direction of h ∈ L 2 (Q) is defined by [11]: If a minimizer of L exists, then the equation δL(f, w, p, q)(h) = 0 must hold for every h ∈ L 2 (Q) .
The second term can be written: With integration by parts, (9) with p(T ) = 0. The second term can be written: With the divergence free condition, div(w) = 0, and div(pw) = pdiv(w) + grad(p).w = grad(p).w. We obtain the following equation for p: β( dp dt + grad(p).w) = R t θ(t) (R θ(t) f − g) p(., T ) = 0 This adjoint equation is a backward equation.  3 . We obtain: we obtain, with the variation with respect to the speed: The ground truth solution f * is obtained from the scalar transport equation with a uniform velocity along the z axis, w * = (u * , v * , w * ) = (0, 0, 15). The initial density of the contrast agent is f (x, z = 0, t = 0) = 1 for the pixels inside a disk of diameter 17 around the z axis. A Gaussian white noise with peak-to-peak signal-to-noise ratio (PPSNR) of 48 dB was added to the projections. Our algorithm is not globally convergent algorithms. The scalar diffusion transport equation (Eq.2) was then used to obtain a first guess of the density profile f 0 (x, t). This a priori guess of the solution ensures the convergence of the algorithms.
The following loop was used for the numerical solution the optimal control problem: 1) Solve the forward flow equation to obtain f k (x, t) from the current approximation w k (Eq.2). For the integration of the partial differential equation, the problem is discretized by finite difference and we use an explicit scheme of order 2 with numerical diffusion. For this numerical scheme, the CF L condition is verified, CF L = w dt/dx < 1. For the simulation of the transport equation, the time step dt = 0.005 is used. 2) Solve the backward adjoint equation to obtain the adjoint variable p k (Eq.12). The adjoint equation is solved by the same method as described before by reversing time.
3) Update the flow field w k . Following the work of Tai et al. [9], for the TV-Stokes equation, rather than solving the Eq.17 directly, we use a "gradient descent" strategy to determine the flow w k . Given an initial estimate w 0 k = w k−1 , the gradient flow equation for w k is given by: The search for a minimizer is done using the following iterative procedure to update w k and q k with time step t: The asymptotic state of w k when n → ∞ is chosen as the solution of the Eq.18. t is chosen very small to have a stable scheme. With large t, the iterations will converge faster but if they are two large, the numerical scheme is unstable. For our numerical experiments, the value t = 0.005 leads to convergence [9]. The regularization parameters α, β and γ are chosen to have the best decrease of the discrepancy term. The values α, β, γ used in this work are 2.10 −1 , 2.10 −2 and 2.10 −2 respectively.

Results
In order to validate the method, we have calculated the error term for the velocity field and the reconstructed density with L 2 norms Sections of the true and reconstructed w maps obtained at the end of the iterations for z = 10 and t = 100s are displayed in Figure 2. Sections of the true and reconstructed contrast agent density maps obtained at the end of the iterations for z = 10 and t = 100s are displayed in Figure 3. Some plots of the error decreases are displayed in Figure 4 for the discrepancy terms D(f ) and D(w). A large reduction of the error is obtained for the reconstructed density. The convergence speed for the velocity field is slower.

Conclusion
In this work, we consider the CT measurement of contrast agent flow velocity by using projection data obtained during the spiral acquisition of projections data around the sample as the contrast material flows into the phantom. An optimal formulation of the problem is used in which the scalar diffusion equaition is used. The optimality system for the density of the contrast agent and the flow field are solved iteratively. The efficiency of the method is shown wit a dynamic flow phantom. The reconstruction eror for the contrast agent density is much decreased but the convergence of the velocity field is slow.