Keywords

1 Introduction

Urinary bladder cancer (BC) is the fourth most common cancer and ranks ninth as a cause of cancer related deaths among males globally [1, 2]. Accurate identification of tumor involvement, especially muscle-invasive depth (MID), is of paramount significance for treatment strategy [1, 3,4,5]. In recent years, optical cystoscopy (OCy) with transurethral biopsies is the standard reference for monitoring BC, muscle invasiveness prediction, MID determination as well as staging diagnosis [1,2,3,4,5]. However, it is an invasive manner with high expenses and large time consumption, bringing great pain to the patients and ascending the risk of urinary tract bleeding and infection [1,2,3,4,5]. Non-invasive imaging-based biomarkers may provide a direct evaluation of the disease. For example, the local-part thickening of bladder wall in medical images has become a useful indicator for non-surgical detection and extraction of suspected BC regions [1, 6,7,8,9].

Excellent resolution of soft tissues, non-invasive procedures and easy performance of three-dimensional (3D) magnetic resonance imaging (MRI) make it ideally suitable for in vivo bladder imaging, as well as observing and detecting BC with non-surgical imaging-based biomarkers [1, 2, 6,7,8,9]. In order to quantitatively compute such biomarkers from bladder images, it is required to accurately delineate the outer-wall boundary (OB) and inner-wall boundary (IB) of bladder to separate the whole bladder wall from urine and background regions [2, 6,7,8,9]. At present, the segmentation of bladder wall in MRI is manually performed by experts [6, 9]. However, it is a tedious and time-consuming process. An accurate and efficient segmentation algorithm, other than the manual segmentation, is in desirable need to extract the whole bladder wall region from bladder MR images.

The main objective of this paper is to explore an accurate 3D multi-region segmentation method for separating bladder wall from urine and outer background regions on T2-weighted MR images (T2WI) efficiently and robustly, which may greatly help the computer assisted detection and diagnosis of bladder cancer in clinical researches and applications.

Imaging-based segmentation of bladder wall is a difficult work full of challenges, including bladder shape variation, artifacts in the urine region, weak boundary of bladder wall and complicated background intensity distributions [6, 9]. So far, very a few bladder wall segmentation methods have been proposed. Duan et al. developed a coupled level-set framework, which adopted a modified Chan-Vese model to locate the inner and outer surfaces of bladder wall, and then generated a coupled model for segmenting both surfaces of bladder wall on T1-weighted MR images (T1WI) [6, 9]. However, it is difficult to locate the outer-wall boundary (OB) due to the complicated outside intensity distribution. Qin et al. proposed an adaptive shape prior (ASP) constrained coupled directional level-set algorithm, to fully exploit the shape prior information like directional gradient, inter-slice and regional shape variation, minimum thickness and so on, to obtain a robust segmentation result [2, 6,7,8,9]. Cha et al. incorporate the deep learning convolutional neural network (DL-CNN) with conjoint level-set analysis and segmentation system (CLASS) to realize the bladder segmentation on CT images [10]. These methods are either easily influenced by the initialized contours of OB and inner-wall boundary (IB) manually delineated on the middle bladder slice, or limited by time-consuming and complicated procedures [11]. Moreover, their core mathematical schemes are based on local optimization techniques, which are susceptible to image quality [11].

Apart from the level-set based algorithms mentioned above, Ukwatta and Yuan proposed a novel 3D multi-region segmentation algorithm based on coupled continuous max-flow (CMF) and convex relaxation optimization methods to simultaneously segment the wall region of carotid from the outer background region and inner lumen region [11]. Distinct from the discrete global optimization methods like graph-cuts, the convex relaxation methods are employed in a spatially continuous setting, providing sub-pixel accuracy while avoiding metrication artifacts [11,12,13,14,15,16].

Inspired by their work, in this study we try to develop an efficient convex optimization based 3D multi-region segmentation algorithm for bladder T2-weighted MR images (T2WI), basing on the continuous max-flow (CMF) model with anatomical inter-surface constraint and bladder wall thickness (BWT) prior for efficient convex relaxation optimization.

2 Methods

In this section, an accurate multi-region continuous max-flow (MR-CMF) segmentation method is proposed, which separates a dataset of 3D bladder image \( V \) into three regions: the urine region \( R_{u} \), the bladder wall region \( R_{w} \) and the outer background region \( R_{b} \), by simultaneously evolving two coupled surfaces including the inner-wall boundary (IB) \( S_{IB} \) and the outer-wall boundary (OB) \( S_{OB} \) to their globally optimized positions, as illustrated in Fig. 1.

Fig. 1.
figure 1

Major regions and surfaces to be delineated in bladder image \( V \)

Full size image

These three regions should strictly meet such conditions:

$$ \left\{ {\begin{array}{*{20}l} {R_{b} \cup R_{w} \cup R_{u} = V,} \hfill \\ {R_{b} \cap R_{w} = \varnothing } \hfill \\ {R_{w} \cap R_{u} = \varnothing} \hfill \\ {R_{b} \cap R_{u} = \varnothing} \hfill \\ \end{array} } \right. $$
(1)

Considering the prior knowledge of the anatomic information, the inner-wall surface \( S_{IB} \) always remains incapsuled by the outer-wall surface \( S_{OB} \), namely the anatomic inter-surface constraint:

$$ S_{IB} \subset S_{OB} $$
(2)

In the process of segmentation, when the two surfaces \( S_{IB} \) and \( S_{OB} \) are evolved in each iteration, the generated three regions \( R_{u} \), \( R_{w} \) and \( R_{b} \) are determined as:

$$ R_{u} : = \varOmega_{{S_{IB} }} ,\quad R_{w} : = \varOmega_{{S_{OB} }} \backslash \varOmega_{{S_{IB} }} ,\quad R_{b} : = V\backslash \varOmega_{{S_{OB} }} $$
(3)

where \( \varOmega_{{S_{IB} }} \) and \( \varOmega_{{S_{OB} }} \) are the regions enclosed by \( S_{IB} \) and \( S_{OB} \), respectively.

Then, we introduce a convex optimization based approach with a novel and efficient continuous max-flow (CMF) based algorithm to accurately segment the whole bladder wall region \( R_{w} \) from the input 3D bladder MR images. Distinct from the previous bladder segmentation algorithms, such as level-set that bases on the techniques of local optimization for contour evolution, in this method, the two surfaces \( S_{IB} \) and \( S_{OB} \) can be evolved to their globally optimal locations, while enforcing the anatomic inter-surface constraint (2). Moreover, the intensity probability density functions (PDF) of the three regions \( R_{u} \), \( R_{w} \) and \( R_{b} \) are exploited as global descriptors to guide the evolvement of the coupled surfaces \( S_{IB} \) and \( S_{OB} \) [11, 13,14,15,16].

2.1 Optimization Model and Convex Relaxation

Let \( L_{r} : = \{ u,w,b\} \) be the label set representing the three independent regions of image \( V \), and \( L_{s} : = \{ IB,OB\} \) represent the label set of two surfaces to be calculated in this work. Here, we adopted PDF of intensity distributions in the region \( R_{k} \) to describe the probabilities of allocating each element \( x \in V \) to the optimal region.

The PDF of intensities like intensity histogram is a kind of global descriptors of the objects of interest, such that matching the intensity PDF models of the object regions provides a robust mechanism to guide the evolvement of the surfaces [11, 13,14,15,16]. Let \( I(x) \in Z \) be a given 3D bladder MR image, where Z is the set of image intensity values. According to the Parzen method, the intensity PDF \( h_{k} (z) \), where \( z \in Z,k \in L_{r} \), for the respective three regions can be computed as [11, 13, 16]

$$ h_{k} (z) = \frac{{\int_{V} {K(z - I(x))u_{k} (x)dx} }}{{\int_{V} {u_{k} (x)dx} }},k \in L_{r} $$
(4)

where \( K( \cdot ) \) is the Gaussian kernel function described in Eq. (5) and \( u_{k} (x) \) is the labeling function of the corresponding region \( R_{k} \) which is enclosed by the surface \( S_{m} \), \( m \in L_{s} \), respectively.

$$ K(x) = \frac{1}{{\sqrt {2\pi \sigma^{2} } }}\exp ( - x^{2} /2\sigma^{2} ) $$
(5)
$$ u_{k} (x) : = \left\{ {\begin{array}{*{20}l} {1,x \in R_{k} } \hfill \\ {0,otherwise} \hfill \\ \end{array} } \right.,k \in L_{r} $$
(6)

In computation, we define \( u_{1} (x) \) and \( u_{2} (x) \) to be the binary region labeling functions of the regions enclose by \( S_{IB} \) and \( S_{OB} \), respectively.

$$ u_{1} (x) : = \left\{ {\begin{array}{*{20}l} {1,\quad x \in\Omega _{{S_{IB} }} } \hfill \\ {0,\quad otherwise} \hfill \\ \end{array} } \right. $$
(7)
$$ u_{2} (x) : = \left\{ {\begin{array}{*{20}l} {1,\quad x \in \varOmega_{{S_{OB} }} } \hfill \\ {0,\quad otherwise} \hfill \\ \end{array} } \right. $$
(8)

Considering the anatomic inter-surface constraint (2), the binary region labeling functions must satisfy the following equations:

$$ u_{1} (x) \le u_{2} (x),\quad \forall x \in V $$
(9)
$$ \left\{ {\begin{array}{*{20}l} {u_{u} (x) = u_{1} (x),\begin{array}{*{20}c} {} \\ \end{array} } \hfill \\ {u_{w} (x) = u_{2} (x) - u_{1} (x),} \hfill \\ {u_{b} (x) = 1 - u_{2} (x)} \hfill \\ \end{array} } \right. $$
(10)

Therefore, we define the cost function \( C_{k} (x) \), \( k \in L_{r} \), of labeling each element \( x \in V \) to be in the region \( R_{k} \), \( k \in L_{r} \) by log-likelihoods of the PDF value of \( x \) [13, 15, 16], i.e.

$$ C_{k} (x) = - \log (h_{k} (z)),\quad k \in L_{r} $$
(11)

With the combination of the global optimization of the intensity PDF models \( h_{k} (z) \) and their corresponding cost functions \( C_{k} (x) \) with geometric constraint (9)–(10) on the labeling functions, we formulate the following energy function to segment the given image \( V \) by achieving both the minimum total labeling costs and the minimum total regions/volumes of segmented surfaces [15, 16], such that

$$ \begin{array}{*{20}c} {\mathop {\hbox{min} }\limits_{{u_{1,2} (x) \in \{ 0,1\} }} } & {\sum\limits_{{k \in L_{r} }} {\int_{{{R_{k} }}} {C_{k} (x)dx} + \sum\limits_{{k \in L_{s} }} {\int_{{{\partial R_{k} }}} {ds} } } } \\ {\quad \quad \quad s.t.} & {\quad u_{1} (x) \le u_{2} (x),\quad \forall x \in V} \\ \end{array} $$
(12)

Considering the minimum geodesic length should be used to constrain the segmentation results, the proposed optimization model (12) can be equally written as follows:

$$ \begin{array}{*{20}l} {\mathop {\hbox{min} }\limits_{{u_{1,2} (x) \in \{ 0,1\} }} } \hfill & {\sum\limits_{{k \in L_{r} }} { < u_{k} ,C_{k} > } + \sum\limits_{{k \in L_{s} }} {\int_{\varOmega } {g(x)\left| {\nabla u_{k} (x)} \right|dx} } } \hfill \\ {\quad \quad \,\,\,\,s.t.} \hfill & {\quad u_{1} (x) \le u_{2} (x),\quad \forall x \in V} \hfill \\ \end{array} $$
(13)

where the weighted function \( g(x) \ge 0 \) and can be computed as follows [11, 13, 15, 16]

$$ g(x) = \lambda_{1} + \lambda_{2} \exp ( - \lambda_{3} |\nabla \mathop I\limits^{\sim} (x)|),\quad \lambda_{{1,\begin{array}{*{20}c} {} \\ \end{array} 2,\begin{array}{*{20}c} {} \\ \end{array} 3}} \ge 0 $$
(14)

Here, we introduce the image smoothing process with a \( 5 \times 5 \) Gaussian kernel and \( \sigma = 1 \) to filter the original image data \( I(x) \) in the spatial domain and obtain a more smoothing image data \( \widetilde{I}(x) \) for further gradient calculation, which could diminish the influence of complicated background intensity distributions to the segmentation results.

In this work, we solve the challenging combinatorial optimization problem (14) by its convex relaxation [11, 13, 14]

$$ \begin{array}{*{20}l} {\mathop {\hbox{min} }\limits_{{u_{1,2} (x) \in [0,1]}} } \hfill & {\sum\limits_{{k \in L_{r} }} { < u_{k} ,C_{k} > } + \sum\limits_{{k \in L_{s} }} {\int_{\varOmega } {g(x)\left| {\nabla u_{k} (x)} \right|dx} } } \hfill \\ {\quad \quad \,\,\,s.t.} \hfill & {\quad u_{1} (x) \le u_{2} (x),\quad \forall x \in V} \hfill \\ \end{array} $$
(15)

subject to the geometric constraints (9)–(10).

Obviously, the binary-valued labeling functions \( u_{1,2} (x) \in \{ 0,1\} \) are relaxed into the convex constraints \( u_{1,2} (x) \in [0,1] \) in (15). Given the convex energy function of (15) and the geometric constraints (9)–(10), the challenging combinatorial optimization problem (13) is reduced to the convex relaxation optimization problem [11, 13,14,15,16], as shown in (15).

2.2 Dual Formulation: Continuous Max-Flow Model

In this section, we introduce a new spatially continuous flow maximization model, namely the continuous max-flow (CMF) model [11, 13,14,15,16], which is mathematically equivalent or dual to the proposed convex relaxation model (15). Particularly, it results in an efficient algorithm to compute the optimum labeling functions of (15) while avoiding tackling its challenging non-smoothing function terms and constraints directly [11, 15, 16].

The proposed CMF model can be visually explained by means of streaming the maximum flow from a source through the specified graph, as shown in Fig. 2, and mathematically formulated as follows [11, 15, 16]

$$ \mathop {\hbox{max} }\limits_{p,q} \quad \int_{V} {p_{b} (x)dx} $$
(16)

subject to:

Fig. 2.
figure 2

The flow configuration of the proposed CMF model: links between terminals and the image regions, the source flow \( p_{b} (x) \) and the sink flow \( p_{k} (x),\begin{array}{*{20}c} {} \\ \end{array} k \in L_{r} \).

Full size image

  • Flow capacity constraints [11, 15, 16]: the sink flows \( p_{k} (x),k \in L_{r} \) suffice:

    $$ p_{k} (x) \le C_{k} (x),\quad k \in L_{r} , $$
    (17)

    and the spatial flows \( q_{k} (x),k \in L_{s} \) suffice [11, 15, 16]:

    $$ \left| {q_{k} (x)} \right| \le g(x),\quad k \in L_{s} , $$
    (18)

  • Flow conservation constraints [11, 15, 16]: the total flow residue vanishes at each \( x \) within the image domain \( \varOmega_{k} \):

    $$ \begin{aligned} & G_{IB} (x): = (divq_{IB} - p_{w} + p_{u} )(x) = 0; \\ & G_{OB} (x): = (divq_{OB} - p_{b} + p_{w} )(x) = 0; \\ \end{aligned} $$
    (19)

As defined above, the source flow function \( p_{w} (x) \) is free of constraints. Through analysis, we can prove the duality between the CMF formulation (16) and the convex relaxation model (15), i.e.

Proposition 1.

The CMF model (16) and the convex relaxation model (15) are dual (equivalent) to each other, namely

$$ (15) \Leftrightarrow (16) $$

2.3 Duality-Based CMF Algorithm

According to Proposition 1, it is easy to see that the convex relaxation model (15) can be solved equally by computing the CMF model (16), which derives the new duality-based CMF algorithm proposed in this section through the state-of-the-art global convex optimization techniques. In addition, the introduced CMF algorithm enjoys great numerical advantages [11, 13,14,15,16]: (1) it successfully avoids directly tackling non-smoothing total-variation functions in the energy of convex minimization problem (15) by the projections to some simple convex sets instead; (2) it implicitly applies the anatomic inter-surface prior (2) and geometric constraints (9)–(10) into the introduced flow-maximization scheme (as depicted in Fig. 2).

With reference to the theory of augmented multiplier algorithms [11, 13,14,15,16], an efficient hierarchical CMF algorithm can be derived, which iteratively optimizes the following augmented Lagrangian function,

$$ \mathop {\hbox{max} }\limits_{p,q} \quad \mathop {\hbox{min} }\limits_{u} L_{c} (u;p,q): = L(u;p,q) - \frac{c}{2}\sum\limits_{k \in Ls} {\left\| {G_{k} } \right\|^{2} } ,\quad where\quad u = u_{i} (x),i \in \{ 1,2\} . $$
(20)

subject to the flow capacity constraints (17) and (18). The coefficient \( c \) is positive constant [11, 15, 16]. The \( L(u;p,q) \) above is the Lagrangian function associated with the CMF model, which is defined as

$$ L(u;p,q) = \int_{V} {p_{b} (x)dx} + < u_{1} ,divq_{IB} - p_{w} + p_{u} > + < u_{2} ,divq_{OB} - p_{b} + p_{w} > $$
(21)

Therefore, the CMF algorithm explores the following steps at the n-th iteration:

  • Maximize \( L_{c} (u;p,q) \) over the spatial flows \( \left| {q_{k} (x)} \right| \le g(x),\quad k \in L_{s} \), while fixing the other variables \( (u;p)^{n} \), which amounts to

    $$ q_{k}^{n + 1} : = \mathop {\arg \begin{array}{*{20}c} {\hbox{max} } \\ \end{array} }\limits_{{\left| {q_{k} (x)} \right| \le g(x)}} - \frac{c}{2}\left\| {divq_{k} - F_{k}^{n} } \right\|^{2} $$
    (22)

    where \( F_{k}^{n} \)(x), \( k \in L_{s} \), are fixed. This can be computed by the gradient-projection iteration:

    $$ q_{k}^{n + 1} = \Pr {\text{oj}}_{{\left| {q_{k} (x)} \right| \le g(x)}} (q_{k}^{n} + \tau \nabla (divq_{k}^{n} - (F_{k}^{n} ))); $$
    (23)

    where \( \tau > 0 \), represents the step size for convergence [15, 16].

  • Maximize \( L_{c} (u;p,q) \) over the source flows \( p_{b} (x), \) while fixing the other variables \( (u;p_{k} ,q)^{n} \), which amounts to

    $$ p_{b}^{n + 1} : = \mathop {\arg \begin{array}{*{20}c} {\hbox{max} } \\ \end{array} }\limits_{{p_{b} }} \int_{V} {p_{b} dx} - \frac{c}{3}\sum\limits_{{k \in L_{r} }} {\left\| {p_{b} - J_{k}^{n} } \right\|^{2} } $$
    (24)

    where \( J_{k}^{n} (x) \), \( k \in L_{r} \), are fixed. This can be solved exactly by

    $$ p_{b}^{n + 1} (x) = (J_{{_{b} }}^{n} (x) + J_{w}^{n} (x) + J_{{_{u} }}^{n} (x) + 1/c)/3 $$
    (25)

  • Update the labeling functions \( u_{1,2} (x) \)

    $$ u_{{_{1,2} }}^{n + 1} (x) = u_{{_{1,2} }}^{n} - cG_{1,2}^{n} (x) $$
    (26)

    where \( G_{1,2}^{n} (x) \) stand for the representative flow residue functions.

3 Experiments and Results

3.1 Study Population and Image Acquisition

The database consists of five sets of 3D bladder MR images respectively from five patients. They were consecutively selected for this study at Tangdu Hospital. All of them were postoperative pathologically confirmed bladder cancer and underwent preoperative bladder MRI examination.

Subjects were scanned by a clinical whole body scanner (GE Discovery MR 750 3.0T) with a phased-array body coil. A high-resolution 3D Axial Cube T2-weighted MR sequence was used due to its high soft tissue contrast and relatively fast image acquisition [2, 9]. Before MRI scan, each subject was instructed to drink enough mineral water, and waited for an adequate time period so that the bladder was fully distended. The scanning parameters are listed in Table 1.

Table 1. Principal parameters used for MRI acquisition
Full size table

3.2 User Interaction

The user interaction in our method was choosing some nodes of urine region, bladder wall region and outer background region, respectively, on a single transverse slice of 3D MR image, as shown in Fig. 3(a). The nodes were chosen by the user with polygonous box interface tool, as shown in Fig. 3(b). With the nodes of these three regions delineated, the corresponding PDF functions of these regions were derived to conduct the 2D bladder multi-region CMF segmentation. Figure 3(c) shows the PDF curves calculated of the corresponding regions.

Fig. 3.
figure 3

Example of user initialization process for PDF functions calculation. (a) Is the single transverse slice of the 3D MR image. (b) Is the initialization contours manually depicted by the user. (c) Is the calculated PDF functions of the corresponding regions.

Full size image

3.3 2D Multi-region Segmentation

With the PDF functions of three regions, the cost function \( C_{u} (x) \), \( C_{w} (x) \) and \( C_{b} (x) \) were calculated based on (11). Here, we introduce the hard constraints [11, 15, 16]:

$$ \left\{ {\begin{array}{*{20}c} {C_{w} (x) = + \infty ,\quad C_{b} (x) = + \infty ,x \in N_{u} } \\ {C_{l} (x) = + \infty ,\quad C_{b} (x) = + \infty ,x \in N_{w} } \\ {C_{l} (x) = + \infty ,\quad C_{w} (x) = + \infty ,x \in N_{b} } \\ \end{array} } \right. $$
(27)

where \( N_{u} \), \( N_{w} \) and \( N_{b} \) denote the nodes manually delineated from the urine, bladder wall and background regions, respectively.

In order to further reduce the cause of the complicated intensity distributions in the background region, we define a bounding box that roughly removes the sub-background regions with certain range of distance from the bladder wall region, as shown in Fig. 4(a).

Fig. 4.
figure 4

2D multi-region segmentation results. (a) Is bounding box placed around the outer-wall. (b) Shows the segmentation results. (c) Illustrates the convergence curve during iterations.

Full size image
$$ C_{l} (x) = + \infty ,\quad C_{w} (x) = + \infty ,x \in N_{Boundingbox} $$
(28)

The urine region usually appears to be homogenous in intensity distributions and has more differentiable edges at IB, whereas the bladder wall is somewhat heterogeneous in intensity distribution and even has overlapping intensities with the background region. Consequently and empirically, we adopt the minimum bladder wall thickness (BWT) prior \( d_{\hbox{min} BWT} \) as 3 mm, such that the minimum distance between IB and OB should be larger than \( d_{\hbox{min} BWT} \) to maintain the separation of these two surfaces globally.

$$ C_{w} (x) = + \infty ,\quad s.t.\quad \hbox{min} \,d(x,R_{u} ) < d_{\hbox{min} BWT} $$
(29)

The 2D segmentation results exhibit in Fig. 4(b), and the convergence error curve during the entire iterations is shown in Fig. 4(c).

3.4 3D Multi-region Segmentation

With the single slice segmentation results, a region growing segmentation [11, 17] algorithm can be applied to obtain an initial 3D IB surface. To obtain the initial 3D surface of OB, the initial 3D IB surface can be dilated empirically by 8 mm. Finally, the two initial surfaces were simultaneously evolved using the proposed 3D multi-region segmentation approach for 3D bladder wall segmentation. Figures 5 and 6 show the final 3D segmentation results.

Fig. 5.
figure 5

The final segmentation results

Full size image
Fig. 6.
figure 6

The 3D bladder constructed with the segmentation results

Full size image

3.5 Validation

The Dice’s similarity coefficient (DSC) is used to quantitatively evaluate the accuracy of the algorithm with manual segmentation as the ground truth [11]. Besides, the performance of using the 3D coupled CMF algorithm is compared with that of level-set method [6, 9]. The results are shown in Table 2. From the table, we can find that the proposed MR-CMF algorithm outperforms Level-set method in bladder wall segmentation. Figure 7 is the results comparison of these two methods.

Table 2. Accuracy and efficiency of MR-CMF and level-set in bladder segmentation
Full size table
Fig. 7.
figure 7

The result comparison of using Level-set and MR-CMF for bladder segmentation

Full size image

4 Discussions and Conclusions

We described and verified a novel global optimization-based approach for simultaneously segmenting the multiple regions of 3D bladder MR images by constructing two coupled surfaces of IB and OB via 3D MR-CMF algorithm, incorporating the anatomic inter-surface constraints with minimum BWT prior to jointly enhance the segmentation accuracy, efficiency and robustness, especially in the wall region that has overlapping intensity values with the background.

One of the obvious advantages of this new approach over the previous bladder wall segmentation algorithms is that it is not very sensitive to the initialization process, especially the background initialization [11, 15, 16]. The main reason is that the PDF model generated from the background nodes initially delineated from the 2D slice is just used for 2D segmentation of IB and OB. When these two boundaries are obtained, the PDF model corresponding to the current background is iteratively computed, eliminating the initial model.

Another advantage is that the mathematical schemes of MR-CMF are based on globally optimization techniques [11, 13,14,15,16], whereas the previous bladder segmentation algorithms are based on local optimization techniques, which are more sensitive to the user initialization and image quality.

The results in Table 2 and Fig. 7 demonstrate that the proposed 3D MR-CMF algorithm outperforms the Level-set based 3D methods in bladder MR image segmentation, which both has a higher accuracy and a faster calculation speed to meet the need for clinical applications. The proposed method was running on a platform with Intel (R) Core (TM) i5-6400 CPU @ 2.7 GHz, 8.00 GB RAM, equipped with 64-bit Windows 7, and MATLAB 2015 b. The average time expense of five datasets is 6.0 min per each, whereas the average time consumption that of using Level-set approach is 37.9 min.

However, the datasets used to verify this new algorithm is very sparse, more datasets need to be enrolled to revalidate the approach and come up with more objective and convincible results. Besides, in current version of MR-CMF algorithm, we only adopt the manually defined cost functions derived from the PDF functions of the intensity distributions in each targeted region. In the future work, we’ll consider incorporating deep learning convolutional neural network (DL-CNN) into cost functions determination process to dig out some more useful factors as cost functions to further enhance the segmentation results.