Deconvolution-based method for image analysis of cerebral blood perfusion computed tomography has been suggested. This analysis is the important part of diagnostics of ischemic stroke. The method is based on total generalized variation regularization algorithm. The algorithm was tested with generated synthetic data and clinical data. Proposed algorithm was compared with singular value decomposition method using Tikhonov regularization and with total variation based deconvolution method. It was shown that the suggested algorithm gives better results than these methods. The proposed algorithm combines both deconvolution and denoising processes, so results are more noisy resistant. It can allow to use lower radiation dose.
Cerebral perfusion computer tomography is an important method for ischemic stroke diagnostics. This test can allow to localize the area of brain damaged by stroke. It is very important for urgent surgery. Usually this information is obtained using method based on the iodinated contrast agent injection and CT scans of its concentration [1]. Example of such CT scan is shown in Fig. 1.
The basic characteristics of cerebral blood dynamics are cerebral blood flow (CBF), cerebral blood volume (CBV) and mean transit time (MTT), calculating at each point of brain.
There are nondeconvolution-based methods of calculation of these characteristics: moment method and method of maximum slope [2, 3].
Better results can be obtained using deconvolution-based methods. Total variation based deconvolution method was suggested in [4, 5]. Nevertheless total variation approach leads to some artifacts, because it tends to make solution piece-wise constant. In this paper we suggest a deconvolution approach using total generalized variation regularization. In this case solution has not piece-wise constant artifacts.
Quality of CT scans depends on the radiation dose: the higher the radiation, the higher the image quality. Thus, the noisy resistant method can allow to use lower radiation dose.
One of the approaches to reducing the noise impact is the denoising of perfusion scans before deconvolution or perfusion maps after deconvolution [6]. The method suggested in our paper includes denoising stage inside the deconvolution procedure. It decreases the number of method parameters and make the results more stable.
Fig. 1. An example of CT scan at one point in time.
We consider some area of brain, where function ctissue(t; x; y) is defined. This function is the concentration of contrast agent obtained from intensity of CT scans at point (x; y). Concentration in artery point is considered as an artery input function (AIF). We denote it by cartery(t).
These functions are connected by the relation given by the convolution equation [2, 3]: ctissue(t; x; y) = (cartery k)(t)
∫ +1 (
; x; y) d ; 1 where k(t; x; y) is a residual function at point (x; y).
Characteristics CBF, CBV and MTT at point (x; y) can be represented via residual function k(t) [3]:
CBF = CBV =
1 tissue 1 max k(t); ∫ 1
k( ) d ; tissue 0
1 ∫ 1 M T T = k( ) d ;
max k(t) 0 where tissue is a constant density of tissue.
We consider all functions on finite uniform grid, so
equation (
B cartery(t2) A = ∆t BB@ ...
Ak = c;
0 cartery(t1) . . .
: : : : : : . . .
cartery(tT ) cartery(tT 1) : : : cartery(t1) There t1; t2; : : : ; tT are the nodes of the uniform grid with step ∆t.
Solution of the deconvolution problem for
equation (
An example of function k(t) obtained without
regularization procedure for system (
There are diferent methods to regularize the problem. In this paper we consider method of total generalized variation and compare it with some state-ofthe-art methods.
Singular value decomposition (SVD) We represent matrix A as a product of three matrices [7]:
A = U VT;
where U and V are orthogonal matrices, and
(
This representation is called singular value decomposition (SVD) of matrix A.
Using SVD, we can write the solution of (
kls = V 1 UT c; (
It is the solution of following minimization
problem:
kls = arg min (jjAk cjj22): (
k2RT
Small singular values make a huge impact on the
values of k. It is the reason of ill-conditioning of
matrix A. These values can be suppressed using
smoothing factor :
(
i i2 + 2 : Using of matrix
1 = diag [ 1(t;ikh); 2(t;ikh); : : : ; r(t;ikh)]
instead of 1 in (
Total variation (TV)
Let consider matrices composed by values of k and c in all considering points of space:
K = [k1; k2; : : : ; kN] ; C = [c1; c2; : : : ; cN] : Generally, regularization method in this paper can be written as a minimization problem of functional:
J (K) = F (K; C) + R(K; ); (
The most common used data fidelity functional is the Frobenius norm of the residual:
F (K; C) = jjAK Cjj22: (
In total variation method we use following regularization functional [4, 5]:
R(K; ) = jjKjjT V = ∑ 1 j Kei+1;j;t Kei;j;tj where Ke 2 RN1 N2 T is the reshaped matrix K, = ( 1; 2) is a regularization parameter.
We can use diferent weights for spatial and temporal derivatives.
Total variation approach may lead to some artifacts, because it makes solution a piece-wise constant. 3.3
Total generalized variation (TGV)
Total generalized variation uses also the second order derivative. In this work we use following form of TGV stabilizer [8]:
T GV 2(z) = 1jj∇zjj1 + 2jj∇(∇z)jj1: (12) Approximation on regular grid for one-dimensional case can be written as
T GV 2(z) = 1 ∑
i + 2 ∑ i jzi+1 jzi+1 zij 2 zi + zi 1j; (13) where z = (z1; z2; : : : ; zT ) is the grid function.
Finally, the regularization functional has the form: R(K; ) = ∑ 1 j Kei+1;j;t
Kei;j;tj
We use Nesterov accelerated gradient descent [9]
for functional minimization (
ϵ ∇F (yk); yk+1 = zk + k(zk+1 zk); where k = 1 3 / (k + 1), and ϵ is the learning rate. In this paper we used ϵ = 10 8.
Clinical perfusion data does not have ground truth values of residue function k(t) and perfusion parameters. Therefore, synthetic data was generated. As a base for generation we take perfusion maps for phantom [10]. Then we generate function k(t) = C e (at)2 . (14) (15)
To evaluate the stability of the method, we add a white additive gaussian noise to the synthetic data.
Residue function k(t) — the result of applying of described methods — is shown in Fig. 4.
We do not have ground truth values for clinical data, but we can compare methods in Fig. 4. We can see that residue function k(t) with TGV does not tends to be a piece-wise constant as it is with TV method.
The obtained mean transit time (MTT) map for synthetic data by diferent method is given in Fig. 5. MTT map is most informative for diagnosis purposes.
0.004 0.003 ic 0.002 t e th 0.001 n yS 0.000 0.0010 0.015 0.010 0.015 0.010 0.005 0.000 SVD
TV TGV Fig. 4. k(t) for a brain tissue point. Ground truth (d) TGV