Renal compartment segmentation in DCE-MRI images

Medical Image Analysis 32 (2016) 269–280

Contents lists available at ScienceDirect

Medical Image Analysis journal homepage: www.elsevier.com/locate/media

Renal compartment segmentation in DCE-MRI images Xin Yang a,∗, Hung Le Minh a, Kwang-Ting (Tim) Cheng b, Kyung Hyun Sung c, Wenyu Liu a a

Huazhong University of Science and Technology, Wuhan, 430074, China University of California, Santa Barbara, Santa Barbara, 93106, USA c University of California, Los Angeles, Los Angeles, 90095-7350, USA b

a r t i c l e

i n f o

Article history: Received 26 September 2015 Revised 16 April 2016 Accepted 13 May 2016 Available online 16 May 2016 Keywords: DCE-MRI Kidney segmentation Renal compartment segmentation Image registration PCA k-means clustering

a b s t r a c t Renal compartment segmentation from Dynamic Contrast-Enhanced MRI (DCE-MRI) images is an important task for functional kidney evaluation. Despite advancement in segmentation methods, most of them focus on segmenting an entire kidney on CT images, there still lacks effective and automatic solutions for accurate segmentation of internal renal structures (i.e. cortex, medulla and renal pelvis) from DCE-MRI images. In this paper, we introduce a method for renal compartment segmentation which can robustly achieve high segmentation accuracy for a wide range of DCE-MRI data, and meanwhile requires little manual operations and parameter settings. The proposed method consists of five main steps. First, we pre-process the image time series to reduce the motion artifacts caused by the movement of the patients during the scans and enhance the kidney regions. Second, the kidney is segmented as a whole based on the concept of Maximally Stable Temporal Volume (MSTV). The proposed MSTV detects anatomical structures that are homogeneous in the spatial domain and stable in terms of temporal dynamics. MSTVbased kidney segmentation is robust to noises and does not require a training phase. It can well adapt to kidney shape variations caused by renal dysfunction. Third, voxels in the segmented kidney are described by principal components (PCs) to remove temporal redundancy and noises. And then k-means clustering of PCs is applied to separate voxels into multiple clusters. Fourth, the clusters are automatically labeled as cortex, medulla and pelvis based on voxels’ geometric locations and intensity distribution. Finally, an iterative refinement method is introduced to further remove noises in each segmented compartment. Experiments on 14 real clinical kidney datasets and 12 synthetic dataset demonstrate that results produced by our method match very well with those segmented manually and the performance of our method is superior to the other five existing methods. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The kidney maintains whole-body homeostasis by filtering and excreting metabolic waste products, regulating acid-based balance, and moderating blood pressure and fluid volume (Koeppen and Stanton, 2012). Decreased kidney function, caused by disorders such as diabetes mellitus and hypertension, is a growing worldwide problem, e.g. up to 5% of the world’s population suffers from end-stage renal disease (ESRD), who may rely on dialysis or kidney transplantation for therapy (European Kidney Health Alliance (EKHA), n.d.). Moreover, recent study has shown that chronic kidney diseases could be an important risk factor leading to cardiovascular diseases (Menon et al., 2005). Therefore, there is a strong

∗

Corresponding author. E-mail address: [email protected], [email protected] (X. Yang).

http://dx.doi.org/10.1016/j.media.2016.05.006 1361-8415/© 2016 Elsevier B.V. All rights reserved.

demand for methods of precisely monitoring renal function, enhancing the assessment of disease progression, the prognosis, and follow-up therapy. Conventional methods for diagnosis of renal dysfunction are based on clinical chemistry measurements, which usually suffer from low sensitivity to pathologies (e.g. the change in creatinine level is detectable only after 60% or more functional loss) and could overestimate the actual glomerular filtration rate (GFR) by up to 20% (Myers et al., 2006). In addition, these measurements cannot localize disease locations. Recently, DCE-MRI has received considerable attention for accurate assessment of regional renal function (Michaely et al., 2005; Prasad, 2006), providing one-stop comprehensive morphological and functional information, without the utilization of ionizing radiation. In DCE-MRI, as a contrast agent is injected into the bloodstream, the passage of the contrast agent through the organ is reflected by signal intensity changes over time in the MRI images (as shown in Fig. 1(a)-(d)). The patterns of the regional time-intensity curves (Fig. 1(e)) give functional information, while the 3D MRI images provide anatomical

270

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

Fig. 1. (a) – (d) An exemplar slice of DCE-MRI series including both kidneys at four temporal points 1, 3, 7, and 23 respectively. For this patient, 1156 abdomen images are taken at 34 temporal points. Typically, the kidney enhancement shows the following 3 stages: (1) cortical enhancement (20-30 s after administration of contrast) that mostly reflects the contrast within renal vasculature (Fig. 1(b)); (2) medullary enhancement that usually occurs about a minute later and is dominated by the contrast in renal tubules (Fig. 1(c)); and (3) enhancement of the collecting system (i.e. pelvis) several minutes afterwards (Fig. 1(d)). By analyzing the enhancement of the renal tissues as a function of time, i.e. timeintensity curves (Fig. 1(e)), one can determine clinically important single-kidney parameters such as the renal blood flow (RBF), glomerular filtration rate (GFR), and cortical and medullary blood volumes.

information, which helps in distinguishing the diseases that affect different parts of the kidney. Accurate segmentation of renal compartments (i.e. renal cortex, medulla and pelvis) is essential for functional kidney evaluation. In today’s clinical routine, renal compartment segmentation is usually carried out manually by either a radiologist or a technician who empirically selects three 3D volumes at temporal points when three compartments are highlighted respectively, and then delineates a contour of each structure based on his/her knowledge of anatomy. Manual segmentation is time-consuming and fastidious, and meanwhile subject to inter- and intra-observer variability. Even though automatic 3D segmentation methods like thresholding, 3D active contours (Huang, 2009) and region growing based techniques (Lin et al., 2003; Pohle and Toennies, 2001; Yan and Wang, 2010) can be applied to alleviate those problems, it remains challenging to distinguish internal renal structures from a single 3D volume due to the following limitations of DEC-MRI images: 1) low spatial resolution, signal-to-noise ratio and partial volume effects due to fast and repeatedly scanning, 2) inhomogeneous intensity and contrast changes during perfusion in each compartment, especially for disordered kidneys (as shown in Figs. 1(b)-(d), bottom part of the right kidney has different intensity changes from the top part). Prior work: Several papers in the literature addressed the problem of renal compartment segmentation. Li et al. (2011) handled cortex segmentation as a multiple-surface extraction problem, which is solved using the optimal surface search method based on a graph construction scheme (Yin et al., 2010). This method is primarily designed for 3D CT images and evaluated on CT data, thus valuable temporal information embedded in the intensity time courses of DCE-MRI images is not considered (i.e. the temporal intensity evolution is different for each of the three kidney compartments, as shown in Fig. 1(e)). To address this problem, several methods have been proposed. (Sun et al., 2002) proposed an energy function that exploits both the spatial correlation among voxels and the intensity change of every voxel across the image sequences. (Boykov et al., 2001) introduced a temporal Markov model to describe the temporal intensity changes for each pixel. (Rusinek et al., 2007) proposed an automated graph-cuts method for kidney segmentation. In (Zöllner et al., 2009) and ((Chevaillier et al., 2008), the authors employed k-means clustering of temporal intensity evolution to segment the three internal renal structures. In the paper (Song et al., 2008), the authors proposed a 4D level set framework to combine information from spatial anatom-

ical structures and temporal dynamics for kidney segmentation in DCE-MRI images. However, those methods were evaluated only on normal kidneys; in practice they are sensitive to pathologies. Recently, the authors in (Khrichenko and Darge, 2010) presented a program, named CHOP-fMRU, for computer-aided renal compartment segmentation and functional analysis. CHOP-fMRU involves several manual tasks, e.g. manual delineation of a rough kidney contour for initialization, and the quality of the segmentation and analysis results depends heavily on the quality of the manual tasks. There have been numerous dedicated research efforts (Cuingnet et al., 2012; Goceri, 2011; Khalifa et al., 2011; Mory et al., 2012; Spiegel et al., 2009; Yuksel et al., 2007) for automatic segmentation of kidneys as a whole. Some of the most popular ones rely on the shape and appearance models learned from a set of previously segmented kidneys. For instance, in (Spiegel et al., 2009) the authors learned the kidney mean shape and principal modes of variation via an Active Shape Model to constrain the segmentation results. In (Yuksel et al., 2007), the authors modelled the kidney shape and intensity distribution using a signed distance map and Gaussian mixture models, respectively. Similarly, in (Goceri, 2011) the authors modelled the kidney’s intensity distribution using Gaussian mixture model (GMM) and the model parameters are learned based on a given MR image dataset. In (Khalifa et al., 2011) and (Mory et al., 2012), the authors integrated shape priors into a geometric deformable model to extract the kidney region. (Cuingnet et al., 2012) trained regression forests for roughly detecting a kidney’s location and the kidney surface is then extracted using the random forest (Breiman, 2001; Criminisi et al., 2011) and template deformation (Mory et al., 2012) algorithms. A coarse-to-fine strategy is utilized for accelerating the segmentation process. Model-based approaches have demonstrated promising results for segmenting an entire healthy kidney. However, the structure of renal compartments is much more challenging to model than a whole kidney due to its high complexity and variability in shape, in particular for those of disordered kidneys or kidneys with disturbances arising from surrounding organs (as shown in Figs. 1(b)-(d), shapes of the left and right functional compartments are quite different). As a result, it’s challenging for model-based methods to achieve satisfactory results for dysfunctional kidneys and for inner compartment segmentation. In this paper, we propose a method which can segment renal compartments for DCE-MRI images from both healthy subjects and patients with kidney problems. The proposed method consists of five main steps, as illustrated in Fig. 2. First, the input sequence

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

271

Fig. 2. Illustration of the proposed renal compartment segmentation framework.

of temporal volumes is aligned in the spatial domain to correct the motion caused by a patient’s movement during the image acquisition. The intensity of the aligned data is then normalized to enhance the kidney region. Second, a whole kidney is automatically segmented based on the detection of Maximally Stable Temporal Volume (MSTV). The proposed detection of MSTV exploiting both 3D spatial correlation among voxels and temporal dynamics for each voxel to provide a reliable segmentation robust to noises from surrounding tissues and kidney shape variations. Third, voxels in the segmented kidney are described by N principal components (PCs) where N is an empirical parameter. Our extensive experimental results indicate that for all cases the first 10 PCs capture most information needed for renal compartment segmentation. Thus discarding the rest of the components can effectively remove temporal redundancy and suppress noises with little loss of useful information. K-means clustering of the PCs is then leveraged to cluster voxels into multiple compartments. Fourth, the segmented inner compartments from k-means is labeled (i.e. recognized) as the cortex, the medulla and the pelvis respectively based on the feature of the time-intensity curve and geometric location of each segmented compartment. Finally, an iterative refinement method based on the modified Otsu thresholding (Otsu, 1975) is proposed to effectively and efficiently remove noises in each segmented compartment. The proposed segmentation method has been evaluated using 14 real clinical kidney data and 12 synthetic kidney data, and compared with the manual results as well as the results produced by five state-of-the-art methods. The results of the proposed method show excellent agreement with the manual results and superior performance to the comparable methods. 2. Description of the proposed 5-step method 2.1. Step 1 – Pre-processing Regional time-intensity curves are important metrics for kidney disease diagnosis (Lee et al., 2007; Michaely et al., 2006). A slight misalignment of 3D temporal volumes could lead to great variations in time-intensity curves, as significant intensity changes could be a mixture of misalignment-caused changes and contrast agent changes. Figs. 3(a) and (b) compare the average timeintensity curves before (denoted in red color) and after (denoted in green color) motion correction. Significant intensity changes can be observed, indicating great impact of misalignment on the final time-intensity curves. However, during the acquisition some mis-

alignment is inevitable due to the patient’s motion arising from breathing, heartbeat and bowel peristalsis. To overcome this problem, we apply motion correction via image registration before further processing. Several studies have been introduced in the literature for registration in DCE-MRI images. For instance, (Gerig et al., 1991) proposed to use Hough transform to register the edges in an image to the edges of a mask. In (Giele et al., 2002), the authors introduced a phase difference movement detection method to correct kidney displacements. However, these studies required building a mask manually by drawing the kidney contour on a 3D DCE-MRI volume. In this work, we employed an automatic non-rigid registration method based on mutual information (Yuksel, 2005; Yuksel et al., 2005). Specifically, we gradually deformed the target image according to a Free Form Model (FFD) described by B-splines basis functions until the mutual information between the target image and a reference image is maximized, indicating an optimal alignment between the two images. A multi-resolution strategy is applied to reduce the computational cost and to increase the robustness of the algorithm. An important issue for DCE-MRI image registration is the selection of reference and target images. In our method, we choose the volume in which the kidney region is mostly highlighted as the reference and the remaining volumes as targets. To demonstrate the effectiveness of our registration method, we manually delineate the entire kidney contour and pelvis contour at temporal point T = 3 and T = 19 respectively. We overlay the manual segmentation of the entire kidney on the volume at T = 25 (as shown in Fig. 3(c)) and observe that there is obvious misalignment between the manual segmentation and true kidney contour at different temporal points. Similar results can be observed when overlaying the pelvis manual segmentation at T = 19 on the volume at T = 31 (as shown in Fig. 3(e)). Fig. 3(d) and (f) show the superposition results after registration. Clearly, themanually segmentation can be well aligned with the true contour, indicating the effectiveness of our registration method. The goal is to obtain the concentration changes of the contrast agent for each renal compartment with respect to the time, based on which we can perform renal function assessment. However, the relationship between the contrast concentrations and the obtained MR signal intensities is in general non-linear. High contrast concentration could lead to susceptibility which may cause the signal intensity to decrease with the increase of the concentration. Moreover, MR signal intensity varies with the pulse sequence parameters, the pre-contrast relaxation times, and the blood flow velocity.

272

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

Fig. 3. (a) and (b) compare the average time-intensity curves for the blue region before (denoted in red curve) and after (denoted in green color) registration respectively. The comparison results demonstrate that there are large intensity changes before and after registration, indicating great negative impact of misalignment among different temporal volumes on the final time-intensity curves. (c)–(f) evaluates the effectiveness of our registration methods. (c) and (d): manual segmentation of an entire kidney obtained at temporal T = 3 (denoted in red color) overlaid on exemplar original slices at temporal point T = 25 before and after registration. (e) and (f): manual segmentation of pelvis obtained at temporal point T = 19 (denoted in red color) overlaid on original images at temporal point T = 31 before and after registration. Results show that before registration there are obvious misalignments between the manual segmentation and true image data, and those misalignments can be almost removed after our registration process.

Fig. 4. (a) – (d): average time-intensity curves of four kidneys. The baselines Spre-contrast for these four data are selected at temporal points T = 8, 4, 5, and 3, respectively (denoted by red arrows). (e)-(f): an exemplar slice at different temporal points from the original temporal volumes (e) and after normalization (f), respectively.

The simplest approach to solve this problem is to express the contrast concentration as the relative intensity enhancement via intensity normalization as follows:

E (t ) =

Spost-contrast (t ) − Spre-contrast Spre-contrast

(1)

where Spre-contrastand Spost-contrast(t)are the baseline (i.e. the kidney region has the minimum intensity value and afterwards the kidney starts to be highlighted robustly) and the contrastenhanced signal intensities, respectively. We observe the fact that intensity values of kidney tissues are gradually enhanced during perfusion while most non-kidney tissues remain unchanged with time, thus by calculating the relative intensity enhancement we can also highlight the kidney regions. Proper selection of a pre-contrast volume is a key step in the normalization process, as for different DCE-MRI data the baseline locates at different temporal points. For instance, Fig. 4(a)-(d) show

the average time-intensity curves of four different kidneys. The baselines Spre-contrast for these four data are selected at temporal points T = 8, 4, 5, and 3, respectively. To select Spre-contrast automatically for each input data, we detect the temporal point at which the average intensity achieves a local minimum value and after which the average intensities are higher and all above a predefined threshold Tbaseline . In our experiments, Tbaseline is set to one third of the mean value of the data’s time-intensity curve. Fig. 4(e) and (f) display an exemplar slice at different temporal points before and after normalization. Clearly, after normalization, the background noises are mostly removed and the inner compartments of kidney are highlighted. 2.2. Step 2 – Initial segmentation based MSTV As contrast agent perfuses into a kidney, intensity contrast between the kidney and its surrounding tissues are enhanced

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

273

Fig. 5. (a) An exemplar slice of DCE-MRI series taken at 6 temporal points 3, 6, 13, 19, 33 and 43. For this patient, 2448 abdomen images were taken over 72 temporal points. As time increases, the cortex, the medulla and the pelvis are highlighted sequentially. (b) Segmentation of the 3D volume at a single temporal point T = 19 (indicated by purple rectangle in (a)) using a sequence of thresholds ranked in an increasing order. Figures closer to the left are segmentations using smaller thresholds and those closer to the right are results using larger thresholds. (c) Segmentation of the slice via thresholding at 6 temporal points. Red, green and blue rectangles denote voxels at three kidney locations. Voxels connected using solid lines are temporally connected. We define two segmented volumes are temporally connected if they are temporally adjacent and the voxel overlap between them is greater than 80%. A sequence of temporally connected volumes is denoted as a temporal volume.

gradually (as shown in Fig. 5(a)). As a result, when binarizing a 3D volume via thresholding after contrast injection, the voxels in the segmented kidney possess three characteristics: 1) they are stable across a wide range of thresholds (as shown in Fig. 5(b)); 2) most of them are spatially connected (i.e. one voxel is among the 26 neighbors of another voxel in the segmented kidney); 3) a large portion of them appear in the temporally adjacent segmented volumes (i.e. having large voxel overlap in their temporal dynamics, as shown in Fig. 5(c)). In contrast, intensity enhancement of non-kidney tissues is rare and random; therefore, segmented nonkidney voxels are sensitive to thresholds, usually are not spatially connected and have small overlap in temporal domain. Based on these characteristics, we segment the kidney tissues as a whole by detecting a structure that remains most stable across both a wide range of consecutive threshold values and a wide range of consecutive temporal points. We define such a structure as Maximally Stable Temporal Volume (MSTV) and propose to detect MSTV via connected component trees. In the following, we first provide formal definitions of Temporal Volume and MSTV, and then detail the method of MSTV detection based on connected component trees. (1) Temporal Volume. We denote vt , where 1 ≤ t ≤ T, as a set of spatially connected voxels segmented from the original 3D volume at temporal point t by thresholding. T is the total number of temporal points in a DCE-MRI series. Two temporally consecutive voxel sets vt−1 and vt , where 1 ≤ t − 1 < t ≤ T , are considered temporally connected if the segmented voxels in vt−1 and vt have greater than λ (λ is 80% in this study) in common (i.e. voxel overlap between them is greater than λ). If any two temporally consecutive voxel sets in a sequence ϑ = {v1 , . . . , vt−1 , vt }, where t ≥ 2, are temporally connected, we denote the sequence ϑ as a temporal t  volume. The cardinality of ϑ is defined as | ϑ | = |vi | 1

(2) Maximally Stable Temporal Volume. If ϑ j−1 and ϑj are two temporal volumes obtained using threshold values jj−1 j j−1 j j−1 1 and j respectively, and v1 ⊆ v1 ,…, vt−1 ⊆ vt−1 , vt ⊆

vtj , then ϑ j−1 is a subset of ϑj , i.e. ϑ j−1 ⊆ ϑ j . Let {ϑ1 , . . . , ϑ j−1 , ϑ j . . .}m , where 1 ≤ m ≤ M, be a sequence of nested temporal volumes and M is the total number of

nested sequences detected in DCE-MRI data, the stability of the sequence is defined as:

S ( j ) = β j × |ϑ j |

(2)

Eq. (2) evaluates the stability of a segmentation, where β j is the number of threshold values for which the cardinality of ϑj is still sufficiently close to that of ϑ j−β j , i.e. |ϑ j \ ϑ j−β j |< τ where τ is a

predefined parameter (e.g. 5%). In other words, β j indicates across how many consecutive threshold values that the resulting segmentation remain stable. Similarly, |ϑj | indicates across how many temporally adjacent volumes the resulting segmentation remains stable. Our goal is to search j∗ which provides maximum Sm (j∗) for {ϑ1 , . . . , ϑ j−1 , ϑ j . . .}m . For all M sequences of nested temporal volumes, we identify m∗ which achieves a maximum Sm∗ (j∗) which corresponds to a structure with an optimized characteristics of both large cardinality of ϑj∗ and stability of ϑj∗ across a wide range of parameters. We define the structure which achieves Sm∗ (j∗) as a maximally stable temporal volume (MSTV). The detection of MSTV is based on the concept of connected component trees. For each temporal volume at time t, we construct a tree by binarizing it using a sequence of thresholds in an ascending order (as shown in Fig. 6). The tree root corresponds to the segmentation based on the minimum threshold value, which classify all voxels as foreground. As the threshold increases, some voxels are classified as background, segmenting the volume into several disconnected sub-regions. Accordingly, the tree root is split into several children nodes, each node corresponds to a subset of spatially connected voxels. Temporal volumes can be detected by searching tree nodes whose sizes remain similar across at least two trees (as indicated by blue circles in Fig. 6(b)). Stable temporal volumes can be detected by searching tree nodes whose sizes change little across at least two tree levels (as indicated by red curves in Fig. 6(b)). Once all stable temporal volumes are detected, we search for a MSTV which achieves the maximum Sm∗ (j∗) value. The detection process of MSTV is illustrated in Pseudo Code 1. First, for every volumetric data at time t, where 1 ≤ t ≤ T, a connected component tree Trt is constructed by binarizing the volume vt using all possible thresholds j’s where 1 ≤ j ≤ J. A node at level j of tree Trt consists of a group of spatially connected voxels. Second, for every threshold j (i.e. level j of the trees), we transverse trees Tr

274

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

Fig. 6. (a) An illustrative example of segmentation using threshold values ranked in increasing orders. (b) Two connected component trees at consecutive temporal points T = t and t + 1. Pseudo Code 1: MSTV-based Whole Kidney Segmentation from a DCE-MRI series Input: DCE-MRI data V = {V1 , . . . Vt−1 , . . . VT }, All possible thresholds J = {1,…j,…J}, j is ranked in an increasing order Parameters λ and τ Output: Segmentation of a whole kidney Procedure 1: Connected Component Trees Construction for Vt = V1 , . . . VT for j = 1,…J Binarize Vt using threshold j→ a list of connected voxels, {· · · (vtj )i . . .}; Insert {· · · (vtj )i . . .} → the jth level of Tree Trt ; End Insert Tree Trt → Connected Component Trees T r = {T r1 , . . . T rt } End Procedure 2: Maximum Stable Temporal Volumes Detection for j = 1,…J j , . . .}k (1≤ k ≤ K) Find temporally connected sequence ϑ jk ={. . . , vt−1 K is the total number of temporal volumes in level j End Transverse Tr from root to the level J to detect M sequences of nested temporal volume {ϑ1 , . . . , ϑ j−1 , ϑ j . . .}m (1≤ m ≤ M) Calculate Sm (j∗) for each {ϑ1 , . . . , ϑ j−1 , ϑ j . . .}m Search m∗ which achieves maximum Sm (j∗),and sequence is {ϑ j∗−β j∗ , . . . , ϑ j∗ }m∗ Procedure 3: MSTV-based Whole Kidney Segmentation Construct histogram H={h1 ,…,hN }for voxels in {ϑ j∗−β j∗ , . . . , ϑ j∗ }m∗ , N = |ϑj∗ | for ϑk = {ϑ j∗−β j∗ , . . . , ϑ j∗ }m∗ Update corresponding items of H for voxels appear in ϑk End Select voxels whose votes in H is greater than α → whole kidney segmentation

= {T r1 , . . . T rT } along the temporal axis to search for all temporal volumes. Then we transverse trees Tr vertically from roots to the last level J to detect MSTV. Third, we select voxels which appear in most temporal volumes in MSTV to form the initial segmentation of the whole kidney. A dilation operation is also applied to fill holes inside the kidney segmentation. Our MSTV can be considered as an extension of Maximally Stable Extremal Region (Matas et al., 2004), one of the best 2D interest point detectors in computer vision, from 2D to 4D. Meanwhile, we propose a MSTV detection method based on connected component trees and apply it for kidney segmentation. Our MSTV does not impose any shape constraints on the segmentation results, thus it can adapt to kidney shape variations caused by renal dysfunction. 2.3. Step 3 – PCA – kmeans clustering Time intensity curves associated with the cortex, the medulla and the pelvis voxels are distinguishable from each other (as shown in Fig. 1(e)). Therefore, they could be separated by unsu-

pervised clustering based on time intensity curves. However, the raw temporal data often has a high dimension, which increases the computation cost and induces numerical problems. Moreover, not all voxels of the same tissue are highlighted at the same moment, leading to misalignment of curves which belong to the same tissue and in turn resulting in mis-classification. In addition, most dimensions in the raw temporal data are redundant; the redundancy would dilute the useful information and disturb the true distribution. To address these problems, we employ Principal Component Analysis (PCA) to reduce the dimension of the temporal data. Voxel features in the temporal dimension are then described by a set of principal components (PCs). PCA transforms the raw data into a new coordinate space, such that the greatest variance lies along the first coordinate (denoted as the 1st PC), the second greatest variance along the second coordinate, and so on. By discarding the less significant components, PCA reduces the dimension of dynamic data and suppress noises. In addition, by linearly combining several temporal dimensions to form a new feature dimension,

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

275

Fig. 7. First 10 PCs of a dynamic DCE-MRI data.

misalignment can be avoided. Thorough analysis of our experimental results show that, for all cases we experimented with, greater than 99.4% of the total information is included in the first 10 PCs. As illustrated in Fig. 7, the 1st PC captures most global information of a kidney, and the 2nd to the 10th components encode detailed information of inner structures. For the later PCs, the variances tend to be increasingly affected by noise. Based on our experiment study, we choose 10 PCs for further analysis. Once voxels in the segmented kidney are described by the first 10 PCs, unsupervised clustering is applied to separate voxels into three groups: the cortex, the medulla and the pelvis. Among a number of suitable clustering methods, we choose k-means clustering in this study due to its simplicity, efficiency and effectiveness. 2.4. Step 4 – Renal compartment labeling In this step, we identify three clusters, from the k clusters generated from PCA-kmeans, and label them as cortex, medulla and pelvis. Several features can be used to distinguish voxels in different compartments and we employ the following two features for compartment recognition: (1) The average geometric location flocation : We calculate the average coordinates of the kidney’s contour (cx_avg , cy_avg , cz_avg ) and denote C = (cx_avg , cy_avg , cz_avg ) as the geometric center of the kidney. For each voxel we calculate its L2 distance to the geometric center and calculate the average L2 distance of each cluster i as its feature flocation_i . As cortex, medulla and pelvis voxels locate at the outermost, middle and innermost layer of a kidney, therefore, flocation_cortex > flocation_medulla > flocation_pelvis ; (2) The temporal point when the average time-intensity curve peaks fTIC : cortex voxels reach their maximum intensity values first, followed by medulla voxels and finally pelvis voxels. That is, fTIC_cortex < fTIC_medulla flocation_ ib > flocation_ ic and fTIC_ ia < fTIC_ ib
consistent conclusion or there are more than one consistent labeling results based on the two features, we consider those clusters are challenging to label due to too much noises, over-segmentation of a compartment into multiple small clusters, or lack of some compartment tissues. For such cases, we provide a user-friendly GUI with drop-down menus, slices in coronal plane, and average time-intensity curves for each cluster. Users can manually label a cluster as “cortex”, “medulla”, “pelvis” or “unknown” by viewing slices and the curve of this cluster. We merge clusters which have the same renal compartment manual label into a single one for further processing. In our experiments we observe that, for all healthy kidneys there always exist a unique combination which automatically and correctly labels all three compartments based on features flocation and fTIC . For disordered kidneys (i.e. hydronephrosis) and kidneys undergone operations, some compartments were not highlighted or might not even exist, thus their flocation and fTIC did not always comply with the rules as those of normal kidneys. Consequently, this step could yield inconsistent or more than one consistent labeling results and demand manual labeling.

2.5. Step 5 – Refinement We propose an iterative refinement method for removing noises and for recovering mis-classification due to ambiguous boundaries between clusters in the principal component space. The refinement method starts from the segmented cortex, to the medulla and then to the pelvis. First, we calculate an average time-intensity curve based on the candidate cortex voxels obtained in Step 4. Then we detect a temporal point tmax at which the curve reaches its maximum intensity value. We consider tmax and its neighboring temporal points tmax- 1 and tmax+ 1 as the candidate moments at which the cortex tissues are maximally highlighted. We chose the maximum intensity enhancement (MIE) for each candidate voxel according to Eq. (3):

MIEcortex (i ) = max{Stmax −1 (i ), Stmax (i ), Stmax +1 (i )} i ⊂ {candidate cortex voxels}

(3)

where St (i) is the intensity signal at voxel i temporal point t . At tmax- 1 , tmax and tmax+ 1 , the intensities of medulla and non-kidney tissues do not change too much from their pre-contrast intensities. Thus, at those moments, the MIE of voxels should be much smaller than the cortex voxels and hence can be excluded by thresholding. We chose two thresholds (i.e. Thigh and Tlow ) to classify candidate cortex voxels into three categories. For voxels whose MIE is greater than Thigh or smaller than Tlow , we classify them as confident cortex voxels and noise voxels respectively; for voxels whose MIE is smaller than Thigh while greater than Tlow , we consider them as ambiguous voxels and send them to the spatial filtering module for further verification. Specifically, for every ambiguous voxel, we examine its spatially adjacent voxels. If all of them are labeled as cortex voxels after the noise removal step described above, we relabel it from a non-cortex to a cortex voxel. The two thresholds are automatically selected via a modified Otsu method (Otsu, 1975). We exhaustively search for the two thresholds that minimize the

276

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

Fig. 8. Illustration of the iterative refinement pipeline.

intra-class variance as shown by Eq. (4):

thigh = argmin{ω0 (t )σ (t ) + (1 + λ )ω1 (t )σ (t )} t∈[1,L]

2 0

2 1

tlow = argmin{ω0 (t )σ02 (t ) + (1 − λ )ω1 (t )σ12 (t )} t∈[1,L]

ω0 (t ) =

t−1  i=0

p(i ), ω1 (t ) =

L  i=t

(4)

p( i )

where w0 (t) and w1 (t) are the probabilities of the two classes separated by a threshold t and σ02 (t) and σ12 (t) are variances of the two classes, λ is a weight adjusting the purity of two classes. For λ greater than 0, the probability and variance of the 2nd class plays a more importance role, thus t tends to maximize the purity of the 2nd class. In our experiment, we set λ as 0.2. Once we remove some noises, we iterate the above process to recalculate the average time-intensity curve so that a more accurate tmax can be selected for the following thresholding process. After a few iterations the refinement process converges (i.e. the cortex class does not change any more) and final cortex segmentation is obtained. Fig. 8 illustrates the iterative refinement pipeline. Similarly, this iterative refinement method is also applied to the segmented medulla and pelvis. 3. Experiments and results 3.1. Experimental setup We implemented our method using Matlab (Release 2013a, Natick MA) on a PC running on a single Intel core i7 CPU, Windows 7 OS. For the non-rigid registration step, we described the Free Form Model by B-spline basis functions, which are uniformly placed on a grid of control points. We utilized a 3level multi-resolution strategy to reduce the computational cost of the non-rigid registration and for each level we set the distance between every two control points to 16. The number of iterations for optimization is set to 50 for each hierarchical level. The parameters λ and τ in Pseudo Code 1 are set to 80% and 5% respectively. For the PCA-kmeans clustering step, we set the number of clusters k to 3 for all testing cases. We compare our method with five state-of-the methods: CHOPfMRU (Khrichenko and Darge, 2010), an dedicated software for renal compartment segmentation from DCE-MRI images, two Graph Cut (GC) based methods presented in (Rusinek et al., 2007) and (Boykov et al., 2001) respectively, Region Competition (Zhu and Yullie 1996), a popular active contour method for segmentation, and Grow Cut (Zhu et al. 2014). In addition, to examine the importance of each step in our framework we also evaluate the performance of our method with four different configurations: without registration (denoted as w/o registration), without MSTV-based kidney segmentation (denoted as w/o MSTV), without PCA dimension reduction (denoted as w/o PCA), and without refinement. We used the software provided by (Khrichenko and Darge, 2010) for CHOP-fMRU, the implementation based on the source code provided by authors of (Rusinek et al., 2007) and (Boykov et al., 2001)

for the the two graph cut methods, the implementation in ITKSNAP (Yushkevich et al., 2006) for Region Competition and the implementation in 3D slicer, an open-source platform for medical image analysis, for Grow Cut. Since Region Competition (RC) and Grow Cut are originally designed for 2D or 3D, but not 4D, data, we manually selected 3D volumes at those temporal points when the cortex, the parenchyma (i.e. the cortex and the medulla) and a whole kidney respectively seem maximally highlighted. Region Competition and Grow Cut are applied to each 3D volume to segment the cortex, the parenchyma and a whole kidney; the medulla and the pelvis were obtained by subtracting the cortex from the segmented parenchyma, and by subtracting the parenchyma from the whole kidney, respectively.

3.2. Datasets The imaging protocols used in this study were approved by the local institutional review board. Each enrolled subjected was screen for MRI risk factors and provided informed consent in accordance with institutional policy. This study consists of evaluation of 14 kidney cases: 9 normal cases, 2 cases with potential hydronephrosis problems, and 3 cases with pyeloplasty operations. The MRI data acquisition was performed using a 3.0T GE MR 750 system. To minimize the risk of gadolinium in patients with impaired kidney function, a low dose at 1/5 of Gadovist was used as the contrast agent with the injection rate of 0.3 mL/s, followed by 10 mL saline chaser at the same rate. Bellows respiratory triggering was implemented resulting in a temporal phase every two respiratory cycles. A 3D T1-weighted gradient echo sequence with a dual-echo bipolar readout was used for data acquisition, and we used an in-house variable density Cartesian undersampling scheme called DISCO to perform high spatiotemporal resolution dynamic MRU (Saranathan et al., 2012). A two-point Dixon reconstruction was used for robust fat–water separation. Imaging parameters were: flip angle is 15o , TR = 3.56 ms, matrix size = 256 × 256, FOV = 340 × 340 mm2 , the total number of slices is 34, and slice thickness is 4 mm. To evaluate the segmentation performance under different levels of noises and blurs, which are two common sources of segmentation errors, we generated a synthetic 4D renal dataset with 12 kidney cases as (Rusinek et al., 2007). Specifically, we first generated representative normal cortical, medullary and pelvic timeintensity curves by averaging intensity values from cortex, medulla and pelvis of normal subjects. Then, we generated a 4D mask of renal compartments by manual segmentation of a healthy kidney (as shown in Fig. 9(a)). We replaced values of 4D mask with values from the averaged time-intensity curves (as shown by Fig. 9(b)) and then we blurred the generated image with a discrete Gaussian kernel and added random and spatial uncorrelated noises to it (as shown in Figs. 9(c)-(e)). We set the width of 3D Gaussian kernel σ as 0 and 1.5 mm, the magnitude of the random noise n as 0% and 10% of the average intensity values of the precontrast cortical signal, resulting in 4 simulated data for normal cases. Similarly, we

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

277

Fig. 9. (a) 4D mask of renal compartments obtained by manual segmentation. White, black, gray pixels denote cortex, medulla and pelvis respectively. (b) – (e) Synthetic data obtained by replacing the 4D mask values by the averaged time-intensity curves with added Gaussian noise and blurs. Table 1 Comparison of different segmentation methods.

Disordered Kidneys

Kidneys with Op Healthy Kidneys

Compart-ment

Our method

CHOP-fMRU [2010]

GC [2007]

GC [2001]

RC[1996]

GrowCut[2014]

Cortex Medulla Pelvis Cortex Cortex Medulla Pelvis

0.95 0.98 0.87 0.97 0.99 0.98 0.96

0.71 0.62 0.63 0.69 0.70 0.68 0.72

0.87 0.85 0.86 0.89 0.96 0.96 0.92

0.82 0.83 0.82 0.85 0.91 0.92 0.89

0.60 0.52 0.59 0.47 0.50 0.51 0.60

0.65 0.54 0.59 0.61 0.64 0.69 0.62

simulate data for kidney with insufficiency and with operations. Thus, a total of 12 synthetic renal examinations were generated. 3.3. Evaluation metric We evaluated the segmentation accuracy using the Dice Similarity Coefficient (DSC), a widely used metric to evaluate segmentation algorithms for different medical image modalities. The DSC is defined as:

DSC (S, G ) =

2 × |S ∩ G| |S | + |G|

(5)

where S and G represent the sets of automatically segmented voxels and manually segmented voxels respectively; |·| denotes the set cardinality. The DSC ranges from 0, if S and G do not overlap at all, to 1, if S and G are identical. 3.4. Experimental results We first compare our method with the five state-of-the-art methods for 14 clinical data. Table 1 summarizes the average DSC for all five methods. Our method outperforms the other four comparison methods. The DSC achieved by our method is over 0.9 for most cases except for the pelvis (0.87) of disordered kidneys. This might be because part of the pelvis tissues was not highlighted due to the renal disease, causing disconnection among voxels in the spatial and temporal domains. Comparing to the CHOP-fMRU, the DSC achieve by our method is 24%∼36% greater than that of CHOPfMRU. Comparing to the two graph cut based methods (GC [2001] and GC [2007]) our method achieves slightly better segmentation accuracy. However, in terms of efficiency and convenience, our method significantly outperforms than those two methods. Specifically, our method requires little user interactions and takes only 26 s for segmentation (i.e. 15 s for MSTV-based whole kidney segmentation, 7 s for PCA-kmean clustering, 2 s for renal compartment recognition and 2 s for refinement). In comparison, the two graph cut based methods takes 21mins for segmentation and demand manual labeling of 15 regions as seeds for each compart-

ment (Rusinek et al. 2007). Comparing to the two 3D segmentation methods Region Competition and Grow Cut, the average DSC of our method is 35%∼50% higher than those of Region Competition and 28%∼44% greater than those of Grow Cut. We believe the poor performance of Region Competition and Grow Cut is mainly due to the highly varying contrast during perfusion. Although we manually selected the most relevant volumes in different phases of perfusion, it remains challenging to distinguish internal renal structures from a single image sequence. Fig. 10 shows the visualization of a health kidney in our clinical data based on different methods. Table 2 compares the DSC of our methods with different configurations. Motion correction, MSTV, PCA dimension reduction and refinement are essential and complementary to accurate segmentation results, i.e. up to 17%, 13%∼41%, 21%∼46%, and 3%∼17% improvements are achieved by using motion correction, the MSTVbased whole kidney segmentation, by PCA dimension reduction, and by refinement respectively. Fig. 11 shows the visualization of a healthy kidney in our clinical data based on manual segmentation and our method with different configurations. Table 3 shows the DSC of our method for the simulated data. Generally speaking, even for simulated data with serious noises and blurs (as shown in Fig. 9(e)), our method can achieve over 85% and 82% segmentation accuracy for normal and abnormal kidneys respectively. 3.5. Functional analysis Based on the renal compartment segmentation results, the following kidney function parameters (Hadjidekov et al., 2011) were calculated to aid the diagnosis of renal diseases: (1) Calyceal transit time (CTT): The amount of time for the contrast agent to pass from the renal cortex to the calyces (i.e. the collecting system). The CTT is most useful when the contralateral kidney is normal and we can classify CTT as symmetric, delayed or rapid on the hydronephrotic side. If the CTT is symmetric (i.e. its value is similar to that of a

278

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

Fig. 10. Segmentation results of a healthy kidney obtained by: (a) our method, (b) CHOP-fMRU [2010], (c) Graph Cut based method [2007], (d) Graph Cut + Markov Model based method [2001] (e) Region Competition [1996], and (f) Grow Cut [2014]. The red, yellow and blue labels indicate segmentation of the cortex, medulla and pelvis, respectively. Table 2 Comparison of our method with different configurations.

Disordered Kidneys

Kidneys with Op. Healthy Kidneys

Compartment

Our method

w/o Registration

w/o MSTV

w/o PCA

w/o Refinement

Cortex Medulla Pelvis Cortex Cortex Medulla Pelvis

0.95 0.98 0.87 0.97 0.99 0.98 0.96

0.86 0.96 0.70 0.93 0.98 0.97 0.96

0.74 0.85 0.58 0.56 0.79 0.83 0.81

0.54 0.57 0.56 0.56 0.60 0.52 0.75

0.92 0.95 0.75 0.84 0.91 0.88 0.79

Fig. 11. Segmentation results of a healthy kidney obtained by: (a) manual delineation, i.e. ground truth, (b) our method, (c) our method without registration, (d) our method without MSTV, (e) our method without PCA, and (f) our method without refinement. The red, yellow and blue labels indicate segmentation of the cortex, medulla and pelvis, respectively.

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280 Table 3 Segmentation results for various levels of Gaussian blurs and noises. Parameters

Cortex

Medulla

Pelvis

Disordered Kidneys

σ =0, n=0% σ =0, n=10% σ =1.5, n=0% σ =1.5, n=10%

0.95 0.91 0.86 0.83

0.97 0.92 0.88 0.83

0.93 0.89 0.87 0.82

Kidneys with Op.

σ =0, n=0% σ =0, n=10% σ =1.5, n=0% σ =1.5, n=10%

0.98 0.94 0.89 0.85

-

-

Healthy Kidneys

σ =0, n=0% σ =0, n=10% σ =1.5, n=0% σ =1.5, n=10%

0.98 0.92 0.89 0.85

0.97 0.92 0.88 0.83

0.93 0.89 0.87 0.82

279

on the two methods is about 3.0%; and the average difference between the renal parenchymal volume values calculated based on the two methods is about 5.2%), and the diagnostic outcomes based on these two segmentation methods are consistent for all testing cases. For instance, for Data 1, according to the manual segmentation, the RTTs of the left and right kidneys are 382.01 s and 391.80 s, respectively, based on the manual delineation, and the RTTs of the left and right kidneys are 391.80 s and 391.80 s respectively, based on our proposed segmentation. The RTT values are within the range of 245 s to 490s, thus both left and right kidneys are diagnosed as potentially obstructed. Similar conclusion can be made for Data 2 - 4, which are from patients with pyeloplasty operations on one kidney, and for Data 5 - 7, which are from patients with health kidneys. 4. Conclusion and future work

(2)

(3) (4) (5)

normal contralateral kidney), it indicates that the system remains compensated during the diuretic challenge. If the CTT is delayed, it suggests an acute decrease in GFR and an increased re-absorption of urine by the renal tubules in response to increased intra-pelvic pressure (i.e. the system is de-compensated). Renal transit time (RTT): The amount of time for the contrast agent to pass from the renal cortex to the ureter below the lower pole of the kidney (Jones et al., 2004). If the RTT is less than 245 seconds, the system is considered non-obstructive. If the RTT is greater than 490 seconds the system is probably obstructed. An RTT between 245 and 490 seconds is considered equivocal and should be managed conservatively with close follow-up to ensure a stable renal function. Time-to-peak (TTP): The amount of time to reach maximal parenchyma (i.e. cortex and medulla) enhancement. Whole renal volume: 3D volume of the renal parenchyma and pelvis. Renal parenchymal volume: 3D volume of the contrastenhanced renal parenchyma.

Table 4 compares the functional parameters computed based on the segmentation by our proposed method and manual delineation. Although the segmentation by two methods (our proposed method and manual delineation) does not perfectly match, the differences between the parameters produced by two methods are relatively small (e.g. the average difference between the CTT values calculated based on our method and manual delineation is about 4.2%; the average difference between the RTT values calculated based on the two methods is about 5.8%; the average difference between the whole renal volume values calculated based

We present a method for the renal compartment segmentation of DCE-MRI, which requires little manual delineation and parameter configurations. The proposed method has the following unique characteristics: 1) an automatic method for selection of a pre-contrast volume (i.e. baseline) based on which the intensity of the DCE-MRI data is normalized to effectively enhance the kidney region and remove background noises; 2) an automatic method for robustly and efficiently segmenting a kidney as a whole based on Maximally Stable Temporal Volume (MSTV) and connected component trees. MSTV-based kidney segmentation can well adapt to kidney shape variations due to renal dysfunction. In addition, it does not require any training and complex parameter configurations, which is very convenient for clinicians without image processing and computer vision background; 3) Using PCA and K-means clustering to automatically partition kidney voxels into multiple clusters; 4) a method based on the geometric locations and time-intensity curve of each cluster’s voxels for accurately identifying the compartment that a cluster belongs to; 5) an iterative refinement method based on maximum intensity enhancement for automatically and effectively removing noise and for recovering mis-detected kidney voxels. Experimental results on 14 clinical data demonstrate that: 1) The segmentation results of our method match well with those produced manually and our method achieves better performance than existing state-of-the-art segmentation methods; 2) Four essential steps in our method, namely motion correction, MSTVbased kidney segmentation, PCA dimension reduction and refinement, are complementary to each other, and their combination increases the robustness and accuracy of the segmentation result; 3) For all testing cases, the diagnosis results based on our method are fully consistent with those derived from manual delineation.

Table 4 Comparison of kidney parameters based on our automatic segmentation and manual segmentation. Data

Method

1

Our Manual

2

Our Manual

3

Our Manual

4

Our Manual

5

Our Manual

6

Our Manual

7

Our Manual

CTT (s) Left

Right

RTT (s) Left

Right

TTP (s) Left

Right

Renal volume (ml) Left Right

Parenchymal volume (ml) Left Right

362.42 372.21 312.80 322.56 176.44 186.24 225.29 244.87 195.90 205.70 135.98 135.98

362.42 372.21 176.13 185.92 156.72 166.52 254.67 264.47 210.15 222.52

391.80 382.01 234.60 224.83 166.63 176.44 166.52 176.31 195.90 186.11 148.34 176.31

391.80 391.80 195.70 215.27 156.72 166.52 195.90 205.67 197.79 210.15

117.54 195.90 175.95 205.26 137.23 127.43 190.21 190.21 166.52 166.52 86.53 97.95

264.47 284.06 137.00 146.78 140.15 140.15 127.36 205.67 111.29 98.90

37 38 25 25 102 103 73 73 116 116 27 28

32 32 21 20 97 97 71 70 104 106 22 22

39 40 48 48 75 74 18 22 45 46

31 32 42 41 72 72 12 19 40 41

280

X. Yang et al. / Medical Image Analysis 32 (2016) 269–280

Results on 12 synthetic data shows that even for MRU data with poor quality (i.e. large noises and blurs), our method can provide acceptable segmentation results for clinical usage. Our future work includes further acceleration of the current method, especially for the registration part, building more user friendly GUI for public release of the software and more evaluation on an even broader range of clinical data. References Boykov, Y., Lee, V.S., Rusinek, H., Bansal, R., 2001. Segmentation of dynamic N-D data sets via graph cuts using Markov Models. In: Medical Image Computing and Computer Assisted Interventions - MICCAI, pp. 1058–1066. Breiman, L., 2001. Random forests. J. Mach. Learn. 45 (1), 5–32. doi:10.1023/A: 1010933404324. Chevaillier, B., Ponvianne, Y., Collette, J.L., Mandry, D., Claudon, M., Pietquin, O., 2008. Functional semi-automated segmentation of renal DCE-MRI sequences. In: Acoustic, Speech and Signal Processing - ICASSP, pp. 525–528. doi:10.1109/ ICASSP.2008.4517662. Criminisi, A., Shotton, J., Konukoglu, E., 2011. Decision Forests for Classification, Regression, Density Estimation, Manifold Learning and Semi-Supervised Learning. Microsoft Research Cambridge, Tech. Rep. MSRTR. 5(6), pp. 81–227. doi:10.1561/ 060 0 0 0 0 035. Cuingnet, R., Prevost, R., Lesage, D., Cohen, L.D., Mory, B., Ardon, R., 2012. Automatic detection and segmentation of kidneys in 3D CT images using random forests. In: Medical Image Computing and Computer Assisted Interventions - MICCAI, pp. 66–74. European Kidney Health Alliance (EKHA), n.d. Factsheet: The kidney in Health and Disease. Gerig, G., Kikinis, R., Kuoni, W., von Schulthess, G.K., Kubler, O., 1991. Semiautomated ROI analysis in dynamic MRI studies. part i: image analysis tools for automatic correction of organ displacement. J. Comput. Assist. Tomogr. 15 (5), 725–732. Giele, E., Hasman, A., Engelshoven, J., Boer, J., 2002. Computer methods for semi-automatic MR renogram determination. Technische Universiteit Eindhoven. Goceri, E., 2011. Automatic kidney segmentation using Gaussian mixture model on MRI sequences. In: Electrical Power Systems and Computers, pp. 23–29. doi:10. 1007/978- 3- 642- 21747-0_4. Hadjidekov, G., Hadjidekova, S., Tonchev, Z., Bakalova, R., Aoki, I., 2011. Assessing renal function in children with hydronephrosis - additional feature of MR urography. J. Radiol. Oncol. 45 (4), 248–258. doi:10.2478/v10019-011-0038-z. Huang, C.-L., 2009. Shape-based level set method for image segmentation. In: International Conference on Hybrid Intelligent Systems - HIS, pp. 243–246. doi:10. 1109/HIS.2009.55. Jones, R.A., Perez-Brayfield, M.R., Kirsch, A.J., Grattan-Smith, J.D., 2004. Renal transit time with MR urography in children. J. Radiol. 233 (1), 41–50. doi:10.1148/ radiol.2331031117. Khalifa, F., Elnakib, A., Beache, G.M., Gimel’farb, G., El-Ghar, M.A., Ouseph, R., Sokhadze, G., Manning, S., McClure, P., El-Baz, A., 2011. 3D kidney segmentation from CT images using a level set approach guided by a novel stochastic speed function. In: Medical Image Computing and Computer Assisted Intervention -MICCAI, pp. 587–594. doi:10.1007/978- 3- 642- 23626- 6_72. Khrichenko, D., Darge, K., 2010. Functional analysis in MR urography - made simple. J. Pediatr. Radiol. 40 (2), 182–199. doi:10.10 07/s0 0247-0 09-1458-4. Koeppen, B.M., Stanton, B.A., 2012. Renal Physiology: Mosby Physiology Monograph Series. Elsevier - Health Sciences Division, USA. ISBN: 9780323086912. Lee, V.S., Rusinek, H., Bokacheva, L., Huang, A.J., Oesingmann, N., Chen, Q., Kaur, M., Prince, K., Song, T., Kramer, E.L., Leonard, E.F., 2007. Renal function measurements from MR renography and a simplified multicompartmental model. Am. J. Physiol. 292 (5), F1548–F1559. doi:10.1152/ajprenal.00347. Li, X., Chen, X., Yao, J., Zhang, X., Tian, J., 2011. Renal cortex segmentation using optimal surface search with novel graph construction. In: In: Medical Image Computing and Computer Assisted Intervention - MICCAI, pp. 387–394. doi:10.1007/978- 3- 642- 23626- 6_48. Lin, D.T., Lei, C.C., Hsiung, S., 2003. An Efficient Method for kidney segmentation on abdominal CT Images. In: Australian and New Zealand Intelligent Information Systems Conference, pp. 75–82. Matas, J., Chum, O., Urban, M., Pajdla, T., 2004. Robust wide-baseline stereo from maximally stable extremal regions. J. Image Vis. Comput. 22 (10), 761–767. doi:10.1016/j.imavis.20 04.02.0 06. Menon, V., Gul, A., Sarnak, M.J., 2005. Cardiovascular risk factors in chronic kidney disease. J. Kidney Int.. 68, 1413–1418. doi:10.1111/j.1523-1755.20 05.0 0551.x .

Michaely, H.J., Schoenberg, S.O., Oesingmann, N., Ittrich, C., Buhlig, C., Friedrich, D., Struwe, A., Rieger, J., Reininger, C., Samtleben, W., Weiss, M., Reiser, M.F., 2006. Renal artery stenosis: functional assessment with dynamic MR perfusion measurements–feasibility study. J. Radiol.. 238, 586–596. doi:10.1148/ radiol.2382041553. Michaely, H.J., Schoenberg, S.O., Rieger, J.R., Reiser, M.F., 2005. MR angiography in patients with renal disease. J. Magn. Reson. Imaging Clinics North Am. 13 (1), 131–151. doi:10.1016/j.mric.20 04.12.0 07. Mory, B., Somphone, O., Prevost, R., Ardon, R., 2012. Real-time 3D image segmentation by user-constrained template deformation. In: Medical Image Computing and Computer Assisted Intervention - MICCAI, 15, pp. 561–568. Myers, G.L., Miller, W.G., Coresh, J., Fleming, J., Greenberg, N., Greene, T., Hostetter, T., Levey, A.S., Panteghini, M., Welch, M., Eckfeldt, J.H., 2006. Recommendations for improving serum creatinine measurement: a report from the laboratory working group of the national kidney disease education program. J. Clinical Chem. 52 (1), 5–18. doi:10.1373/clinchem.2005.0525144. Otsu, N., 1975. A threshold selection method from gray-level histograms. J. Autom.. 11, 285–296. Pohle, R., Toennies, K.D., 2001. Segmentation of medical images using adaptive region growing. in medical imaging. J. Int. Soc. Optics Photon. 1337–1346. doi:10. 1117/12.431013. Prasad, P.V., 2006. Functional MRI of the kidney: tools for translational studies of pathophysiology of renal disease. Am. J. Physiol. 290 (5), F958–F974. doi:10. 1152/ajprenal.00114.2005. Rusinek, H., Boykov, Y., Kaur, M., Wong, S., Bokacheva, L., Sajous, J.B., Huang, A.J., Heller, S., Lee, V.S., 2007. Performance of an automated segmentation algorithm for 3D MR renography. J. Magn. Reson. Med. 57 (6), 1159–1167. doi:10.1002/ mrm.21240. Saranathan, M., Rettmann, D.W., Hargreaves, B.A., Clarke, S.E., Vasanawala, S.S., 2012. Differential subsampling with cartesian ordering (DISCO): A high spatio temporal resolution dixon imaging sequence for multiphasic contrast enhanced abdominal imaging. Magn. Res. Imag. 35 (6), 1484–1492. Song, T., Lee, V.S., Rusinek, H., Chen, Q., Bokacheva, L., Laine, A., 2008. Segmentation of 4D MR renography images using temporal dynamics in a level set framework. In: Biomedical Imaging: From Nano to Macro, 20 08. ISBI 20 08. 5th IEEE International Symposium, pp. 37–40. doi:10.1109/ISBI.2008.4540926. Spiegel, M., Hahn, D.A., Daum, V., Wasza, J., Hornegger, J., 2009. Segmentation of kidneys using a new active shape model generation technique based on nonrigid image registration. J. Comput. Med. Imaging Graph. 33 (1), 29–39. doi:10. 1016/j.compmedimag.2008.10.002. Sun, Y.S.Y., Moura, J.M.F., Yang, D.Y.D., Ye, Q.Y.Q., Ho, C.H.C., 2002. Kidney segmentation in MRI sequences using temporal dynamics. In: Biomedical Imaging. IEEE International Symposium, pp. 98–101. doi:10.1109/ISBI.2002.1029202. Yan, G., Wang, B., 2010. An automatic kidney segmentation from abdominal CT images. In: Intelligent Computing and Intelligent Systems - ICIS, pp. 280–284. doi:10.1109/ICICISYS.2010.5658676. Yin, Y., Zhang, X., Williams, R., Wu, X., Anderson, D.D., Sonka, M., 2010. LOGISMOSlayered optimal graph image segmentation of multiple objects and surfaces: cartilage segmentation in the knee joint. J. Med. Imag. IEEE Trans. 29 (12), 2023–2037. doi:10.1109/TMI.2010.2058861. Yuksel, S.E., 2005. Image processing methods for the detection of acute rejection after kidney transplantation. In: MS thesis, Department of Electrical and Computer Engineering. University of Louisville, Louisville, KY, USA, pp. 27–46. Yuksel, S.E., El-Baz, A., Farag, A.A., Abo El-Ghar, M.E., Eldiasty, T.A., Ghoneim, M.A., 2005. Automatic detection of renal rejection after kidney transplantation. In: International Congress Series, 1281, pp. 773–778. doi:10.1016/j.ics.2005.03.146. Yuksel, S.E., El-Baz, A., Farag, A.A., El-Ghar, M., Eldiasty, T., Ghoneim, M.A., 2007. A kidney segmentation framework for dynamic contrast enhanced magnetic resonance imaging. J. Vib. Control 13, 1505–1516. doi:10.1177/1077546307077417. Yushkevich, P.A., Piven, J., Hazlett, H.C., Smith, R.G., Ho, S., Gee, J.C., Gerig, G., 2006. User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability. J. Neuroimage 31 (3), 1116–1128. doi:10.1016/j.neuroimage.2006.01.015. Zhu, L., Kolesov, I., Gao, Y., Kikinis, R., Tannenbaum, A., 2014. An effective interactive medical image segmentation method using fast growcut. In: Interactive Medical Image Computing Workshop, pp. 11–20. Zhu, S.C., Yullie, A., 1996. Region competition: unifying snakes, region growing, and bayes/mdl for multiband image segmentation. J. Pattern Anal. Mach. Intell. 18, 884–900. doi:10.1109/34.537343. Zöllner, F.G., Sance, R., Rogelj, P., Ledesma-Carbayo, M.J., Rørvik, J., Santos, A., Lundervold, A., 2009. Assessment of 3D DCE-MRI of the kidneys using non-rigid image registration and segmentation of voxel time courses. J. Comput. Med. Imag. Graph. 33 (3), 171–181. doi:10.1016/j.compmedimag.20 08.11.0 04.