Automatic Delineation of Kidney Region in DCE-MRI

TL;DR

This paper introduces an automatic method using Dynamic Mode Decomposition (DMD) for delineating kidney regions in DCE-MRI, reducing manual effort and improving accuracy over traditional methods.

Contribution

The paper presents a novel DMD-based framework for automatic kidney delineation in DCE-MRI, validated on synthetic and real datasets, outperforming baseline methods.

Findings

DMD-based method achieves low RMSE of 0.0097

Method performs comparably to human experts

Validated on synthetic and real datasets

Abstract

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.

Equations

$$\text{P}=[p_1,p_2,p_3,\cdots ,p_N]=\left(\begin{array}{cccc}{p}_1^1& p_2^1& ...& p_N^1\\ ⋮& ⋮& ⋮& ⋮\\ p_1^{mn}& p_2^{mn}& ...& p_N^{mn}\end{array}\right).$$

$$\text{P}=[p_1,A{p}_1,A^2{p}_1,A^3{p}_1,\cdots ,A^{N-1}{p}_1].$$

$$P_2=A{P}_1.$$

$$A{P}_1\approx P_{1}H.$$

$$H=\left(\begin{array}[]{ccccc}0&0&\ldots&0&c_{0}\\ 1&0&\ldots&0&c_{1}\\ &\ddots&\ddots&\vdots&\vdots\\ &&1&0&c_{N-2}\\ &&&1&c_{N-1}\end{array}\right).$$

$$P_2\approx P_{1}H.$$

$$\tilde{H}=U^{\ast }{P}_{2}V{\Sigma }^{-1}.$$

$$\tilde{H}\omega =\sigma \omega$$

$$\Psi =P_{2}V{\Sigma }^{-1}\omega$$

Full text

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

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.

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 $N$ images P, let $p_{r}$ be the $r^{th}$ image whose size is $m\times n$. This image $p_{r}$ is converted to $mn\times 1$ column vector, resulting in the construction of a data matrix P of size $mn\times N$ for $N$ images in the DCE-MRI data.

[TABLE]

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 :

[TABLE]

The krylov subspace can then be represented using two matrices $P_{2}$ and $P_{1}$ where $P_{2}\equiv[p_{2},p_{3},\cdots,p_{N}]$ and

$P_{1}\equiv[p_{1},p_{2},\cdots,p_{N-1}]$.

[TABLE]

The mapping matrix $A$ is responsible for capturing the dynamics or the fluctuating intensities within the image sequence. The sizes of the matrices $P_{2}$ and $P_{1}$ are both $mn\times N-1$ each. Therefore, the size of unknown matrix $A$ would be $mn\times mn$. Unfortunately, solving for A is computationally very expensive due to it size. For instance, if an image has a size of $240\times 320$ i.e., $m=240$ and $n=320$, the size of $A$ is then $76800\times 76800$.

Since solving $A$ 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 $H$ from the standard Arnoldi iteration trefethen1997numerical ).

[TABLE]

Here, $H$ a companion matrix also known as a shifting matrix that simply shifts images $1$ through $N-1$ and approximates the last frame $N$ by linearly combining the previous $N-1$ images, i.e., $P_{N}=c_{0}p_{1}+...+c_{N}p_{N-1}=[p_{1},p_{2},p_{3},\cdots,p_{N-1}]c$. $H$ requires the storage of $N-1\times N-1$ data matrix is significantly smaller than $A$ in dimensions.

[TABLE]

Thus, for the last frame $N$, where $N$ is much fewer (dimension) than the dimensionality of $A$ ($mn\times mn$), one can write $P_{2}$ as a linear combination of the previous vectors. Consistent with Equations 3 and 4, we then have:

[TABLE]

From Equations 4 and 6, we have $AP_{1}\approx P_{2}\approx P_{1}H.$

In Schmid2 , the author describes a more robust solution, which is achieved by applying a singular value decomposition (SVD) on $P_{1}$. From Equation 6, SVD decomposition on $P_{1}$ subspace is calculated to obtain $U$, $\Sigma$ and $V^{*}$ matrices that are left singular vectors, singular values and right singular vectors respectively. The inversions of these matrices are then multiplied with $P_{2}$ subspace to obtain the The full-rank matrix $\widetilde{H}$, determined on the subspace spanned by the orthogonal basis vectors $U$ of $P_{1}$. described by:

[TABLE]

Here, $U^{*}\in\mathbb{C}$ and $V\in\mathbb{C}$ are the conjugate transpose of $U$ and $V^{*}$, respectively; and $\Sigma^{-1}\in\mathbb{C}^{N\times N}$ denote the inverse of the singular values $\Sigma$. After obtaining the $\widetilde{H}$ matrix, the eigenvalue analysis is performed to obtain $\omega$ eigenvectors and $\sigma$ a diagonal matrix containing the corresponding eigenvalues.

[TABLE]

It is known that the eigenvalues of $\widetilde{H}$ approximate some of the eigenvalues of the full system $A$. The associated eigenvectors of $H$ provide the coefficients for the linear combination that is necessary to express the dynamics within the image sequence basis. The dynamic modes $\Psi$ are thus calculated as follows:

[TABLE]

Since the eigenvalues associated with $\widetilde{H}$ are complex, it is not possible to establish the order of the dynamic modes $\Psi$ 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.

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 $1$ in DCE-MRI image (Figure 4(a)) covers both background region ($2$) and kidney region ($3$), similarly the synthetically generated images follow the same structure as shown in Figure 4(c).

A series of $100$ 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 $40\times 40$). 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 $0.05$ $mmol/kg$ of Gd-DTPA (Magnevist) contrast agent, on a $1.5T$ 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 $1.5s$ collected for $180s$. 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 $100$ images is given as an input to DMD algorithm which produces $99$ DMD modes (In theory, for a image sequence with $N$ images, we obtain $N-1$ 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.

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 $118$ 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.

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.

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.

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.