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/.