DCE-MRI renography involves injection of an MRI contrast agent into the blood stream for obtaining magnetic resonance images (MRI) of internal body structures such as kidney, liver and spleen as shown in Figure 1. The images are taken sequentially over the time, until the contrast agent is discharged. Delineation of the kidney region from these dynamic image sequence is required in order to quantify various aspects of renal function, such as filtration and perfusion or blood flow. However, Absolute quantification in DCE-MRI is often obfuscated by the issues including region of interest (ROI) definition, which may be influenced by the adjacent regions giving rise to the so called Partial Volume Effect (PVE). PVE is the contamination of any single voxel with two or more signal intensities from different tissues or regions (refer Figure 1 for more explanation). The contamination is caused by the finite resolution of imaging systems that limits the higher frequencies i.e., fine details such as edges and structures. This results in producing inaccurate or blurred images especially in DCE-MRI where the spatial information is sacrificed for better temporal resolution since time is an important factor in the assessment of kidney function. The two primary issues concerning PVE are: (i) the action of the Point Spread function on the image data, and (2) the action of finite voxel sampling (comb function) and the integration of signal within the volume of the voxel (rect function). Further details, albeit for the same effects in 2D can be found in Yip et al, Phys Med Biol 2011, listed under my papers on the uni web page. Therefore, DCE-MRI is required to undergo several steps including: 1) suppression of PVE gutierrez2010partial , selection of kidney ROI zollner2009assessment ; rusinek2007performance and a final stage of modelling kidney function by analysis of the resulting time-intensity plots tofts2012precise .
A template-based approach based on tissue classification proposed in gutierrez2010partial
has demonstrated that PVE is a very significant issue in quantification of DCE-MRI studies. From the literature, eigen based methods such as Independent component analysis (ICA) and its variants have also been used for obtaining the kidney ROIzollner2007assessment
. In addition, k-means clusteringzollner2009assessment has also been employed. Although these methods have shown improvements in the processing of the data, a complete automated solution has not been provided. The major issue with the aforementioned approaches is the need for manual or semi-automatic delineation of the target region, and the need for prior knowledge of the point spread function of the particular MRI sequence used, which will vary between particular sequences and vendors. In addition, this can be labour intensive, time-consuming and inefficient as the human expert has to examine the whole sequence of images to find the most suitable frame for delineation. Therefore, to address this problem in more holistic manner, we in this study, present a framework consisting of DMD, thresholding and blob analysis for automatic delineation of kidney ROI with the PVE compensation within the DCE-MRI sequence.
DMD was originally introduced in the area of computational fluid dynamics (CFD) Schmid2 , specifically for analysing the sequential image data generated by non-linear complex fluid flows Schmid1 ; Schmid2 ; Schmid3 ; tirunagarianalysis . The DMD decomposes a given image sequence into several images, called dynamic modes. These modes essentially capture different large scale to small scale structures (sparse components) including a background structure (low rank model) DBLP:journals/corr/GrosekK14 ; randDMD ; tirunagari2016can . DMD has gained significant applications in various fields 6926317 ; 2015arXiv150804487M ; brunton2016extracting , including for detecting spoof samples from facial authentication video datasets tirunagari2015detection and for detecting spoofed finger-vein images tirunagari2015windowed . The advantage of this method is its ability to identify regions of dominant motion in an image sequence in a completely data-driven manner without relying on any prior assumptions about the patterns of behaviour within the data. Therefore, it is thus potentially well-suited to analyse a wide variation of blood flow and filtration patterns seen in renography pathology.
Our preliminary results show that DMD mode-2 highlights the kidney region due to the presence of injected contrast agent. Therefore, as shown in Figure 2, this mode is then used for delineating the kidney region using thresholding and blob analysis.
Prior works directed towards tackling the issues related to Partial Volume Correction (PVE) and delineation of kidney ROI are often dealt separately using different methods. In particular, to our knowledge, there is no literature of using a single method that is able to tackle both of these issues simultaneously. Therefore, our first contribution is to introduce a single methodology that is based on DMD to tackle the aforementioned issues problem. Our second contribution is to Validate our technique for the first time on medical data with applications to DCE-MRI; and finally, our third contribution is in improving the understanding of the application through our DMD framework.
The remainder of this paper is organised as follows: in Section 2, we consider the theory for DMD. Section 3 presents synthetically generated dataset as well as ten different DCE-MRI datasets used in this study. In Section 4 we present our experiments and results and finally, conclusions are drawn in Section 5.
Our methodological framework (Figure 3) consists of DMD, thresholding and selection of kidney region using connected component (blob) analysis.
In our previous studies Tirunagari2017 , we presented a novel automated, registration-free movement correction approach based on windowed and reconstruction variants of Dynamic Mode Decomposition (WR-DMD) to suppress unwanted complex organ motion in DCE-MRI image sequences caused due to respiration.
2.1 Dynamic Mode Decomposition (DMD)
In a dynamic sequence of images P, let be the image whose size is . This image is converted to
column vector, resulting in the construction of a data matrixP of size for images in the DCE-MRI data.
The images in the DCE-MRI data are collected over regularly spaced time intervals and hence each pair of consecutive images are linearly correlated. It can be justified that a linear mapping A exists between them forming a span of krylov subspace krylov1931numerical ; saad1981krylov ; ruhe1984rational :
The krylov subspace can then be represented using two matrices and where and
The mapping matrix is responsible for capturing the dynamics or the fluctuating intensities within the image sequence. The sizes of the matrices and are both each. Therefore, the size of unknown matrix would be . Unfortunately, solving for A is computationally very expensive due to it size. For instance, if an image has a size of i.e., and , the size of is then .
Since solving is computationally expensive for large image dimensions, we need an alternative solution. Our assumption that the images form a Krylov span, allows us to introduce from the standard Arnoldi iteration trefethen1997numerical ).
Here, a companion matrix also known as a shifting matrix that simply shifts images through and approximates the last frame by linearly combining the previous images, i.e., . requires the storage of data matrix is significantly smaller than in dimensions.
Thus, for the last frame , where is much fewer (dimension) than the dimensionality of (), one can write as a linear combination of the previous vectors. Consistent with Equations 3 and 4, we then have:
, the author describes a more robust solution, which is achieved by applying a singular value decomposition (SVD) on. From Equation 6, SVD decomposition on subspace is calculated to obtain , and matrices that are left singular vectors, singular values and right singular vectors respectively. The inversions of these matrices are then multiplied with subspace to obtain the The full-rank matrix , determined on the subspace spanned by the orthogonal basis vectors of . described by:
Here, and are the conjugate transpose of and , respectively; and denote the inverse of the singular values . After obtaining the
matrix, the eigenvalue analysis is performed to obtaineigenvectors and a diagonal matrix containing the corresponding eigenvalues.
It is known that the eigenvalues of approximate some of the eigenvalues of the full system . The associated eigenvectors of provide the coefficients for the linear combination that is necessary to express the dynamics within the image sequence basis. The dynamic modes are thus calculated as follows:
Since the eigenvalues associated with are complex, it is not possible to establish the order of the dynamic modes directly. Nevertheless, we calculate the absolute value for the phase-angles associated with the complex eigenvalues and the modes with unique phase-angles are selected. Doing this will remove one of the conjugate pairs in the dynamic modes which have same phase-angles but with different signs and capture similar information sayadi2013dynamic . After discarding one of the conjugate pairs, the dynamic modes are then sorted in the ascending order of their phase-angles.
In this section we briefly describe our synthetically generated data as well as the DCE-MRI data.
3.1 Synthetic data
In order to validate our experiments with the ground-truth available, we artificially generated synthetic data corresponding to the DCE-MRI events as that shown in Figure 4(c). That is, as the liver region denoted by in DCE-MRI image (Figure 4(a)) covers both background region () and kidney region (), similarly the synthetically generated images follow the same structure as shown in Figure 4(c).
(b). The kidney region’s temporal evolution is thus modelled as a combination of Poisson and log distributions. The liver region is modelled as a sigmoid distribution and the background region as a random normal distribution. The generated regions are then convoluted using a Gaussian kernel operator (variance = 22 and block size of). This operator thus simulates the point spread function’s (PSF’s) pixel-level mixing of the kidney regions as shown in Figure 4(d).
The time-intensity plot of convolved kidney region thus shows the influence of the liver (sigmoid) and background functions as shown in Figure 6. Similarly the time-intensity plots of liver are influenced by background and kidney regions and vice versa.
3.2 DCE-MRI data
The functional kidney DCE-MRI datasets used in this study have been obtained from ten healthy volunteers as shown in Figure 5. These datasets are acquired after injection of of Gd-DTPA (Magnevist) contrast agent, on a Siemens Avanto MRI scanner, using a 32 channel body phased array coil. The MRI acquisition sequence consisted of a 3D spoilt gradient echo sequence utilising an Echo time, TE of 0.6ms, Repetition time, TR of 1.6ms, and a flip angle, FA = 17 degrees with a temporal resolution of collected for . The acquired DCE-MRI datasets cover the abdominal region, enclosing left and right kidneys and abdominal aorta.
4 Experiments & Results
In this section, we present our experimental procedure as well as the results.
4.1 Evaluating the DMD framework over synthetically generated data
Synthetic data consisting of images is given as an input to DMD algorithm which produces DMD modes (In theory, for a image sequence with images, we obtain DMD modes). Recall that each dynamic mode captures a “principal dynamics” axes of the image sequence. Modes that show pre-dominating liver, kidney and background functions are selected, as shown in Figure 6 (a). These modes are then thresholded to obtain the segmented versions of the respective functions as shown in Figure 6 (c). These segmented versions of the dynamic modes are then projected on to the convolved image sequence and mean pixel intensities of synthetic kidney, liver and background functions are computed. The time-intensity plots are shown in Figure 6 (c). It is clear that the mixing from the other functions and are totally suppressed and the time-intensity plots produced by the proposed framework recovers the original functions.
|(a) Dynamic modes||(b) Thresholded||(c) PVC|
4.2 Evaluating the DMD framework on DCE-MRI
DCE-MRI image sequences are given as an input to the DMD algorithm. The output of the DMD captures large (Figure 7(Top)) and small scale structures (Figure 7(bottom)) in the form of modes. Dynamic mode-1 reveals the background (low-rank) model and the remaining modes capture the sparse representations. The contrast changes are captured in the top modes, particularly, mode-2 capturing kidney region and mode-3 and 4, the spleen and the liver regions respectively, as shown in Figure 7(top) on the dataset-1. Noise and residuals are captured in the lower modes (Figure 7(bottom)).
Therefore, we have selected DMD mode-2 across all the 10 datasets as shown in Figure 8 (a). These modes are then thresholded to obtain a binary of version (Figure 8 (b)). Later using connected component analysis on the left half of the image, the largest blob is selected as a kidney ROI as shown in Figure 8 (c). The mean of intensity values inside this automatically delineated kidney ROI are calculated to produce time-intensity plots.
|(a) Dynamic mode-2 across 10 datasets of DCE-MRI data|
|Thresholding across 10 datasets of DCE-MRI data|
|Delineated kidney region obtained via blob analysis across 10 datasets of DCE-MRI data|
4.3 Comparison with human expert
As a baseline method a minimum bounding box region around the kidney region (3) is selected as shown in Figure 4(a) across all the datasets. To evaluate the performance of the DMD framework we compare it against the performance of the manual delineation from a human expert. From Figure 9 we see that the kidney function produced from the automatic delineation clearly compares favourably with the human expert; where as the minimum bounding box region due to the influence of the PVE deviates from the kidney function quantified by a human expert.
The root mean square error (RMSE) between the quantification produced by the proposed framework and baseline method against the human expert is shown in Figure 10. It is clear that the quantification results obtained from the proposed framework has lower RMSE when compared to the baseline method.
This study shows the significance of the proposed framework based on DMD, thresholding and blob analysis as a viable automatic delineation algorithm to effectively quantify the kidney function by reducing the partial volume effect in the DCE-MRI data.
Our proposed DMD framework is applied on the DCE-MRI datasets collected from 10 healthy volunteers as well as the synthetic data mimicking the DCE-MRI data. We found that DMD can capture functional structures like kidney, liver and spleen in the top ranked dynamic modes. Particularly dynamic mode-2 capturing the kidney structure across all the datasets. This mode-2 is then thresholded and kidney template is obtained using connected component analysis. The kidney function produced from our proposed framework compares favourably with the human expert; while the minimum bounding box region due to the influence of the PVE deviates. Our results demonstrate that our proposed framework is very promising in delineating and quantifying the kidneys in DCE-MRI studies by removing the PVE.
Acknowledgements.The funding for this work has been provided by the Department of Computer Science and the Centre for Vision, Speech and Signal Processing (CVSSP) - University of Surrey. We also would like to express our gratitude towards Kidney Research UK for funding the DCE-MRI data acquisition as part of a reproducibility study.
Berger, E., Sastuba, M., Vogt, D., Jung, B., Amor, H.B.: Dynamic mode decomposition for perturbation estimation in human robot interaction.In: The 23rd IEEE International Symposium on Robot and Human Interactive Communication, pp. 593–600 (2014). DOI 10.1109/ROMAN.2014.6926317
- (2) Brunton, B.W., Johnson, L.A., Ojemann, J.G., Kutz, J.N.: Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition. Journal of neuroscience methods 258, 1–15 (2016)
- (3) Erichson, N., Donovan, C.: Randomized low-rank dynamic mode decomposition for motion detection. Computer Vision and Image Understanding In Press (2016). DOI 10.1016/j.cviu.2016.02.005. N. Benjamin Erichson acknowledges support from the UK Engineering and Physical Sciences Research Council (EPSRC).
- (4) Grosek, J., Kutz, J.N.: Dynamic mode decomposition for real-time background/foreground separation in video. CoRR abs/1404.7592 (2014). URL http://arxiv.org/abs/1404.7592
- (5) Gutierrez, D.R., Wells, K., Diaz Montesdeoca, O., Moran Santana, A., Mendichovszky, I., Gordon, I.: Partial volume effects in dynamic contrast magnetic resonance renal studies. European journal of radiology 75(2), 221–229 (2010)
- (6) Krylov, A.: On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined. Izvestija AN SSSR (News of Academy of Sciences of the USSR), Otdel. mat. i estest. nauk 7(4), 491–539 (1931)
- (7) Mann, J., Kutz, J.N.: Dynamic Mode Decomposition for Financial Trading Strategies. ArXiv e-prints (2015)
- (8) Ruhe, A.: Rational krylov sequence methods for eigenvalue computation. Linear Algebra and its Applications 58, 391–405 (1984)
- (9) Rusinek, H., Boykov, Y., Kaur, M., Wong, S., Bokacheva, L., Sajous, J.B., Huang, A.J., Heller, S., Lee, V.S.: Performance of an automated segmentation algorithm for 3d mr renography. Magnetic Resonance in Medicine 57(6), 1159–1167 (2007)
- (10) Saad, Y.: Krylov subspace methods for solving large unsymmetric linear systems. Mathematics of computation 37(155), 105–126 (1981)
- (11) Sayadi, T., Schmid, P., Nichols, J., Moin, P.: Dynamic mode decomposition of controlled h-and k-type transitions. Annual Research Briefs (Center for Turbulence Research, 2013) p. 189 (2013)
- (12) Schmid, P., Li, L., Juniper, M., Pust, O.: Applications of the dynamic mode decomposition. Theoretical and Computational Fluid Dynamics 25(1-4), 249–259 (2011)
- (13) Schmid, P.J.: Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, 5–28 (2010)
- (14) Schmid, P.J., Meyer, K.E., Pust, O.: Dynamic mode decomposition and proper orthogonal decomposition of flow in a lid-driven cylindrical cavity. In: 8th International Symposium on Particle Image Velocimetry, pp. 25–28 (2009)
- (15) Tirunagari, S., Poh, N., Bober, M., Windridge, D.: Windowed dmd as a microtexture descriptor for finger vein counter-spoofing in biometrics. In: Information Forensics and Security (WIFS), 2015 IEEE International Workshop on. IEEE (2015)
- (16) Tirunagari, S., Poh, N., Bober, M., Windridge, D.: Can dmd obtain a scene background in color? In: Image, Vision and Computing (ICIVC), International Conference on, pp. 46–50. IEEE (2016)
- (17) Tirunagari, S., Poh, N., Wells, K., Bober, M., Gorden, I., Windridge, D.: Movement correction in dce-mri through windowed and reconstruction dynamic mode decomposition. Machine Vision and Applications 28(3), 393–407 (2017). DOI 10.1007/s00138-017-0835-5. URL http://dx.doi.org/10.1007/s00138-017-0835-5
- (18) Tirunagari, S., Poh, N., Windridge, D., Iorliam, A., Suki, N., Ho, A.T.: Detection of face spoofing using visual dynamics. Information Forensics and Security, IEEE Transactions on 10(4), 762–777 (2015)
- (19) Tirunagari, S., Vuorinen, V., Kaario, O., Larmi, M.: Analysis of proper orthogonal decomposition and dynamic mode decomposition on les of subsonic jets. CSI Journal of Computing 1, 20–26 (2012)
- (20) Tofts, P.S., Cutajar, M., Mendichovszky, I.A., Peters, A.M., Gordon, I.: Precise measurement of renal filtration and vascular parameters using a two-compartment model for dynamic contrast-enhanced mri of the kidney gives realistic normal values. European radiology 22(6), 1320–1330 (2012)
- (21) Trefethen, L.N., Bau III, D.: Numerical linear algebra, vol. 50. Siam (1997)
- (22) Zöllner, F.G., Kocinski, M., Lundervold, A., Rørvik, J.: Assessment of renal function from 3d dynamic contrast enhanced mr images using independent component analysis. In: Bildverarbeitung für die Medizin 2007, pp. 237–241. Springer (2007)
- (23) Zöllner, F.G., Sance, R., Rogelj, P., Ledesma-Carbayo, M.J., Rørvik, J., Santos, A., Lundervold, A.: Assessment of 3d dce-mri of the kidneys using non-rigid image registration and segmentation of voxel time courses. Computerized Medical Imaging and Graphics 33(3), 171–181 (2009)