Total Generalized Variation Method for Deconvolution-based CT Brain
Perfusion
D.A. Lyukov1 , A.S. Krylov1 , V.A. Lukshin2
d@lyukov.com|kryl@cs.msu.ru|wlukshin@nsi.ru
1
Laboratory of Mathematical Methods of Image Processing, Faculty of Computational Mathematics and Cybernetics,
Lomonosov Moscow State University, Moscow, Russia;
2
National Medical Research Center of Neurosurgery named after academician N.N. Burdenko, Moscow, Russia
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.
Keywords: computed tomography, cerebral perfusion, deconvolution method, total generalized variation.
1. Introduction
Cerebral perfusion computer tomography is an im-
portant 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 dynam-
ics are cerebral blood flow (CBF), cerebral blood vol- Fig. 1. An example of CT scan at one point in time.
ume (CBV) and mean transit time (MTT), calculating
at each point of brain.
2. Theoretical model
There are nondeconvolution-based methods of cal- We consider some area of brain, where function
culation of these characteristics: moment method and ctissue (t, x, y) is defined. This function is the con-
method of maximum slope [2, 3]. centration of contrast agent obtained from intensity
of CT scans at point (x, y). Concentration in artery
Better results can be obtained using
point is considered as an artery input function (AIF).
deconvolution-based methods. Total variation based
We denote it by cartery (t).
deconvolution method was suggested in [4, 5]. Nev-
These functions are connected by the relation given
ertheless total variation approach leads to some ar-
by the convolution equation [2, 3]:
tifacts, because it tends to make solution piece-wise ctissue (t, x, y) = (cartery ∗ k)(t)
constant. In this paper we suggest a deconvolution ∫ +∞
approach using total generalized variation regulariza- (1)
= cartery (ξ)k(t − ξ, x, y) dξ,
tion. In this case solution has not piece-wise constant −∞
artifacts. where k(t, x, y) is a residual function at point (x, y).
Characteristics CBF, CBV and MTT at point
Quality of CT scans depends on the radiation dose: (x, y) can be represented via residual function k(t)
the higher the radiation, the higher the image qual- [3]:
ity. Thus, the noisy resistant method can allow to use 1
lower radiation dose. CBF = max k(t),
ρtissue
∫ ∞
One of the approaches to reducing the noise im- 1
CBV = k(τ ) dτ, (2)
pact is the denoising of perfusion scans before de- ρtissue 0
∫ ∞
convolution or perfusion maps after deconvolution [6]. 1
MTT = k(τ ) dτ,
The method suggested in our paper includes denoising max k(t) 0
stage inside the deconvolution procedure. It decreases where ρtissue is a constant density of tissue.
the number of method parameters and make the re- We consider all functions on finite uniform grid, so
sults more stable. equation (1) turns to linear system:
Copyright © 2019 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
3.1 Singular value decomposition (SVD)
Ak = c, (3)
We represent matrix A as a product of three ma-
where: trices [7]:
cartery (t1 ) 0 ... 0 A = U Σ VT , (4)
cartery (t2 ) cartery (t1 ) ... 0 where U and V are orthogonal matrices, and
A = ∆t .. .. .. .. .
. . . . Σ = diag [σ1 , σ2 , . . . , σr ]
cartery (tT ) cartery (tT −1 ) ... cartery (t1 ) is a diagonal matrix composed of the singular values
of matrix A, r = rang A.
There t1 , t2 , . . . , tT are the nodes of the uniform grid This representation is called singular value decom-
with step ∆t. position (SVD) of matrix A.
Using SVD, we can write the solution of (3) in the
3. Deconvolution problem form:
Solution of the deconvolution problem for equa- kls = V Σ−1 UT c, (5)
tion (3) is an ill-posed problem, and matrix A is ill- where Σ = diag [σ1−1 , σ2−1 , . . . , σr−1 ].
−1
conditioned. If we solve system (3) directly, solution k It is the solution of following minimization prob-
will be unstable. Small noise in input data can com- lem:
pletely change solution and the solution can not be kls = arg min (||Ak − c||22 ). (6)
used for medical diagnosis. k∈RT
An example of function k(t) obtained without reg- Small singular values make a huge impact on the
ularization procedure for system (3) is shown in Fig. values of k. It is the reason of ill-conditioning of ma-
2–3. The unstability is so large that we can not find trix A. These values can be suppressed using smooth-
any real information on the solution. ing factor λ:
(tikh) σi
σi, λ = 2 . (7)
σi + λ2
10 Using of matrix
Σ−1
(tikh) (tikh) (tikh)
λ = diag [σ1, λ , σ2, λ , . . . , σr, λ ]
5
instead of Σ−1 in (6), we get another method called
Tikhonov regularization. Here λ is a regularization
0 (tikh)
parameter. Vector kλ obtained by this method is
the solution of another minimization problem:
5
(tikh) ( )
kλ = arg min ||Ak − c||22 + λ2 ||k||22 . (8)
k∈RT
10
0 5 10 15 20 25 30 35
3.2 Total variation (TV)
Fig. 2. k(t). Range is limited to interval [−10, 10], so
Let consider matrices composed by values of k and
bigger oscillations are seen as vertical strips.
c in all considering points of space:
K = [k1 , k2 , . . . , kN ] , C = [c1 , c2 , . . . , cN ] .
Generally, regularization method in this paper can
10 9 be written as a minimization problem of functional:
10 8
10 7
J(K) = F (K, C) + R(K, λ), (9)
10 6 where F (K, C) is the data fidelity functional, R(K, λ)
10 5 is regularization term, and λ is a regularization pa-
10 4
rameter.
10 3
10 2 The most common used data fidelity functional is
10 1 the Frobenius norm of the residual:
10 0
10 -1
F (K, C) = ||AK − C||22 . (10)
0 5 10 15 20 25 30 35
In total variation method we use following regular-
ization functional [4, 5]:∑
Fig. 3. |k(t)| on a logarithmic scale. R(K, λ) = ||K||λT V = e i+1,j,t − K
λ1 |K e i,j,t |
i, j, t
∑
There are different methods to regularize the prob- + e i,j+1,t − K
λ1 |K e i,j,t |
(11)
lem. In this paper we consider method of total gen- i, j, t
∑
eralized variation and compare it with some state-of- + e i,j,t+1 − K
λ2 |K e i,j,t |
the-art methods. i, j, t
where K e ∈ RN1 ×N2 ×T is the reshaped matrix K, To evaluate the stability of the method, we add a
λ = (λ1 , λ2 ) is a regularization parameter. white additive gaussian noise to the synthetic data.
We can use different weights for spatial and tem- Residue function k(t) — the result of applying of
poral derivatives. described methods — is shown in Fig. 4.
Total variation approach may lead to some arti-
We do not have ground truth values for clinical
facts, because it makes solution a piece-wise constant.
data, but we can compare methods in Fig. 4. We
can see that residue function k(t) with TGV does not
3.3 Total generalized variation (TGV) tends to be a piece-wise constant as it is with TV
Total generalized variation uses also the second or- method.
der derivative. In this work we use following form of The obtained mean transit time (MTT) map for
TGV stabilizer [8]: synthetic data by different method is given in Fig. 5.
T GVλ2 (z) = λ1 ||∇z||1 + λ2 ||∇(∇z)||1 . (12) MTT map is most informative for diagnosis purposes.
Approximation on regular grid for one-dimensional
case can be written as 6. References
∑
T GVλ2 (z) =λ1 |zi+1 − zi | [1] Axel Leon. Cerebral blood flow determination by
i rapid-sequence computed tomography: theoretical
∑ (13)
+λ2 |zi+1 − 2 zi + zi−1 |, analysis. // Radiology. — 1980. — Vol. 137, no.
i 3. — P. 679–686.
[2] Theoretic basis and technical implementations of
where z = (z1 , z2 , . . . , zT ) is the grid function.
CT perfusion in acute ischemic stroke, part 1: the-
Finally, the regularization functional has the form:
∑ oretic basis / A.A. Konstas, G.V. Goldmakher, T-
R(K, λ) = λ1 |K e i+1,j,t − Ke i,j,t |
Y Lee, M.H. Lev // American Journal of Neurora-
i, j, t
∑ diology. — 2009. — Vol. 30, no. 4. — P. 662–668.
+ e i+1,j,t − 2 K
λ2 |K e i,j,t + K
e i−1,j,t | [3] Deconvolution-based CT and MR brain perfu-
i, j, t sion measurement: theoretical model revisited and
∑ practical implementation details / Andreas Fiesel-
+ e i,j+1,t − K
λ1 |K e i,j,t |
mann, Markus Kowarschik, Arundhuti Ganguly et
i, j, t
∑ al. // Journal of Biomedical Imaging. — 2011. —
+ e i,j+1,t − 2 K
λ2 |K e i,j,t + K
e i,j−1,t | Vol. 2011. — P. 14.
i, j, t [4] Tensor total-variation regularized deconvolution
∑
+ e i,j,t+1 − K
λ3 |K e i,j,t | for efficient low-dose CT perfusion / Ruogu Fang,
i, j, t
Pina C Sanelli, Shaoting Zhang, Tsuhan Chen
∑ // International Conference on Medical Image
+ e i,j,t+1 − 2 K
λ4 |K e i,j,t + K
e i,j,t−1 |.
Computing and Computer-Assisted Intervention /
i, j, t Springer. —2014. — P. 154–161.
(14) [5] Robust low-dose CT perfusion deconvolution
where λ = (λ1 , λ2 , λ3 , λ4 ) is a regularization parame- via tensor total-variation regularization / Ruogu
ter. Fang, Shaoting Zhang, Tsuhan Chen, Pina C
Sanelli // IEEE transactions on medical imaging.
4. Optimization algorithm — 2015. — Vol. 34, no. 7. — P. 1533–1548.
We use Nesterov accelerated gradient descent [9] [6] Kadimesetty V. S. et al. Convolutional neural
for functional minimization (9): network-based robust denoising of low-dose com-
puted tomography perfusion maps // IEEE Trans-
z0 = 0,
actions on Radiation and Plasma Medical Sci-
y0 = 0, ences. – 2018. – Vol. 3. – no. 2. – P. 137-152.
(15)
zk+1 = yk − ϵ · ∇F (yk ), [7] Golub Gene H., Reinsch Christian. Singular value
yk+1 = zk + βk (zk+1 − zk ), decomposition and least squares solutions // Nu-
merische mathematik. — 1970. — Vol. 14, no. 5.
where βk = 1 − 3 / (k + 1), and ϵ is the learning rate. — P. 403–420.
In this paper we used ϵ = 10−8 . [8] Nasonov A., Krylov A. An improvement of BM3D
image denoising and deblurring algorithm by gen-
5. Method testing
eralized total variation // 2018 7th European
Clinical perfusion data does not have ground truth Workshop on Visual Information Processing (EU-
values of residue function k(t) and perfusion parame- VIP). — 2018. — P. 1–4.
ters. Therefore, synthetic data was generated. As a [9] Bottou L., Curtis F. E., Nocedal J. Optimization
base for generation we take perfusion maps for phan- methods for large-scale machine learning // Siam
tom [10]. Then we generate function k(t) = C ·e−(at) .
2
Review. – 2018. – Vol. 60. – no. 2. – P. 223-311.
0.004 0.004 0.004 0.004
0.003 0.003 0.003 0.003
Synthetic
0.002 0.002 0.002 0.002
0.001 0.001 0.001 0.001
0.000 0.000 0.000 0.000
0.001 0.001 0.001 0.001
0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50
0.015 0.015 0.015 0.015
0.010 0.010 0.010 0.010
Clinical
0.005 0.005 0.005 0.005
0.000 0.000 0.000 0.000
0.005 0.005 0.005 0.005
0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35
SVD TV TGV Ground truth
Fig. 4. k(t) for a brain tissue point.
(a) Ground truth (b) SVD (c) TV (d) TGV
Fig. 5. MTT results for synthetic data by different methods.
[10] Digital brain perfusion phantom. —
https://www5.cs.fau.de/research/data/digital-
brain-perfusion-phantom/.