Automatic Delineation of Kidney Region in DCE-MRI

05/26/2019 ∙ by Santosh Tirunagari, et al. ∙ University of Surrey 0

Delineation of the kidney region in dynamic contrast-enhanced magnetic resonance Imaging (DCE-MRI) is required during post-acquisition analysis in order to quantify various aspects of renal function, such as filtration and perfusion or blood flow. However, this can be obfuscated by the Partial Volume Effect (PVE), caused due to the mixing of any single voxel with two or more signal intensities from adjacent regions such as liver region and other tissues. To avoid this problem, firstly, a kidney region of interest (ROI) needs to be defined for the analysis. A clinician may choose to select a region avoiding edges where PV mixing is likely to be significant. However, this approach is time-consuming and labour intensive. To address this issue, we present Dynamic Mode Decomposition (DMD) coupled with thresholding and blob analysis as a framework for automatic delineation of the kidney region. This method is first validated on synthetically generated data with ground-truth available and then applied to ten healthy volunteers' kidney DCE-MRI datasets. We found that the result obtained from our proposed framework is comparable to that of a human expert. For example, while our result gives an average Root Mean Square Error (RMSE) of 0.0097, the baseline achieves an average RMSE of 0.1196 across the 10 datasets. As a result, we conclude automatic modelling via DMD framework is a promising approach.



There are no comments yet.


page 3

page 4

page 5

page 6

page 7

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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) (b)
Figure 1: (a) DCE-MR image of a healthy volunteer where the kidney is contrast enhanced. Quantification of kidney function requires computation of mean intensities of the pixel values. However, the quantification is contaminated by the mixing of intensity values from liver region and background artefacts due to complex patient movements. (b) An example of blurring or partial volume effect using a segment of a DCE-MR image. The image on the middle shows a magnified region from the kidney-liver boundary. Clearly there is no definitive boundary between the tissues due to the limited resolution of the image, which indicates that the voxels may contain more than one signal intensity. The blurring effect, however, is not limited to the boundaries such as above since the kidney has a complex structure that consists of different tissues as shown in (a). The image on the right shows the delineation of kidney region by the proposed framework clearly suppressing the mixing from the liver region.

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 ROI 


. In addition, k-means clustering  

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

Figure 2: (a) DMD mode 2 showing the pre-dominant kidney region by suppressing the background and liver regions. (b) Thresholded version of the DMD mode-2. (c) Selection of the kidney ROI using the 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.

2 Methodology

Our methodological framework (Figure 3) consists of DMD, thresholding and selection of kidney region using connected component (blob) analysis.

Figure 3: Flow chart showing the steps involved in the methodological framework. First, a DCE-MRI sequence consisting of images is processed using the DMD algorithm in order to output dynamic mode images. From which, we select a dynamic mode image corresponding to mode =

. Second, thresholding is performed on this mode and converted to binary image. After binarization, the largest area with connected pixels is selected as template. Finally, the produced template is projected onto DCE-MRI data to quantify the kidney function.

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 matrix

P 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:


From Equations 4 and 6, we have

In Schmid2

, 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 obtain

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

3 Dataset

In this section we briefly describe our synthetically generated data as well as the DCE-MRI data.

(a) (b)
(c) (d)
Figure 4: (a) A random image from a DCE-MRI sequence showing the regions of liver (1), background (2) and kidney (3) regions. (b) Graph showing the time-intensity evolution of idealised kidney, liver and background functions for synthetically generated data. (c) Simulated image representing kidney, liver and background functions corresponding to the DCE-MRI image as shown in (a). (d) Gaussian PSF Convolved idealised representation of the image showing the effects of PVC i.e., by mixing up of regions especially near the borders.

3.1 Synthetic data

(1) (2) (3) (4) (5)
(6) (7) (8) (9) (10)
Figure 5: Exemplars of Dynamic MR images from 10 healthy volunteers’ kidney slice produced by DCE-MRI sequence considered as 10 different datasets in this study. The images here show the central kidney slice at time aortic peak enhancement after the contrast agent is injected. The yellow boundary on the kidneys are a results of manual delineation from a human expert. The mean intensity values are calculated in this region across the time producing time-intensity plots.

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

A series of images are then produced representing kidney, liver and background regions/functions as labelled in Figure 4(c) with respect to their temporal evolution shown in Figure 4

(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
Figure 6: (a) Selection of DMD modes 1, 2 and 7. (b) Thresholding and extraction of synthetic regions. (c) Quantification of kidney, liver and background regions showing the partial volume correction (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)).

Figure 7: (Top) The top 10 ranked DMD modes on dataset-1. (Bottom) 10 images showing low ranked DMD modes on dataset-1.

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
Figure 8: (a) Images showing the kidney regions in Dynamic mode-2 across the 10 datasets used in this study. (b) Images showing the thresholding effect on dynamic mode-2 across the datasets. (c) Images showing the kidney template selected by the largest area of connected components.

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.

(1) (2) (3) (4)
(5) (6) (7) (8)
(9) (10)
Figure 9: Quantification of kidney region in DCE-MRI datasets 1 to 10.

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.

Figure 10: RMSE for the time-intensity plots produced by DMD framework and baseline against human expert across 10 datasets.

5 Conclusions

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.

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.


  • (1)

    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
  • (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
  • (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)