3D kidney segmentation from abdominal diffusion MRI using an appearance-guided deformable boundary.

A new technique for more accurate automatic segmentation of the kidney from its surrounding abdominal structures in diffusion-weighted magnetic resonance imaging (DW-MRI) is presented. This approach combines a new 3D probabilistic shape model of the kidney with a first-order appearance model and fourth-order spatial model of the diffusion-weighted signal intensity to guide the evolution of a 3D geometric deformable model. The probabilistic shape model was built from labeled training datasets to produce a spatially variant, independent random field of region labels. A Markov-Gibbs random field spatial model with up to fourth-order interactions was adequate to capture the inhomogeneity of renal tissues in the DW-MRI signal. A new analytical approach estimated the Gibbs potentials directly from the DW-MRI data to be segmented, in order that the segmentation procedure would be fully automatic. Finally, to better distinguish the kidney object from the surrounding tissues, marginal gray level distributions inside and outside of the deformable boundary were modeled with adaptive linear combinations of discrete Gaussians (first-order appearance model). The approach was tested on a cohort of 64 DW-MRI datasets with b-values ranging from 50 to 1000 s/mm2. The performance of the presented approach was evaluated using leave-one-subject-out cross validation and compared against three other well-known segmentation methods applied to the same DW-MRI data using the following evaluation metrics: 1) the Dice similarity coefficient (DSC); 2) the 95-percentile modified Hausdorff distance (MHD); and 3) the percentage kidney volume difference (PKVD). High performance of the new approach was confirmed by the high DSC (0.95±0.01), low MHD (3.9±0.76) mm, and low PKVD (9.5±2.2)% relative to manual segmentation by an MR expert (a board certified radiologist).

Introduction

The prevalence of chronic kidney disease (CKD) in the U.S. (2017) was about 15% of the adult population and 662,000 were of end-stage kidney disease (ESKD) [1]. Although transplantation is the definitive therapy for ESKD, approximately 17,000 transplants were annually performed during this time due to limited availability of transplantable kidneys [2]. Thus, all efforts should be employed to prolong the survival rate of these kidneys. However, acute rejection (AR) is one of the most serious barriers to transplanted kidneys. Therefore, clinicians are in a bad need of a fast, accurate, and reliable diagnostic tool to early detect AR for a higher chance of rescuing the transplanted kidney [3], [4]. An imaging modality, such as diffusion-weighted magnetic resonance imaging (DW-MRI) with the ability to provide both anatomical and functional information, expose patients to no radiation or contrast agent like other imaging modalities (e.g., computed topography (CT) and dynamic MRI), is critically needed to develop a computer-aided diagnostic (CAD) system for the early detection of AR post-transplantation [5]. DW-MRI allows mapping of water molecules’ diffusion process in biological tissues, in-vivo and non-invasively. These water molecule diffusion patterns are quantified by apparent diffusion coefficients (ADCs) and reveal details about tissue (e.g., kidney) status, either normal or diseased. In particular, healthy renal allografts provide higher ADC values than those with AR ones and thus facilitate the transplant status classification process [3], [4], [6], [7]. Segmentation of transplanted kidney from the surrounding abdominal structures and tissues is the key step to developing such a CAD system [8]. Thus, in literature several studies have focused their research on kidney segmentation. However, kidney segmentation from DW-MRI is very rare and most of it is performed by clinical research either manually or by an ROI tool [9], [10], [3], [4]. To the best of our knowledge, our group is the only group who started to automatically segment kidneys from DW-MRIs. Therefore, we will overview the related work devoted to kidney segmentation on other image modalities. Namely, we will review some of these studies that used different segmentation methods (e.g., threshold, region growing and graph cuts, and evolving deformable boundary) and different imaging modalities (e.g., CT and dynamic MRI).

With threshold methods, the kidney is segmented through studying the pixel intensity distribution in a certain region of interest (ROI). The ROI can be placed in a manual, semi-automatic or automatic manner. Giele et al. [11] presented a method for segmenting and registering the kidney from dynamic contrast-enhanced MRI (DCE-MRI). Their method involves manual insertion of a contour surrounding the kidney in one image with high contrast, then the kidney motion in other images is handled using a phase difference movement detection technique. However, it had a low accuracy and it compensated for the translation resulting from the motion without referring to the rotation, which was revisited in [12]. Mavromatis and Sequeira [13] presented a segmentation algorithm for medical tissues based on calculating texture directional maximums to segment a cancerous kidney from CT images. Priester et al. [14] segmented the renal transplant in MR images using two groups of images; acquired before and after the injection of the gadolinium-diethylenetriamine pentaacetic acid (DTPA) contrast agent. The average images of each group are subtracted and the output is thresholded resulting in a binary mask. Post-morpholigical processing is conducted to obtain the renal cross section. Giele [15] used a similar subtraction technique followed by thresholding for renal segmentation. Post-processing using morphological operations was used to close kidney contours. Separating the cortex from the medulla involves generating two ROIs (outer and inner) on the shape of onion rings using erosion. The outer ROI contains the cortex tissues and the inner ROI also contains the medulla tissues. They separated the calix manually. Koh et al. [16] used 3D H-maxima morphological transform for renal segmentation, where training data and prior information are being excluded by using edge information and rectangular masks.

Pohle and Toennies [17] segmented the kidney cortex using a region growing method. Their method is capable of self-learning of the homogeneity criterion based on the properties of the region under consideration. Boykov et al. [18] presented a globally optimal segmentation method for registered renal dynamic multidimensional data using temporal Markov based graph cuts. A vector having intensity values over time is used to characterize every voxel. Some seed points are initially inserted over the objects and the background, which are used to estimate a 2D histogram to be used with energy minimization. The properties of regions and boundaries are compromised to meet the imposed constraints. In spite of the good results of this technique, it requires a manual interaction from user. Rusinek et al. [19] used rigid registration with graph cut for renal segmentation in order to handle displacements. Farag et al. [20] presented a segmentation method for kidney that combined shape constraints with boundary properties and regions using graph cuts. Chevaillier et al. [21] segmented the internal kidney components in DCE-MRI images. They categorized the pixels based on the contrast evolution with time, making use of vector quantization algorithm. As in [18], their method requires interaction from the user. Freiman et al. [22] proposed an automatic, non-parametric graph min-cut model based technique for segmentation of kidney from CT scans. In their technique both shape and intensity information are integrated in the model. The latter is iteratively estimated using expectation maximization, which meanwhile performs segmentation through estimating the maximum a posteriori Markov random field using the graph min-cut. Li et al. [23] automatically segmented the kidneys using wavelet based clustering. Yang et al. [24] classified kidney tissues using fuzzy c-mean clustering.

Furthermore, kidney segmentation was explored using evolving deformable boundary techniques. Leventon et al. [25] combined the shape and deformable model by attracting the level-set function to the likely shapes from a training set specified by principal component analysis (PCA). Wang et al. [26] discussed using constrained optimization with deformable contours and applied it for kidney segmentation. Their system handled noise by using region information as constraints in addition to boundary information. Tsagaan et al. [27] presented a technique for 3D kidney segmentations from CT scan using a deformable model, which is described by a non-uniform rational B-spline (NURBS) surface and a priori shape. They used the principal curvature as a shape feature. Several research has been performed to segment kidneys using DCE-MRIs in humans and rats [28], [29]. For human studies, large scale movement for kidney registration and segmentation was roughly handled by Sun et al. in [30] using translational gradient based similarity registration. This was followed by the subtraction of an image with high contrast from a pre-contrast one, where the resulting difference was used to find the kidney contour by applying the level-set method. This contour was transferred across other frames to estimate the parameters of the registration. Concerning rats, kidney contours were found by Sun et al. [28] using a level-set variational method, integrating a model for subpixel movement and smoothness temporal constraints. The medulla and the cortex were segmented using the level-set method [31]. Using hybrid region- and edge-based models added an improvement to the segmentation techniques utilizing deformable models [32]. Moreover, multiple object segmentation was performed using multiphase level-sets [33], [34]. Based on both intensity and shape prior information, Abdelmunim et al. [35], [36] implemented a variational level-set segmentation framework to segment medical shapes (e.g., kidney). Spiegel et al. [37] presented a 3D kidney segmentation approach from CT scans using an active shape model in which non-rigid registration is used to find the correspondence between input training data points. A deformable model based parametric framework for segmenting kidneys was presented by Yuskel et al. [38]. The evolution of the contour had a shape prior constraint using signed distance maps, in addition to an intensity based distribution constraint calculated using linear combination of discrete Gaussians (LCDG) [39]. Campadelli et al. [40] proposed a framework for automatic segmentation of kidneys in CT images. Their framework is intensity based and depends on a multiplanar fast marching technique. The framework is generic and can be used for segmentation of various organs due to its robustness to intensity and shape variations. Gloger et al. [41] used Bayesian statistics in addition to shape prior to build a level-set kidney segmentation framework. Cuingnet et al. [42] presented an automatic method for detecting and segmenting kidneys in CT images. After kidneys localization, probabilistic kidney maps were calculated using random forests, utilizing both intensity and first/second order derivatives of each voxel and its neighbors. Finally, template deformation algorithm was performed on these maps to extract kidney surface.

Although the aforementioned segmentation techniques have relatively succeeded in segmenting kidneys with high contrast from CT and dynamic MRI, these techniques were not designed to handle challenges that exist in DW-MRIs such as low signal-to-noise ratio (SNR), low contrast, and diffused boundaries due to the high intensity similarities between the kidney and its background, especially at high b-values, which hinder the aforementioned methods from accurate segmentation of kidneys from diffusion MRIs.

To overcome these limitations, this paper presents a DW-MRI geometric deformable kidney segmentation framework, shown in Fig 1, that is robust to noise and low contrast. This is achieved by combining higher-order Markov-Gibbs random field (MGRF) model parameters and adaptive shape model, in addition to the first-order visual appearance model, into a joint MGRF model. Our framework includes preprocessing using bias correction [43] and histogram equalization, estimation of the joint MGRF model parameters, and extraction of the kidney volume from the surrounding tissues using the level-sets guided by the estimated joint MGRF model. This paper presents the following specific contributions. (i) Higher order MGRF appearance model (up to the 4-order), which takes into account the spatial dependencies between each voxel and its nearest neighbors in the DW-MR images. This better accounts for low SNR and low contrast, in addition to intra-kidney variabilities. (ii) An adaptive shape prior model of the expected kidney shape, which has the advantage over the fixed shape prior model not only to increase the robustness against noise and low contrast, but also to handle kidney motions caused by breathing and heart beating and to account for kidney variability due to inter-patient anatomical differences. This is achieved by adapting the shape probabilities of the kidney volume to be segmented based on its visual appearance. It is worth mentioning that the kidney-background visual appearance, shape prior, and statistical spatial dependencies are adaptable to kidney labels, which is an advantage of the presented framework.

Fig 1

The proposed framework for kidney segmentation.

Materials and methods

Participants

In total, 72 patients who underwent kidney transplantation provided verbal consents to participate in this study. All scans and biopsies were performed from July 2014 to February 2017, and all kidney transplants were done at Mansoura University (Egypt), with the donated kidneys having been obtained from live donors. However, eight patients were excluded due to a later on refusal of the study and/or contraindications for the MRI such as metallic prostheses, technical problems, artificial valves, or claustrophobia. The remaining 64 patients, (52 males and 12 females) and range in age from 12 to 54 years (the mean age of 27.2 ± 10.2 years), were included in the study. Patients were divided into two groups: non-rejection (NR) and acute rejection (AR). As a part of routine medical care after transplantation, all patients of both groups were assessed with serum creatinine laboratory values with a normal (basal) level of ≤ 1.3 mg ⋅ dl−1. The NR group (17 patients) included patients with healthy graft function. Most of the NR group patients only underwent the DW-MRI scans two weeks after transplantation. The AR group (47 patients) included patients with acute renal rejection, based on renal biopsy histology. All patients of the AR group underwent the DW-MRI scans two weeks after transplantation and just before the renal biopsy. All patients were asked to hold respiration (breath) during the study to reduce respiratory effect. Both the DW-MRI scans and biopsy were examined by a nephrologist and a board certified radiologist.

Imaging protocol

The DW-MR images were acquired before any biopsy procedure by using a 1.5T SIGNA Horizon scanner (General Electric Medical Systems, Milwaukee, WI, USA). Coronal DW-MR images have been obtained by using a body coil and a multi-shot spin-echo echo-planar sequence (TR/TE, 8000/61.2 ms; bandwidth, 142 kHz; 1.28 × 1.28 mm2 matrix; section thickness of 4 mm; intersection gap of 0 mm; FOV of 36 cm; 7 acquired signals; water signals acquired at different b-values of (b0, b50, b100, b200, b300, b400, b500, b600, b700, b800, b900, and b1000) s/mm2) using a single-direction from right to left. Approximately 50 sections have been obtained in 60–120 seconds to cover the whole kidney.

Detailed methods

Given an input (3D + b-value) DW-MRI, the presented segmentation technique in Fig 1 performs the following steps: (i) preprocessing of the DW-MRI kidney volume to be segmented; (ii) estimating of the joint Markov-Gibbs random field (MGRF) model parameters, namely, the adaptive shape model and the DW-MRI visual appearance features; and (iii) extracting the kidney volume from the surrounding tissues using the level-sets guided by the joint MGRF model estimated in the previous step.

Achieving an accurate kidney segmentation is a challenging task [44], [38] because of kidney motions due to breathing and heart beating; kidney shape changes due to inter-patient anatomical differences; low contrast between the kidney and other abdominal structures, especially at higher gradient strengths and duration, or b-values (Fig 2); low SNR and artifacts that complicate image alignment [45]; and geometric distortions due to long acquisition time [46]. To overcome these challenges, our segmentation relies on multiple image features to accurately delineate the kidney and thus facilitates analysis of transplant status. Details of our segmentation pipeline and its basic notations are outlined below.

Fig 2

Typical coronal cross-section DW-MRI samples showing (a) low contrast between the kidney and surrounding abdominal tissues at b0; (b) inter-patient anatomical differences at b0, (c) low signal-to-noise ratio (SNR), especially, at higher b-values (e.g., b1000); (d) image artifacts; and (e) geometric distortion/diffused boundaries.

: For describing processing steps, let p = (x, y, z) denote a voxel at 3D position with discrete Cartesian coordinates (x, y, z) and let R = {(x, y, z): 0 ≤ x ≤ X − 1; 0 ≤ y ≤ Y−1; 0 ≤ z ≤ Z − 1} be a finite 3D arithmetic lattice of unit voxels. The lattice has the size of XYZ and supports both grayscale images and their parametric or region (segmentation) maps. A grayscale image, g = {g: p ∈ R; g ∈ Q}, takes voxel-wise values from a finite set, Q = {0, 1, …, Q − 1}, of Q integer gray levels, i.e. g: R → Q. A region map, m = {m: p ∈ R; m ∈ L}, takes voxel-wise values from a binary set of region labels, L = {0, 1}, where 0 and 1 indicate the background and kidney, respectively, i.e. m: R → L.

: The overall segmentation pipeline starts with a preprocessing step that combines histogram equalization with non-parametric bias correction [43] in which the noise and inconsistencies due to low-frequency non-uniformity, or inhomogeneity of intensities, are partially suppressed.

A 3D geometric (level-set-based) deformable boundary, is employed for the DW-MRI kidney segmentation due to successful results in a wide range of applications, including medical imaging, (e.g., for segmenting brain, prostate, liver, kidney etc.) [47], [48], [41], [5]. Due to simplicity, flexibility, and ability to handle complex shapes and topological changes independently of surface parameterizations, these deformable boundaries are more popular than the alternative parametric ones. Points of an object-background boundary at each time instant t are specified implicitly as a zero-level set, B = {p: p ∈ R; Φ(p, t) = 0}, of arguments of a specific higher-dimensional function, Φ(p, t), on the lattice R. The function is often a signed distance map: where is the distance from the point p to the boundary B, and d(p, b) is the Euclidean distance between two lattice points p and b, as shown in Fig 3. The function Φ(p, t) evolves in discrete time t = nτ with a fixed step, τ > 0, as [49]: where n = 0, 1, 2, …, is the time index; ∇Φ(p, nτ) is the spatial gradient of Φ(p, nτ): |a| denotes the magnitude of the vector a, and F(p) is a speed function guiding the evolution of an initial boundary B0, defined at the starting instant t = 0, i.e., for n = 0.

Fig 3

3D zero-level set of a function Φ(p = [x, y, z], t).

Most of the conventional speed functions quantify visual appearance differences between the object and its background in terms of mean values and variances of image intensities, intensity edges, or gradient vector flow, and similar regional signal characteristics. Thus, their guidance may fail if images to be segmented are noisy and/or object-background contrast is low. To accurately segment the kidneys from noisy and low-contrast DW-MRI, our guiding function accounts for not only regional kidney-background appearance, but also for a kidney shape prior and high-order spatial relations in the goal region map. To provide voxel-wise guidance for the evolving boundary, all appearance and shape descriptors are combined into a joint Markov-Gibbs random field (MGRF) model of a DW-MR image, g, and its binary kidney-background region map, m. The model is specified by a joint probability distribution P(g, m) = P(g|m)P(m), where P(g|m) and P(m) denote a conditional probability distribution of images, given a map, and an unconditional distribution of region maps, respectively. The latter distribution is factored into two terms: P(m) = Psp(m)P(m), where Psp(m) denotes an appearance-based adaptive shape prior, and P(m) is a high-order Gibbs probability distribution with potentials V. The potentials evaluate strengths of not only the nearest-neighbor pairwise dependencies, but also of triple- and quadruple-dependencies, which specify a higher-order spatially homogeneous MGRF model of region maps. These components of the joint image-map model are outlined below.

First-order kidney/background appearance model

To accurately model DW-MRI appearance, we approximate the empirical marginal probability distribution of intensities with a linear combination of discrete Gaussians (LCDG) [39]. The LCDG with two positive dominant components (one each for the kidney and background) and multiple sign-alternate subordinate components allow for separating the mixed marginal of the DW-MRI voxel-wise intensities into the two distinct LCDGs, each associated with the kidney or background label. This LCDG model adapts the segmentation to changing appearance, such as non-linear intensity variations caused by patient weight and the data acquisition system, and it separates individual submodels of the kidney and background intensities more accurately than a conventional mixture of only positive Gaussians. This adaptation yields a better initial region map after the voxel-wise classification of only the image intensities with no account for the kidney shape.

Higher-order spatial interactions model

Compared to other imaging modalities, lower SNRs and frequent artifacts [45], together with geometric distortions due to long acquisition time [46] and larger inhomogeneities of internal structures, such as cortex and medulla, in the DW-MRI hinder the kidney segmentation. To better account for intra-kidney variabilities, spatial dependencies between each voxel and its nearest neighbors in the DW-MR images have been incorporated into our segmentation. Incorporated spatial relationships not only reduce noise impacts, but also reveal homogeneities and thus enhance the overall segmentation accuracy. Unlike the conventional pairwise-spatial homogeneity descriptors (e.g., in [50]), we use a 4th-order MGRF with analytically estimated potentials to describe those relationships. To find the potential estimates, an initial kidney map, m, is constructed by a simple Bayes classification using joint voxel-wise shape and intensity probabilities. Then, inter-label spatial dependencies in this map, m, are modelled by the 4th-order spatial MGRF with the nearest 26-neighborhood shown in Fig 4(a). This model adds triple and quadruple clique families to the more conventional 2nd-order Potts MGRF [50].

Fig 4

The nearest 26-neighborhood of a voxel for the 4th-order spatial model and examples of its 2nd-order cliques (upper raw), 3rd-order cliques (middle raw), and 4th-order cliques (lower raw).

Note that the central voxel is shown in yellow, while its neighbors are shown (i) in red for the same plane and (ii) in blue and purple for the adjacent planes.

Let C be a family of s-order cliques of the interaction graph with nodes in the lattice sites p ∈ R and edges connecting interdependent pairs of the sites. Let A clique families describe spatial geometry of interdependencies of region labels in the kidney maps, m. Then the model is specified by the Gibbs probability distribution: where is a collection of potential functions for the families C, ν is the clique size (ν ∈ {2, 3, 4}) for the family C; is a label configuration on the clique, (i.e., a pair, triple, or quadruple of binary numbers 0 and 1), and Z is the normalizing factor, called the partition function, over the entire population of the maps:

For equiprobable binary labels, m ∈ {0, 1}, the marginal co-occurrence probabilities over the 2nd-, 3rd-, and 4th-order cliques are , , and , respectively. Provided the cardinalities of the clique families are close to the lattice cardinality for all the families a = 1, …, A and only the equality (“eq”) or inequality (“ne”) of all the clique-wise labels are taken into account, the corresponding estimates of the 2nd-, 3rd-, and 4th-order potentials are as follows: where Fa:eq(m) denote relative empirical frequencies of the equal binary labels in the cliques of each family C over a given training map m; “eq” and “ne” denote two equal or non-equal labels, respectively, for a 2nd-order clique; “eq” denote i equal labels for a 3rd- and 4th-order clique; ; ; ; and This approximation is used for computing the higher-order spatial probabilities Pr(l) of each label; l ∈ L.

Adaptive shape prior model

In addition to the distinct visual appearances, the well-known geometric shapes of medical structures can enhance the segmentation accuracy. Relying on this fact, we use an adaptive model of the expected kidney shape to both handle kidney motions, (e.g., due to breathing and/or heart beating), and account for the kidney’s variability due to inter-patient anatomical differences. In addition, the kidney DW-MR images are very noisy, especially at high b-values.

To build the shape prior, all the b0-scans (after excluding the test subject) of kidneys formed a training database, and these images have been manually delineated by an MRI expert (a board certified radiologist) to get their binary kidney/background region maps. One of the images was chosen as a database reference. All other images were aligned to the reference by a non-rigid 3D registration [51] minimizing the sum of squared voxel-wise intensity differences between the two images. Then, the kidney/background labels of the co-aligned region maps were used to learn the shape prior. It is worth noting that to select the best reference subject, we performed a normalized cross-correlation between the test subject and every other subject in the shape prior training database. The subject with the maximum correlation was selected to be our reference. The choice of the reference has minimal effect on the shape prior as our presented shape prior depends on both the mapped spatial location in addition to signal appearance. In our presented shape prior, an adaptive search space around the mapped spatial location, which maps each voxel from the test subject to the database, is used in searching for voxels within a predefined tolerance range to the test voxel appearance. This means that any misalignment errors in the registration step will be overcome by the adaptive process. Moreover, the final segmentation takes into account the first-order appearance and the higher order spatial interaction that will overcome any errors from the shape prior segmentation. Fig 5 illustrates the co-alignment of the training DW-MRI. Adapting the shape prior to each input DW-MR image to be segmented is guided by the visual appearance of the latter image. The probabilistic shape prior is built as a spatially variant independent random field of region labels Psp(m) = ∏ Pr(m), where Pr(l) is the marginal empirical probability of the label l ∈ L in the voxel p; ∑ Pr(l) = 1.

Fig 5

3D co-alignment of training DW-MRI datasets (S1:S) to a single reference: The first and second rows present the overlapped 3D kidney volumes before and after the alignment, respectively.

Note that the reference subject appears in magenta, while the targets are shown cyan.

As the DW-MR images are challenged by the motion, which might lead to different kidney masks at different b-values, each b-scan is segmented separately. The shape prior is built by segmenting manually and co-aligning the baseline scans (at b0 s/mm2) of all subjects. Then, the shape prior is applied to each b-scan in combination with the other estimated probabilistic models for that scan namely, the 1-order LCDG model of current kidney appearance, in terms of voxel-wise intensities, and the 4-order MGRF model of spatial interactions. This joint probabilistic model provides a stochastic force that guides the evolution of a deformable boundary to segment the kidneys at that b-scan. So, the segmentation of all b-value scans uses the same shape prior, but own intensity and spatial interactions models.

Algorithm 1 summarizes estimating and updating the appearance-guided shape prior for each test DW-MR image to be segmented (the test images are first removed from the training set).

Algorithm 1 Creating/Updating the Shape Prior

1. Preprocess the training DW-MR images by bias correction and histogram equalization.

2. Construct the shape database by applying the co-alignment [51] to the preprocessed DW-MR images.

3. Preprocess the DW-MR image for a test subject and co-align with the shape database.

4. For each voxel, p ∈ R, in the test DW-MR image, gtest, calculate its prior shape probabilities as follows:

(a) Use the co-aligning deformation field to relate the voxel p of the test image to the shape database lattice.

(b) Initialize a 3D window of size N1 × N2 × N3, centered at the related voxel in the shape database lattice.

(c) Find within the window all the voxels with the corresponding intensity, gtest:, in all the training images.

(d) If necessary, increase the window size and repeat Steps 4b to 4d until a non-empty set of such corresponding training intensities is found.

(e) Estimate label probabilities based on relative occurrences of each label in all the training voxels found.

Appearance- and shape-guided deformable model

Adaptation to the kidney-background visual appearance, shape prior, and statistical spatial dependencies between kidney labels is one of the main advantages of our segmentation framework. Estimated directly from the input image and a given shape database, these properties guide the evolving deformable boundary by defining, for each voxel p with intensity g = q, the speed function [5] of Eq (2), F(p) = κϑ, where κ is the mean contour curvature and ϑ specifies the magnitude and direction of moving that voxel: Here, Pr(0) and Pr(1) are the voxel-wise background and kidney probabilities, respectively: where the variables Ωkd: and Ωbg: for the kidney and background, respectively, depend on the voxel-wise probabilities Pr(q|l); l ∈ L, for the LCDG submodels of the kidney (l = 1) or background (l = 0) appearance and on the kidney label probability in the MGRF spatial region map model, Pr(1), and in the adaptive shape prior, Prsp:(1), respectively: Algorithm 2 summarizes the basic steps of the 3D level-set-based kidney segmentation.

Algorithm 2 DW-MRI Segmentation by Geometric Deformable Boundary

1. Approximate the marginal of DW-MRI intensities with the LCDG [39] with two dominant components.

2. Update the shape prior probability using Step 4 of Algorithm 1.

3. Form an initial region map, mini, using the estimated shape prior and LCDG submodels of kidney and background appearances.

4. Estimate the Gibbs potentials for the 4-order spatial MGRF map model from mini.

5. Find the above speed function [5], F(p), using results of Steps 1 to 4.

6. Segment the input image, g, by evolving the level-set function, Φ(p, nτ), of Eq (2) with the speed function found in Step 5.

Results and discussion

Segmentation results: The performance of the proposed segmentation approach was tested on the collected DW-MRI data (a total of 64 subjects). Fig 6 shows segmentation results for different kidney cross-sections (coronal, axial, and sagittal) for two subjects at b0. The segmentation accuracy was evaluated by two volumetric measures, namely, the Dice similarity coefficient (DSC) [52] and percentage kidney volume difference (PKVD), and one distance-based metric—the 95-percentile modified Hausdorff distance (MHD) [53], which characterize the spatial overlap and distribution of the surface to surface distances between the segmented and ground truth kidneys, respectively. The ground truth kidney maps were manually delineated by an MRI expert (a board certified radiologist).

Fig 6

Our segmentation (red) with respect to the expert’s manual ground truth (green): The coronal (left), axial (middle), and sagittal (right) cross-sections for two different subjects in the first and second rows.

To show the effect of adding the higher-order MGRF model to our segmentation, we compared the current results with our previous segmentation using the 2 order- MGRF [54], [55], as shown in Fig 7. Moreover, Table 1 compares both segmentation methods in terms of the DSC, MHD, and PKVD metrics. As shown in Fig 7 and documented in Table 1, our segmentation has been notably enhanced after adding the higher-order MGRF model. This can be explained in part by abilities of the latter model to capture more intricate inhomogeneities of grey levels between different structures, (e.g. cortex and medulla.)

Fig 7

Our segmentation (red) with respect to the expert’s manual ground truth (green) using the 4-order MGRF (first row) compared to our previous segmentation using the 2-order MGRF (second row) [54], [55] for three different subjects (columns), where the first, second, and third columns show large, moderate, and small differences in yellow regions (false positive (FP)) and blue regions (false negative (FN)), respectively.

Table 1

Our segmentation accuracy by the DSC, MHD (mm), and PKVD (%).

All metrics are represented by the minimum (Min), maximum (Max), and mean±standard deviation (SD) values.

	Evaluation Metrics
Our method	DSC	MHD (mm)	PKVD (%)
Min	0.92	2.5	5.7
Max	0.97	6.4	15.7
Mean±SD	0.95±0.01	3.9±0.76	9.5±2.2
Previous [54], [55]	DSC	MHD (mm)	PKVD (%)
Min	0.89	4.0	9.4
Max	0.96	8.0	21
Mean±SD	0.92±0.02	5.7±2.0	14±3.1
P-value	<0.0001	<0.0001	<0.0001

Moreover, the advantages of the presented segmentation approach are highlighted by comparing our results with those obtained from other three well-known segmentation approaches, namely, the vector level-sets [35] (VLS), the segmentation using the random forest classifier [42] (RFC), and the level-set by Chan and Vese [31] (CV). Fig 8 visually assesses the better segmentation accuracy of the proposed segmentation method compared to others [31], [35], [42] for two different subjects. Table 2 compares our segmentation method’s results with the results obtained from the previously mentioned methods in terms of the DSC, MHD, and PKVD metrics. As documented in Table 2, our segmentation performs much better than other segmentation methods by providing higher DSC, lower MHD, and lower PKVD values. Furthermore, our approach is confirmed, by the paired t-tests between the DSC, MHD, and PKVD values for our approach and the compared segmentation approaches, to have a statistically more significant performance than these approaches and evidenced by the P-values less than 0.05.

Fig 8

Comparative cross-sectional segmentation results for our approach (a), the vector level sets [35] (b), the segmentation using the Random Forest classifier [42] (c), and the traditional CV [31] level set (d) for two independent subjects (rows).

The model segmentation is shown in red with respect to the manual ground truth (green) from an expert.

Table 2

Segmentation performance comparison between the presented approach against three other well-known segmentation methods (vector level-sets (VLS), random forest classifier (RFC), and Chan and Vese (CV)) using DSC, MHD (mm), and PKVD (%).

All metrics are represented by the mean± standard deviation (SD).

	Evaluation Metrics
	DSC	MHD (mm)	PVKD (%)
	Mean±SD	Mean±SD	Mean±SD	P-value
Our method	0.95±0.01	3.9±0.76	9.5±2.2	——
VLS [35]	0.90±1.8	5±0.85	10.8±4.7	<0.0001
RFC [42]	0.84±1.43	10.4±3.4	29.5±2.1	<0.0001
CV [31]	0.71±11	76±12	46±13	<0.0001

Fig 9 shows 3D segmentation results for three different subjects with their associated evaluation metrics. In particular, the developed segmentation technique proved its ability to precisely segment the kidney, especially at higher b values. Shown in Fig 10, favorable comparative results between the proposed segmentation approach and the previously mentioned state-of-the-art segmentation methods for one subjects at b0 and higher b values (b500 and b1000), which also confirm the high accuracy and robustness of the proposed segmentation method to low contrast and noisy images.

Fig 9

Our 3D segmentation (red) with respect to the expert’s manual ground truth (green) for three subjects with the associated accuracy scores.

Fig 10

Comparative coronal cross-sectional segmentation results for the proposed approach (a), the vector level sets [35] (b), the segmentation using the Random Forest classifier [42] (c), and the traditional CV [31] level set (d) for one subjects at b0 (first row) and higher b values (b500 (second raw) and b1000 (third raw)).

The model segmentation is shown in red with respect to the manual ground truth (green) from an expert.

Performance of the proposed segmentation technique

After testing the proposed segmentation approach on the available DW-MRI datasets (a total of 64 subjects) in a leave-one-subject-out scenario, our system showed a notable better segmentation performance than the aforementioned well-known segmentation methods confirmed by the higher DSC of 0.95±0.01, lower MHD (mm) of 3.9±0.76 and lower PKVD (%) of 9.5±2.2.

Currently, the entire processing pipeline is implemented in MATLAB using a Dell Precision WorkStation T7500 (4 Intel Xeon W5590 CPUs (2 Cores/CPU) @ 3.33 GHz with 48 GB RAM and 4 TB RAID hard drive). For an input DW-MRI dataset of size 256 × 256 × 50 voxels of each b-value scan, the whole segmentation pipeline takes around 1.83 min. For processing all b-value scans and obtaining the final segmentation, the entire processing takes about 20.13 min on average.

Clinical value of the contributions

Preliminary segmentation results outlined in this paper have shown that the proposed segmentation approach can segment kidneys from surrounding tissues with high reliability and accuracy. Taking into consideration that kidney segmentation is the first step in developing a fully automated CAD system for the early detection of acute renal transplant rejection, the more accurate the segmentation is, the more accurate estimating discriminatory features from the segmented kidneys will be. This will in turn be suitable for clinical applications to differentiate between non-rejection and acute rejection renal transplant.

Moreover, the proposed fully automated segmentation technique provides fast segmentation per b-value scan (in approximately 1.83 min with this MATLAB version) including all the pre-processing steps and the non-rigid registration, which is compatible with clinical routine. These abilities of the proposed segmentation technique encourages us to start developing a complete CAD system that will help clinicians to initiate timely interventions with appropriate treatments, which in turn can improve the delivery of healthcare in the USA and worldwide as a new, automated, fast, and relatively inexpensive diagnostic tool for early detection of acute renal transplant rejection.

Conclusions and future work

In summary, this paper presented a segmentation technique is very promising to provide the guidance for the level-set geometric deformable model to accurately segment kidneys from the other abdomen structures and surrounding tissues using (3D + b-value) DW-MRI data. Our segmentation has been notably enhanced in terms of robustness to noise and low contrast, especially at higher b-values, and accuracy after adding the higher-order MGRF model combined with the adaptive shape model and the first-order visual appearance model. This can be explained in part by abilities of the higher-order MGRF model to capture more intricate inhomogeneities of grey levels between different structures, (e.g. cortex and medulla.) Moreover, the segmentation results hold promise to be incorporated in a complete CAD system for the early detection of acute renal transplant rejection, thus; the appropriate physiological parameters will be extracted from the segmented kidneys (e.g., apparent diffusion coefficient (ADC)) to be used as our future discriminatory features.

In the future, we plan to test the proposed segmentation technique on more datasets collected from different locations with different scanning parameters to prove the reliability of the presented technique in dealing with diverse datasets. In addition, we are currently working to improve and optimize the efficiency of the segmentation algorithm to reduce the total running time as recommended by our medical collaborator, by converting our codes to C++ or using the GPU. It is worth mentioning that we plan to investigate the effect of adding a 3D B-spline non-rigid registration step –after segmenting the kidneys– on handling the possible potential motion and image distortion between the individual b-values caused by eddy current and magnetic susceptibility changes, which might lead to a more accurate estimation of the ADCs to be used as our future discriminatory features between the acute rejection and non-rejection renal transplants. Furthermore, the presented segmentation technique can be extended to be applied on segmenting different soft tissues (e.g., prostate, liver, and spleen).