Unsupervised 3D Prostate Segmentation Based on Diffusion-Weighted Imaging MRI Using Active Contour Models with a Shape Prior

the


Introduction
Prostate cancer is the second leading cause of cancer-related deaths and most frequently diagnosed cancer in American men [1].Therefore, there is a significant interest in improvements of prostate cancer diagnosis and treatment.Imaging methods that could provide reliable information about the location, size, and shape of prostate gland would greatly useful to localize cancer foci, guide biopsies, and radiotherapy.To this date, the most widely used modality for prostate cancer diagnosis is trans rectal ultrasound (TRUS) because of its low cost and short acquisition time.However, its false negative rate is high [2], and prostate cancer visualization is poor.As an alternative, high-resolution MRI allows physicians to better evaluate the prostate diseases that may not be assessed adequately with other imaging methods such as X-ray, TRUS, and computed tomography (CT).Recent studies have shown that MRI has higher accuracy in the detection of prostate cancer [3].Because of the advances in MRI technology, diffusion-weighted (DWI) MRI is also now commonly applied to the prostate along with other MRI techniques.Prostate volume is routinely asked as part of imaging evaluation, as it helps in clinical decision making when combined with serum prostate-specific antigen (PSA) to derive PSA density.Knowledge of prostate boundaries is useful in the planning of conformal radiation therapy and computer-aided prostate cancer localization.Although the identification of prostate boundary is a crucial step in these clinic applications, manual segmentation prostate boundaries on 3D MR images slice by slice is a tedious and laborious job.Moreover, the manual segmentation is subjective and produces different results among different observers.Therefore, accurate automated prostate segmentation based on 3D MR images is extremely useful.However, this task is very challenging, because of the noise and inhomogeneity of MR images and the complex anatomical structures of the prostate and surrounding organs.
In the literature, previous work on automated prostate segmentation is primarily focused on TRUS images.Pathak et al. proposed an edge-based boundary delineation scheme to detect prostate edges as a visual guide to the observer doing manual editing in [4].Because of the visual guide, the accuracy of the detected prostate edges was as good as those of the human observers.Shape information of the prostate were also incorporated in the literature to improve the segmentation performance.In [5,6], ellipses were used to model the prostate shape.In [5], the prostate shape was modeled by parametric deformable superellipses, and Bayesian segmentation algorithm was then applied to 2D TRUS images.In [6], an elliptical level set algorithm was proposed to segment prostate with TRUS images.Due to the shape information, accurate and consistent segmentation results could be obtained in 2D without any manual intervention.Statistical shape information extracted by Gabor filter from training dataset was employed in [7].Besides the global populationbased shape statistics, Yan et al. proposed to combine patientspecific local shape information to segment prostate with TRUS video in [8] and further improved the results proposed in [5,7].The shape-based prostate segmentation method was also applied to CT images.Freedman et al. presented a method based on matching probability distributions of photometric variables that incorporates learned shape and appearance models for the prostate and applied it to 3D CT images in [9].
Compared with research in automated prostate segmentation using TRUS images, attempts on MR images are limited.For most of the literature available algorithms developed for prostate segmentation with MRI, the shape information were widely considered.One way to incorporate the shape information is to learn the shape statistics from training dataset.Tsai et al. derived a model-based, implicit representation of the segmentation curve evolution by applying principle component analysis (PCA) to a set of signed distance representations of the training data.This method is applied for the segmentation of medical images containing known object types [10] and could obtain satisfactory visual results of prostate volume.In [11], authors did not only extract the shape information, but also the texture information of the prostate region by PCA to further improve the segmentation performance.Statistical atlas was also used to incorporate the shape information.Klein et al. developed an atlas matching method to segment prostate from MR images based on prelabeled and registered atlas images [12].In [13,14], two semiautomatic prostate segmentation methods were proposed.In [13], a method using wavelet multiscale products to detect the prostate boundaries were developed.This method requires the user to specify four reference points around the prostate.In [14], the prostate contour of one slice was manually refined and used as initial estimation in the neighboring slices.The contour was propagated in 3D through steps of refinement in each slice.Template matching was also used to fuse the prostate shape information.Based on reasonable initials, these two algorithms could successfully segment the prostate.Although majority of the literature available methods for prostate segmentation with MRI are based on T2 MRI which provides details of the prostate structures, there are a few attempts to consider other MR techniques to perform prostate segmentation.In [15], T2 MRI and MR spectroscopy were combined, and an active shape model segmentation scheme was developed.An appropriate initialization is essential to the accurate segmentation for this method.In [16], a framework based on maximum aposterior (MAP) estimation was proposed to segment prostate from dynamic contrast enhanced MRI (DCE MRI) by fusing appearance and spatial and shape information of the prostate learned from training data.
To this date, all these available prostate segmentation methods with MRI are either supervised or semiautomatic, and the supervised methods have difficulties to handle the large variety of the prostate size, shape, and texture of different patients.Since the prostate appearance varies significantly between patients, we develop an unsupervised segmentation algorithm in this work based on level set framework introducing a shape prior to region-based active contour model.The proposed method is fully automatic, and it segments prostate from apparent diffusion efficient (ADC) images derived from diffusion-weighted imaging (DWI) MRI.The major contribution of this paper is the development of a 3D automated prostate segmentation method which does not need training data and is the first attempt to make use of DWI MRI to differentiate the prostate region with other tissues.
Implicit level set-based representations of a contour have become a popular framework for medical image segmentation.The question of how to fuse higher-level shape prior information into level set-based contour evolutions has been addressed by a number of researchers.In many of the previous study, either the shape prior is extracted from training data [5,7,8,10,11,[16][17][18], or the exact shape of the object is assumed to be known [19,20].In this paper, we propose a novel approach to obtain the shape information from the 3D MR images.The proposed method is a level set-based active contour model which incorporates shape information by adding a shape penalty term.The idea of adding a shape penalty term is given in [19].However, instead of having the exact shape of the object as in [19,20], or learning the shape from training data, we use a three-step strategy [21]: (i) coarsely segment the prostate volume, (ii) then based on the coarse segmentation result, the prostate shape is modeled by a series of deformable ellipses slice by slice to constrain the level set evolution as close to an ellipsoid as possible, (iii) estimate the prostate volume by region-based contour model and shape prior defined by the previous step.To incorporate the shape prior into active contour model, we introduce a shape penalty term to the energy functional and propose a method to automatically select the shape penalty weight.Finally, a series of morphological operations are applied to further refine the prostate boundary.The proposed method is based on our previous study [21], where the prostate was segmented from ADC images in 2D.In this paper, we improve that method in the following aspects: (i) we extend it to 3D from 2D, (ii) we propose a coarse segmentation step and use a stack of parametric deformable ellipses to extract the shape information, (iii) we develop a method to automatically select the weighting parameter of the shape term, and (iv) we apply morphological operations to refine the segmentation results.In our experiment, ADC images were calculated from diffusion-weighted data on a pixel-wise basis, according to ADC = −1/b ln S/S 0 , where S 0 is the signal intensity without the diffusion weighting (b value of 0 sec/mm 2 ) and S is the signal intensity with the gradient (b value of 600 sec/mm 2 ).ADC images measure the diffusion coefficient and use diffusion as contrast.
This paper is organized as follows.In Section 2, the proposed prostate segmentation scheme for 3D MR images is explained in detail, including the basic concepts of regionbased active contour model, the parametric deformable ellipsoid model, the proposed segmentation algorithm with shape information, and the automated shape penalty weight selection method.Section 3 provides the experimental results of applying the proposed method to 3D prostate MR images and comparison with the manual segmentation results.A summary and conclusion of our prostate segmentation method is given in Section 4.

Segmentation Method
In this section, we explain the proposed segmentation method in detail.The proposed method is based on a level set formulation of the Mumford-Shah functional developed by Chan and Vese.We extend this framework by introducing a shape penalty term to constrain the level set evolution.Our input data is a 3-D apparent diffusion coefficient (ADC) maps calculated from diffusion-weighted (DWI) MR prostate dataset.As shown in Figure 1, our method consists of four main steps: (i) a coarse segmentation step to roughly obtain the prostate shape; (ii) a shape information extraction step to estimate the shape of the prostate, (iii) segmentation step to estimate the prostate volume by region-based active contour model combining the shape prior, and (iv) a refining step to smooth the prostate surface and remove the isolated components in the segmentation result.Each of these steps is described in detail next.

Coarse Segmentation by Region-Based Active Contour
Model.In the first step, we use a region-based active contour model to 3D ADC images to obtain a coarse prostate shape to further extract the shape information in the next step.For medical images, including prostate MR images, the tissue of interest may not have complete boundaries, or have complex anatomical structures.The edge-based segmentation method, including active contour with edges model [22], also named geodesic active contours, largely depends on the nearby edges, is sensitive to local minimum and noise and cannot deal with topological changes.Because of these shortages, we consider region-based active contour model in our application.Compared with edge-based models, regionbased models consider the pixel intensities within the entire image dataset [10,23,24].The image dataset is segmented into a certain number of regions based on the regional statistics (sample mean and variance) of the corresponding region.Therefore, region-based active models are more robust to noise and can handle topological changes.In [23], Chan and Vese proposed a pure region-based model to segment image where u is the segmentation image and c 1 and c 2 are the average intensities of the two regions partitioned by the curve C.During the minimization of (1), the image is divided into two regions: inside and outside of the curve.Level set framework is combined to minimize the energy function shown in (1).Level set method first introduced in [23] is a numerical technique that can follow the evolution of interfaces.It has been applied to various image processing applications including image segmentation, reconstruction, and denoising.Chan and Vese's region-based active contour model combines the Mumford-Shan functional and level set framework.The level set formulation of variational active contour model is based on a higher-dimensional level set function ϕ, whose zero level set segments the image into several intensity homogeneous regions.The energy function of Chan-Vese model in level set framework in 3D is:  where u is the segmentation image, H(•) is the Heaviside function, and c 1 and c 2 is the average intensities (sample means) of the two regions segmented by zeros level set ϕ(x, y, z, 0) Figure 2 shows the coarse segmentation results obtained by the region-based active contour model for one patient with different slices in mid-gland, near apex, and base.These results are very poor and not acceptable.We can also see that the intensity information is not sufficient to distinguish the prostate gland from surrounding organs and tissues.Therefore, it is necessary to combine the shape information of the prostate to improve the performance of automated segmentation method.

Shape Information Extraction.
Medical image segmentation in general faces difficulties because of noise, missing boundaries, and complex anatomical structures.Under such conditions, introducing some prior information, such as the general shape, location, intensity, and curvature profile of the tissue of interest could help the segmentation algorithm perform better.In prostate MR images, the prostate gland and surrounding tissues, such as the bladder, rectum, and muscles, have overlapping intensity and texture.For some patients, at certain slices, the prostate boundaries may be missing or blended with those surrounding tissues.However, the prostate has a walnut-like shape in general.Combining this shape information that the prostate is close to an ellipsoid in 3D could constrain the segmentation algorithm evolution and help it extract the prostate more accurately.

Parametric Deformable Ellipsoid Model.
After the coarse segmentation step, we model the prostate shape by a stack of parametric ellipses.In the literature, several methods have been proposed for fitting superellipses [25].In this study, to obtain the prostate shape information, we use a twostep scheme.First, we roughly model the prostate volume by a parametric ellipsoid.Then, we model the apex as a stack of parametric ellipses.Since the prostate volume is not an ideal ellipsoid, we use a stack of parametric ellipses to model the apex to fit the prostate more accurately.In the first step, the prostate volume is fitted by a parametric ellipsoid roughly as follows: where (a, b, c) is the center of the ellipsoid, (r 1 , r 2 , r 3 ) the lengths of the semiaxes, and θ the orientation.Shape parameters v = (a, b, c, r 1 , r 2 , r 3 , θ) define an ellipsoid.In (4), only one rotation is considered.By observing the axial MR images, we can see that the rotation of the prostate in the axial plane is ignorable, so we assume there is no rotation in x-y plane in (4).To obtain an ellipsoid which best fit the prostate shape, we borrow the idea of least-square minimization and superquadric inside-outside function presented in [26].
Based on the implicit representation of the parametric ellipsoid, we define the function called the inside-outside function.When F(x, y, z) = 0, the corresponding voxels are on the surface of the ellipsoid.
When F(x, y, z) > 0, the corresponding voxels are inside of the ellipsoid, and vice versa.To find a function with a minimum corresponding to the ellipsoid that best fits the given prostate shape, we define a shape fitting function: where Ω is the image domain and H(ϕ) is the prostate region obtained by the coarse segmentation step.For voxels inside of the ellipsoid, we have H[F(x, y, z)] = 1; for voxels outside of the ellipsoid, we have H[F(x, y, z)] = 0.The shape parameters are obtained as If the prostate is perfectly segmented, H(ϕ) is the ideal prostate mask, the deformable ellipse converges to the smallest ellipsoid which best fits the prostate volume.In practice, we simplify the estimation of the ellipsoid by estimating c ≈ r 3 = N/2, where N is the number of slices containing prostate region and is predefined by a radiologist.That means, a radiologist first selects the slices which belong to the prostate region.If 18 slices are selected, then we have c = r 3 = 9.In this way, the shape parameters of the ellipsoid which best fits the prostate shape in 3D can be approximately estimated in 2D by finding an ellipse which best fit the prostate shape in the central slice, where That is for the ellipse v ct = (a, b, N/2, r 1 , r 2 , N/2, θ), we have The cost function is minimized by iterative gradient descent method, and the gradient descent with respect to the unknown shape parameters v ct = (a, b, r 1 , r 2 ) is where By observing the ellipsoid fitting results, we can see that the ellipsoid as shown in Figure 3 (corresponding to the same prostate images in Figure 2) is able to roughly model the prostate shape in 3D.To further improve the shape result, we apply this ellipse fitting method to the slices of the apex again to obtain a stack of ellipses fitting the prostate shape more accurately as shown in Figure 4 (corresponding to the same prostate images in Figure 2).That is, we update the shape information by finding an ellipse that best fits each slice, and for each slice, we obtain a set of shape parameters defining an ellipse as follows: and those ellipses are combined with the active contour model.It is worth to mention that at this step, the apex slices need to be predefined.Although the identification of the apex is difficult, it is not crucial in this step.If the slices of the mid-gland are misidentified as apex, the results of the deformable ellipses will not change, because for the midgland slices, the shape has already been fitted by the ellipsoid very well.

Initial Estimates of the Shape Parameters.
It is worth to mention that the gradient descent minimization may   converge to a local minimum instead of a global minimum.Therefore, initial estimates of the set of shape parameters v determines to which local minimum the minimization procedure will converge.In the proposed method, we use a rough estimation of the prostate's true position, orientation and size obtained based on the profile (the intensity values) across the centroid of the coarse segmentation result along vertical and horizontal direction of the ADC images.These initial estimates suffice to assure convergence to the minimum that corresponds to the actual shape of the prostate.After the coarse segmentation step described in Section 2.1, we calculate the number of pixels of the coarse segmentation result of the central slice along the horizontal and vertical direction.We can see that the profile image has roughly a rectangular shape, and we can detect the rectangular edges which corresponding to the prostate boundary by calculating the first derivative of the profile image.By detecting the left, right, top, and bottom boundary of the prostate in the central slice, we can calculate the center and radius of the prostate as initial estimates of the shape parameters to assure a more robust shape extraction.

Active Contour Model with Shape
Prior.There are several ways to incorporate the shape information into level setbased variational approaches.In [10], a number of training shapes are implicitly represented in the segmentation curve using signed distance functions.In [19,20], authors proposed two models to introduce shape priors into Chan-Vese models, but their models are both based on the exact shape of the object is known and segment the known shape or object from the background, where there are several objects.However, in our application, the exact shape of the prostate is unknown and considering the large variety of the prostate shapes and sizes, we propose an unsupervised shapebased segmentation algorithm.We assume the prostate shape is close to an ellipsoid which is estimated by the method discussed in Section 2.2.To combine the shape prior and constrain the level set evolution, we add the shape fitting function (6) to the level set energy function (2) as a penalty term as follows: where β is a weighting parameter of the shape fitting term.The shape penalty term forces the segmented region H(φ) close to the shape prior H(F).We update the level set function by gradient descent method, and the gradient descent with respect to the segmentation function φ is 2.3.2.Shape Weighting Parameter Selection.Usually, the shape weighting parameter is selected manually based on previous experience.In this paper, we present a method to select the shape weight β automatically based on the correlation between the segmentation result and shape prior.The correlation R between the segmentation result and shape prior is defined as where W is the segmentation result with shape weight β, F is the shape prior, and i represents the ith voxel.By varying the shape weighting parameter β, we plot the curve correlation R versus shape weight β shown as Figure 5.An appropriate β should be small enough so that the segmentation result will be able to capture the prostate real boundary.Meanwhile, the appropriate β should be large enough so that the shape prior could constrain the level set function evolution, and segmentation result will close to the shape prior.Considering these two points, we compare all the correlation values, and select the one closest to the upper left corner in the plot.That β corresponds to the smallest shape weight with high correlation between the segmentation result and shape prior.
where d is the Euclidean distance, β max is the smallest shape weight for the correlation R ≥ 0.98.
Figure 6 (corresponding to the same prostate images in Figure 2) shows the different segmentation results obtained by the active contour model with different weight of shape prior.This figure demonstrates the efficacy of our method of selecting appropriate shape parameter β automatically.

Prostate Volume Refinement.
After the segmentation with shape prior, the prostate volume is obtained.However, certain surrounding tissues are also labeled as prostate and appear as some isolated components in the image data.To remove those tissues, a morphological opening operation is firstly applied, and then, only the largest component in the image domain which corresponds to the prostate volume is selected.Finally, a morphological closing is used to restore the prostate boundaries detected.The results are shown as in Figure 7 (corresponding to the same prostate images in Figure 2), which is our final segmentation result.
The main steps of the presented approach can be summarized as (1) Apply a region based active contour model to the 3D ADC image data to obtain a coarse estimation of the prostate mid-gland and apex.
(2) Estimate the prostate shape by using a parametric deformable ellipse model based on the coarse segmentation of the prostate mid-gland and apex.
(3) Apply region-based active contour model again with shape prior obtained in the previous step to further segment the prostate volume in 3D with an automatically selected weighting term.
(4) Apply morphological processing step to refine the prostate volume result.

Experimental Results
In this study, MR image data obtained from ten patients with biopsy-confirmed prostate cancer are used.After the prostatectomy, the prostate was removed and weighted by a pathologist.Because ADC maps provide better anatomical shape and contrast between the prostate gland and other tissue, we apply our method to 3D ADC images.[27], and dice measure (DSC) to evaluate our segmentation scheme.We denote the manual delineated boundary as Q = {q 1 , q 2 , . . ., q η } and automated segmentation results as W = {w 1 , w 2 , . . ., w(ι)}, where each element of W and Q is a point on the corresponding contour.We find the distance of every point in W from all points in Q.We define the distance to the closest point for w j to the contour Q as where • is the 3D Euclidean distance between any two points.The Hausdorff distance is defined as the maximum d(w j , Q) over all j.The MAD is the average of d(w j , Q) over all j.The Hausdorff distance measures the worst possible disagreement between the two boundaries, while the MAD estimates the disagreement averaged over the two outlines.
On the other hand, the DSC value is defined as where W is the automatic segmentation result, Q the manual segmentation by an expert radiologist, and | • | denotes the number of voxels contained in the set.
Figure 8 provides a comparison between the proposed method and the manual segmentation results and Table 1 shows the DSC, MAD and Hausdorff distance values of ten patients.The weights of the prostates are also provided in Table 1, and we can see that the size of the prostate varies significantly among patients.In Table 2, the segmentation results of the base, mid-gland and apex are provided separately.We can see that the majority of mis-segmentation occurs at the base and apex where the surface between prostate and surrounding tissues are very weak.Comparison of segmentation results with other prostate segmentation schemes in the literature show that our system performs at least as well as, or better than other systems.Klein et al. [12] have reported a mean DSC value of 0.82, Zhu et al. [28] have DSC values ranging from about 0.15 to about 0.85, and Toth et al. [15] have DSC values ranging from 0.746 to 0.826, while our DSC values range from 0.738 to 0.871.Note that our method does not need a training dataset.

Conclusions
In this study, we have developed and applied an unsupervised automated segmentation method to the problem of prostate segmentation with 3D DWI MR image data.Accurate segmentation of prostate from MR datasets is useful in many applications.Although many researchers have proposed algorithms for prostate segmentation, attempts on MR prostate segmentation are very limited with only supervised techniques that require a training dataset.Currently, the level set framework is a popular approach for medical image segmentation, and shape information is also considered in most prostate segmentation method presented in the literature.To this date, in most shape-based prostate segmentation methods, either for TRUS or MR images, the shape information is obtained from training data or by compared with atlas images, but the prostate shape, size, and texture vary widely between patients.Besides, in the literature, the majority of MR-based prostate segmentation algorithms are based on T2 MRI.In this paper, we present an unsupervised and automated method to segment prostate volume based on DWI MRI by a shape-based active contour model with level set framework without the need of a training dataset.
We extend the region-based active contour model proposed by Chan and Vese and apply it to MR images by fusing a shape penalty term to the cost function.We firstly apply a coarse segmentation step to 3D ADC image data, and we model the prostate shape by a stack of parametric deformable ellipses to extract the shape prior information.Then, we introduce a shape fitting function to force the active contour evolution close to the shape prior for further segmentation, and we select the shape weighting parameter automatically, as explained in Section 2. The experimental results on 3D MR prostate images show the effectiveness of the proposed method.
Because of the high variability of the prostate appearance between patients, future work will include applying our method to a larger MR dataset.Because of the nonuniformity of the texture and the lack of clear edge of the prostate apex and base, our method performs poor at certain slice for certain patients, future work will also attempt to overcome these limitations.
(a) near base (b) in mid-gland (c) in mid-gland (d) near apex

Figure 2 :
Figure 2: Coarse segmentation results by active contour model overlaid on the ADC images.The prostate was removed after surgery, and the specimen weighs 45 g and measures 4.5 cm SI × 4.0 cm ML × 3.2 cm AP, part(a) is 3 mm to the actual base, and part(d) is 3 mm to the actual apex.

Figure 3 :
Figure 3: The prostate shape can be roughly fitted by an ellipsoid.

Figure 4 :
Figure 4: The prostate shape is fitted by a stack of ellipses.The first row are the ellipses, and the second row is the prostate outlined by a radiologist.

Figure 6 :Figure 7 :
Figure6: Segmentation results comparison with different shape weight.The first row is the segmentation results obtained without shape prior (β = 0), the second row is the segmentation results obtained with a very large shape weight (β = 10 6 ), the third row is the segmentation results obtained with an automatic selected shape weight.The last column shows the segmentation results in 3D.

Figure 8 :
Figure8: 3D segmentation results of three patients.The first row is the manual segmentation by an expert human reader, and the second row is the automated segmentation by proposed method.The prostate specimen corresponding to the first column weighs 45 g and measures 4.5 cm superior to inferior (SI) × 4.0 cm medial to left (ML) × 3.2 cm anterior to posterior (AP), the prostate specimen corresponding to the second column weighs 72.1 g and measures 5.5 cm SI × 5.0 cm ML × 3.6 cm AP, the third column's prostate specimen weighs 56.3 g and measures 4.5 cm SI × 5.0 cm ML × 3.2 cm AP.
The segmentation typically takes about 7 minutes for each patient (about 256 × 256 × 18) on an Intel Core2 Quad PC running at 2.4 GHz.We use edge-based and volume-based metrics measurements for quantitative analysis of segmentation results: Hausdorff distance, mean absolute distance (MAD)

Table 1 :
DSC, MAD, and Hausdorff distance values of ten patients.

Table 2 :
The mean and standard deviation values of DSC, MAD, and Hausdorff distance of the whole prostate and at the base, mid-gland and apex of ten patients.