Printer Friendly

The Time-Domain Integration Method of Digital Subtraction Angiography Images.

1. Introduction

Circulatory system diseases, such as aortic dissection, have been the focus of medical research [1, 2] for their dangerousness and high incidence. Improving the clarity of the captured medical image helps us to diagnose more accurately, which is an important research field.

Digital subtraction angiography (DSA), as a real-time approach, is commonly employed in the clinical diagnosis of circulatory system disease [3, 4], especially in the realtime surgical monitoring and the medical examination among small branches of blood vessels which are difficult to be measured by other methods. In order to protect the patient, it is important to shorten the shooting time and to reduce the dosage of contrast media when capturing the images.

A lot of image denoising methods have been proposed in recent years, but they are all problematic when applied to the DSA images. For example, the image reconstruction method based on the level set theory [5], wavelet decomposition and reconstruction method [6, 7], Bayesian method [8], and image denoising method based on anisotropic diffusion [9] generally need a long operation time, which cannot meet the real-time requirements of DSA image processing. Moreover, images processed by these approaches are usually not clear enough to show the details such as edges and textures. In 2004, Candes et al. [10] proposed an image denoising method based on the sparse decomposition. On this basis, Needell and Vershynin [11] proposed the regularized orthogonal matching pursuit (ROMP) method; Scholefield and Dragotti [12] used a sparse quadtree decomposition representation to remove the noise in images; Adler et al. employed the shrinkage learning approach to acquire the high-resolution reconstruction image [13-17]. However, the operation of these approaches is also very complicated. Therefore, it is necessary to find an image processing method which is more suitable for the real-time analysis of DSA images.

In this work, an image time-domain integration method based on blood flow periodicity has been proposed. In this algorithm, the DSA images of the first cardiac cycle after the injection of the contrast agent are extracted denoised by the wavelet reconstruction method firstly, and then these images are integrated to obtain the time-domain integration image, which is named after the TDI image in this paper. This method contributes to the diagnosis of circulatory system diseases.

2. Materials and Methods

2.1. The Noise Model. The theoretical gray-scale [w.sub.j]([x.sub.0],[y.sub.0]) at a certain pixel ([x.sub.0], [y.sub.0]) on the j-th frame of DSA images can be obtained from the following equation by the Lambert-Beer law [18, 19]:

[mathematical expression not reproducible] (1a)

where [I.sub.0] and [I.sub.r] are the X-ray transmission amount before and after the addition of the contrast agent, respectively. [V.sub.j]([x.sub.0], [y.sub.0]) is the volume of blood vessels at pixel ([x.sub.0], [y.sub.0]). Nj ([x.sub.0], [y.sub.0]) and c ([V.sub.j]) are the number and amount of substance concentration of contrast agent particles at pixel ([x.sub.0], [y.sub.0]), respectively. k donates the absorption coefficient.

The image quality degrades in the original DSA image, as a result of the limitations on the imaging system's resolution and the influence of additive noise such as Gaussian noise, which is donated by [w.sub.j] and can be expressed in as follows [6]:

[w.sub.j] = [w.sub.j] x [h.sub.j] + [n.sub.j], (1b)

where [h.sub.j] and [n.sub.j] represent the point spread function and the additive noise, respectively. Operator "*" is the convolution operator. Equation (1b) can be rewritten to matrix form using the block Toeplitz matrix [H.sub.j], as shown in the following equation:

[W.sub.j] = [w.sub.j] x [H.sub.j] + [n.sub.j]. (1c)

It is difficult to solve Equation (1c) when only [W.sub.j] is given. However, since the gray-scale level of a certain pixel is proportional to the number of contrast agent particles in that pixel, the number of contrast agent particles follows the motion pattern of blood. And as for the blood motion pattern, on consideration of the periodicity of human heartbeat, the blood flow rate in human body is also cyclical, which can be expressed in the following equation:

[V.sub.b](nT + t)= [bar.[v.sub.b]](t)+ [PSI]([bar.[V.sub.b]), t [member of] Z, (2)

where [V.sub.b] (nT +1) is the velocity field of blood at time (nT + t). [bar.[V.sub.b]](t) denotes the average velocity field of blood at time t, which is the mean flow velocity at the same time in multiple cardiac cycles. T is the length of the cardiac cycle. [PSI]([bar.[V.sub.b]]) characterizes the changes in flow velocity owing to factors such as the instability of human blood pressure. [PSI]([bar.[V.sub.b]) can be regarded as a zero-mean-value distribution with a small variance, since patients are under the total anesthesia during the shooting process and their vital signs remain stable. According to the Wilke-Chang equation [20], the free diffusion rate of the contrast agent in the blood is much smaller than the blood flow rate, and thus the contrast agent obeys the same movement law as the blood. Therefore, the periodicity of the blood flow rate can be employed to improve the clarity of DSA images.

2.2. Image Integration. In order to decrease the shoot time and the contrast agent's injection quantity, images in the first cardiac cycle (donated by the cardiac cycle S) after the injection of the contrast agent are analyzed in this study. Firstly, the time at which the cardiac cycle begins is set to be t = 0. Subsequently, the velocity of the i-th contrast agent particle in this cardiac cycle is expressed as [u.sub.i](t) = ([u.sub.i](t) x x, [u.sub.i](t) x y, [u.sub.i](t) x z). Since the motion of the contrast agent particles is consistent with that of the blood, and on consideration of the velocity stability shown in Equation (2), the velocity of the particle at each position on its trajectory can be regarded as a sample of the blood flow field at that location. Therefore, once the substantial number of particles is extracted, the average velocity field of the contrast agent at the pixel ([x.sub.0], [y.sub.0]) in the entire cardiac cycle can be estimated by the mean velocity field value of particles which flowed through that pixel during the calculated cardiac cycle, as shown in the following equation:

[bar.[v.sub.c]]([x.sub.0],[y.sub.0]) = 1/N([x.sub.0],[y.sub.0]) [N([x.sub.0],[y.sub.0]).summation over (i=1)][u.sub.i]([t.sub.i]), (3a)

where [t.sub.i] represents the time when the i-th particle approached pixel ([x.sub.0], [y.sub.0]).

The total time length that the i-th particle appears in pixel ([x.sub.0], [y.sub.0]) during one cardiac cycle T satisfies Equation (3b), where [L.sub.i]([x.sub.0], [y.sub.0]) represents the distance of the i-th particle in the range of pixel ([x.sub.0], [y.sub.0]) and [u.sub.i,//]([x.sub.0], [y.sub.0]) is set to be the magnitude of velocity component which is parallel to the image plane in pixel ([x.sub.0], [y.sub.0]) during that cardiac cycle since the photographing each DSA image can be regarded as a sample of each particle's location:

[T.sub.i]([x.sub.0],[y.sub.0]) = min{[L.sub.i]([x.sub.0],[y.sub.0]/[u.sub.i,//]([x.sub.0], [y.sub.0]), T}. (3b)

A variate [lambda] is set to represent the absorption capacity of light in the unit time of a single contrast agent particle. After that, the time-weighted gray-scale value of the i-th particle at pixel ([x.sub.0], [y.sub.0]), [p.sub.i]([x.sub.0], [y.sub.0]), can be expressed by the following equation:

[P.sub.i]{[x.sub.0],[y.sub.0]) = [lambda] x [T.sub.i]{[x.sub.0],[y.sub.0]). (3c)

On combination of Equations (3a)-(3c), [bar.[v.sub.b.//]]([x.sub.0], [y.sub.0]) can be characterized by the sum of the time integral intensities of particles which have appeared in pixel ([x.sub.0], [y.sub.0]) during the cardiac cycle, as shown in the following equation:

[mathematical expression not reproducible] (3d)

where [bar.L] is the average moving distance of the contrast agent particles within that pixel. Since the size of a pixel is small, [bar.L] is approximately equal to the length of each pixel, [L.sub.pixel]. According to Equations (1a) and (3d), when the frame rate M tends to infinity, the following equation can be obtained:

[mathematical expression not reproducible] (3e)

Equation (3e) demonstrates that the overall time-domain integration value of pixel ([x.sub.0], [y.sub.0]), b([x.sub.0], [y.sub.0]), can be expressed as the integral of each picture's gray value at that position in the entire cardiac cycle. On consideration that the time step [DELTA]t is short in the actual case, Equation (3e) can be employed in the calculation of the captured DSA images. Therefore, the relationship shown in Equation (3f) can be established. And the image b is named after the time-domain integration image or the TDI image:

[mathematical expression not reproducible]. (3f)

Furthermore, to strength the denoising effect, the images are denoised by the median filter before the time-domain integration since Ling's work [21] shows that the noiseless image is usually insensitive to a median filter. Table 1 shows the specific steps of our method.

3. Results

Figure 1 shows a group of DSA images from a patient with aortic dissection. Figures 1(a1)-1(c1) are the DSA images at t = 1/6T, 1/2T, and 5/6T, respectively. Figure 1(d1) is the TDI image of the dissecting aneurysm extracted by our algorithm. For the sake of comparison, the image of dissecting aneurysm region in Figures 1(a1)-1(c1) are extracted and then the normalized gray-scale histograms of dissecting aneurysm regions among Figures 1(a1)--1(d1) are obtained, which are shown in Figures 1(a2)-1(d2), respectively. Table 2 shows the mean value, standard deviation, and coefficient of variation of DSA images of the dissecting aneurysm region in one cardiac cycle of the patient in Figure 1.

Table 2 shows that the coefficient of variation of the TDI image is higher than all the DSA images in that cardiac cycle, which means that our TDI image can enhance the details. Moreover, the gap between the peaks in Figure 1(d2) is clearer than those in Figures 1(a2)-1(c2), which means that our TDI image has higher resolution than the original images.

4. Discussion

Since the blood flow velocity has a certain degree of uncertainty in the actual situation, which directly influenced the gray-scale value of the shot DSA image, Equation (1c) can be rewritten as follows:

[W.sub.j] - [w.sub.j] x [H.sub.i] + [n.sub.j] = (([q.sub.j] + [[PHI].sub.j]) x [H.sub.j]) + [n.sub.j], (4)

where [q.sub.j] donates the gray-scale value after the removal of motion randomness.

Equations (3f) and (4) illustrate that the time-domain integration image b can be obtained by the following equation:

[mathematical expression not reproducible]. (5a)

where matrix [G.sub.j] is defined to characterize the differences between [W.sub.j] and [q.sub.j]. A new symbol [b.sub.x] is defined here, which represents the noise-free time-domain integration image. [b.sub.x] = [[summation].sup.M-1.sub.j=0][q.sub.j]. Then, the error analysis of Equation (5a) is implemented below.

[mathematical expression not reproducible], (5b)

[mathematical expression not reproducible], (5c)

[mathematical expression not reproducible]. (5d)

According to Reference [22], when the cardiac cycle length is 0.8 seconds, the blood flow rate at the aortic inlet, denoted by [v.sub.inlet], follows Equation (5b). The first part of Equation (5b) denotes the ejection phase, which is followed by a brief closure of blood after the closure of the aortic valve. The flow rate in the rest time is 0. Since the time varies when each particle enters the view field, the time when they arrive at the same position on the image also varies. Therefore, the number of images which satisfy [N.sub.j](x, y) [not equal to] 0 is greater than one in most positions. Thus, Equation (5c) can be obtained.

According to Equation (5c), the signal-to-noise ratio of image b ([b.sub.x] is defined as the noise-free image in this calculation) is higher than the mean signal-to-noise ratio of original DSA images at each pixel. Therefore, the result of our method contains lower noise in the entire view field, and our method can suppress the noise.

Furthermore, Equation (5d) can be deduced from Equation (3e), where [[parallel][parallel].sub.2] stands for the two-norm of matrix. According to Equation (5d), image b is closer to the noisefree TDI image [b.sub.x] than all of the DSA images, which means that image b has the highest clarity. Thus, the image b can improve the clarity of the original DSA images.

5. Conclusion

In summary, this study presents a DSA image denoising and enhancement method based on the periodicity of blood flow. Firstly, the DSA images are reconstructed through the median filter, and then DSA images in a cardiac cycle are integrated and the overall time-domain integration image b is obtained. According to the mathematical derivation as well as the verification using aortic dissecting aneurysm images, this study demonstrates that the TDI image of the contrast agent has a lower overall noise than original DSA images, and it is also clearer than the original image. This method can contribute to the feature location extraction and disease prophylaxis in circulatory disease, such as the first break's position extraction of aortic dissecting aneurysms and the analysis of stress distribution on vessel wall [23], which are also our future work.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors' Contributions

Shuo Huang and Le Cheng contributed equally to this work.


This work was supported by the National Key R&D Program of China under grant number 2017YFC0112801, the National Natural Science Foundation of China under grant numbers 61127002, 11572087, 61771130, and 31300780, the Research Center for Learning Science of Southeast University under grant number 3207038391, and the School of Biological Sciences and Medical Engineering under grant number 3207037434.


[1] R. Erbel, F. Alfonso, C. Boileau et al., "Diagnosis and management of aortic dissection: recommendations of the task force on aortic dissection, European society of cardiology," European Heart Journal, vol. 22, no. 18, pp. 1642-1681, 2001.

[2] A. White, J. Broder, J. Mando-Vandrick, J. Wendell, and J. Crowe, "Acute aortic emergencies-part 2: aortic dissections," Advanced Emergency Nursing Journal, vol. 35, no. 1, pp. 28-52, 2013.

[3] J. O. Balzer, M. Doss, A. Thalhammer, H.-G. Fieguth, A. Moritz, and T. J. Vogl, "Urgent thoracic aortal dissection and aneurysm: treatment with stent-graft implantation in an angiographic suite," European Radiology, vol. 13, no. 10, pp. 2249-2258, 2003.

[4] E. Sueyoshi, H. Nagayama, I. Sakamoto, and M. Uetani, "Carbon dioxide digital subtraction angiography as an option for detection of endoleaks in endovascular abdominal aortic aneurysm repair procedure," Journal of Vascular Surgery, vol. 61, no. 2, pp. 298-303, 2015.

[5] S. Osher and J. A. Sethian, "Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations," Journal of Computational Physics, vol. 79, no. 1, pp. 12-49, 1988.

[6] S. Wan, B. I. Raju, and M. A. Srinivasan, "Robust deconvolution of high-frequency ultrasound images using higher-order spectral analysis and wavelets," IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 50, no. 10, pp. 1286-1295, 2003.

[7] S. G. Mallat, "Multifrequency channel decompositions of images and wavelet models," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 12, pp. 2091-2110, 1989.

[8] M. Lebrun, A. Buades, and J. M. Morel, "A nonlocal bayesian image denoising algorithm," SIAM Journal on Imaging Sciences, vol. 6, no. 3, pp. 1665-1688, 2013.

[9] P. Perona and J. Malik, "Scale-space and edge detection using anisotropic diffusion," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, no. 7, pp. 629-639, 1990.

[10] E. J. Candes, J. Romberg, and T. Tao, "Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information," IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489-509, 2004.

[11] D. Needell and R. Vershynin, "Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit," Foundations of Computational Mathematics, vol. 9, no. 3, pp. 317-334, 2009.

[12] A. Scholefield and P. L. Dragotti, "Image restoration using a sparse quadtree decomposition representation," in Proceedings of 2009 IEEE International Conference on Image Processing, pp. 1473-1476, IEEE, Cairo, Egypt, November 2009.

[13] A. Adler, Y. Hel-Or, and M. Elad, "A shrinkage learning approach for single image super-resolution with overcomplete representations," in Proceedings of 2010 European Conference on Computer Vision, pp. 622-635, Springer Berlin Heidelberg, Crete, Greece, September 2010.

[14] M. Elad, M. A. T. Figueiredo, and Y. Ma, "On the role of sparse and redundant representations in image processing," Proceedings of the IEEE, vol. 98, no. 6, pp. 972-982, 2010.

[15] D. Needell and R. Vershynin, "Signal recovery from incomplete and inaccurate measurements via regularized orthogonal matching pursuit," IEEE Journal of Selected Topics in Signal Processing, vol. 4, no. 2, pp. 310-316, 2010.

[16] J. L. Starck, J. M. Fadili, and C. Chesneau, "Image deconvolution by stein block thresholding," in Proceedings of 2009 IEEE International Conference on Image Processing, pp. 1329-1332, Cairo, Egypt, November 2009.

[17] D. Tzikas, A. Likas, and N. Galatsanos, "Variational bayesian blind image deconvolution with student-T priors," in Proceedings of 2007 IEEE International Conference on Image Processing, pp. I-109-I-112, San Antonio, TX, USA, September 2007.

[18] M. Elizabeth, Concise Medical Dictionary, Oxford University Press, Oxford, UK, 2015.

[19] W. D. Jeans, "The development and use of digital subtraction angiography," British Journal of Radiology, vol. 63, no. 747, pp. 161-168, 1990.

[20] C. R. Wilke and P. Chang, "Correlation of diffusion coefficients in dilute solutions," AIChE Journal, vol. 1, no. 2, pp. 264-270, 1955.

[21] J. Ling and A. C. Bovik, "Smoothing low-SNR molecular images via anisotropic median-diffusion," IEEE Transactions on Medical Imaging, vol. 21, no. 4, pp. 377-384, 2002.

[22] C. Y. Wen, A. L. Yang, and J. W. Chai, "Investigation of pulsatile flow field in healthy thoracic aorta models," Annals of Biomedical Engineering, vol. 38, no. 2, pp. 391-402, 2010.

[23] J. Charonko, S. Karri, J. Schmieg, S. Prabhu, and P. Vlachos, "In vitro comparison of the effect of stent configuration on wall shear stress using time-resolved particle image velocimetry," Annals of Biomedical Engineering, vol. 38, no. 3, pp. 889-902, 2010.

Shuo Huang (iD),(1,2) Le Cheng, (3) Bin Zhu, (3) Ping Zhou, (1) Yu Sun (iD),(1,4) Bing Zhang (iD),(3) and Suiren Wan (1)

(1) School of Biological Sciences and Medical Engineering, Southeast University, Nanjing 210096, China

(2) Shanghai United Imaging Healthcare Co., Ltd., Shanghai 201807, China

(3) Department of Radiology, The Affiliated Drum Tower Hospital of Nanjing University Medical School, Nanjing 210008, China

(4) Institute of Cancer and Genomic Science, University of Birmingham, Birmingham B15 2TT, UK

Correspondence should be addressed to Yu Sun;

Received 21 May 2018; Accepted 28 August 2018; Published 30 September 2018

Academic Editor: Dominique J. Monlezun

Caption: Figure 1: A set of DSA images from a patient.
Table 1: The image time-domain iteration method.

1. Delineate the region-of-interest

2. Extract the images in the first cardiac cycle after the contrast
agent enters the region-of-interest

3. Image denoising with the median filter

4. Extract the time-weighted image, [p.sub.i] : [p.sub.i] = [W.sub.j]
x [DELTA][t.sub.j]

5. Image time-domain integration, b = [[summation].sup.M-1.sub.j=0]

6. Extract the TDI image in the region-of-interest

Table 2: The mean value, standard deviation, and coefficient of
variation of DSA images of the aortic dissection for the patient in
Figure 1 in the first cardiac cycle after the contrast agent enters
the aortic dissection.

Time                              0       1/6T     1/3T     1/2T

Mean                            0.6794   0.6340   0.6948   0.6507
Standard deviation              0.3585   0.3452   0.3473   0.3554
Coefficient of variation (%)    52.77    54.44    49.98    54.62

Time                             2/3T     5/6T      T      TDI image

Mean                            0.6958   0.5965   0.5978    0.6025
Standard deviation              0.3430   0.3784   0.3835    0.4425
Coefficient of variation (%)    49.30    63.44    64.15      73.44
COPYRIGHT 2018 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2018 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Huang, Shuo; Cheng, Le; Zhu, Bin; Zhou, Ping; Sun, Yu; Zhang, Bing; Wan, Suiren
Publication:Computational and Mathematical Methods in Medicine
Date:Jan 1, 2018
Previous Article:Stimulation Strategies for Tinnitus Suppression in a Neuron Model.
Next Article:Mathematical Modelling of Human African Trypanosomiasis Using Control Measures.

Terms of use | Privacy policy | Copyright © 2021 Farlex, Inc. | Feedback | For webmasters