A diffusion-weighted imaging based diagnostic system for early detection of prostate cancer

A new framework for early diagnosis of prostate cancer using Diffusion-Weighted Imaging (DWI) is proposed. The proposed diagnostic approach consists of the following four steps to detect locations that are suspicious for prostate cancer: 1) In the first step, we isolate the prostate from the surrounding anatomical structures based on a Maximum A Posteriori (MAP) estimate of a new log-likelihood function that accounts for the shape priori, the spatial interaction, and the current appearance of prostate tissues and its background (surrounding anatomical structures); 2) In order to take into account any local deformation between the segmented prostates at different b-values that could occur during the scanning process due to local motion, a non-rigid registration algorithm is employed; 3) A KNN-based classifier is used to classify the prostate into benign or malignant based on three appearance features extracted from registered images; and 4) The tumor boundaries are determined using a level set deformable model controlled by the diffusion information and the spatial interactions between the prostate voxels. Preliminary experiments on 28 patients (17 malignant and 11 benign) resulted in 100% correct classification, showing that the proposed method is a promising supplement to current technologies (biopsy-based diagnostic systems) for the early diagnosis of prostate cancer.


INTRODUCTION
Prostate cancer is a major health problem, and the most frequently diagnosed malignancy in the American male population [1].Recent prostate cancer studies reported an estimated 241,740 new cases and a mortality rate of close to 28,170 in 2012 [2].Fortunately, early diagnosis of prostate cancer increases the survival rate of the patients [3].

Current Imaging Modalities for Prostate Cancer Diagnosis
Currently, there are different techniques that are used for early diagnosis of prostate cancer.However, the accuracy of these techniques are clearly unsatisfactory.For example, Prostate Specific Antigen (PSA) screening doesn't offer accurate information about the location and extent of the lesion(s) [4].In addition, PSA is associated with a high risk of over diagnosis of prostate cancer.On the other hand, imaging tests using different imaging modalities, such as Transrectal Ultrasound (TRUS) [5], Computed Tomography (CT) [6], MR Spectroscopy (MRS) [7], Dynamic-Contrast Enhanced Magnetic Resonance Imaging (DCE-MRI) [8], and Diffusion-Weighted Imaging (DWI) [9] are still critically needed.TRUS imaging [10] is widely used for guided needle biopsy due to the real time nature of the imaging system, ease of use, and portability.However, TRUS images have low Signal-to-Noise Ratio (SNR) making it difficult to detect malignant tissues [11].Another traditional imaging modality is CT.It is widely used for diagnosis and followup of prostate cancer [12].However, it has poor softtissue contrast resolution that does not allow precise distinction of the internal or external anatomy of the prostate and thus CT images have shown limited specificity for prostate diagnosis [13].
On the other hand, MR image-based modalities, such as T2-weighted MR, MRS, DCE-MRI, and DWI, have also been widely employed for early detection of prostate cancer [14].Despite widely use of T2-weighted MR im-aging for prostate cancer, the technique is limited by unsatisfactory sensitivity and specificity for cancer detection and localization [15].To improve the diagnostic performance of MR imaging in evaluations for prostate cancer, various other techniques have been applied.MRS provides metabolic information about prostate tissue by demonstrating the relative concentration of chemical compounds.However, MRS has its own limitations, such as the need of an additional software and longer acquisition time [16], which lead to increased costs and decreased throughput.Furthermore, MRS suffers from lack of spatial resolution.In addition, signal from periprostatic fat and seminal vesicles can distort spectral waveforms [14].DCE-MRI has been recently suggested for improved visualization and localization of the prostate cancer [17].It provides valuable pathologic and anatomical information.However, DCE-MRI has the drawback of intravenous contrast agent (e.g., gadolinium) administration which is harmful to the kidney [18] and requires a longer setup time.
Recently, DWI has emerged as an imaging modality that has shown more capabilities in determining the size and the shape of the prostate gland and localizing the cancer foci [19].DWI is non-contrast functional imaging technique, whereby the image contrast is determined by the random microscopic motion of water protons, i.e., the Brownian motion [19].Moreover, DWI has the distinct advantage of being acquired very rapidly, without the use of any intravenous contrast material or specialized hardware, and this is the main motivation behind this work.

Clinical Studies for Prostate Cancer Diagnosis Using DWI
In recent years, a growing number of clinical studies [19][20][21][22][23][24][25][26][27][28][29][30] have evaluated the utility of DWI, either in combination with or in comparison with other MRI techniques, for the detection of prostate cancer.These studies have reported various sensitivities and specificities of cancer diagnosis.
Earlier studies [19,20] have investigated the abilities of DWI for prostate cancer diagnosis using an endorectal coil.However, the reported results demonstrated low diagnostic sensitivity.To increase the sensitivity of diagnosis, Shimofusa et al. [21] suggested the addition of strong magnetic field gradient pulses (b-values) to the pulse sequence instead of using endorectal coil.In their study [21], they detected the tumor in the central zone of the prostate in five of eight total patients using DWI with strong magnetic field gradient pulses.Alternatively, the compared diagnostic results with T2-weighted imaging, detected the tumor only in one of the eight patients.Since then, DWI was used for the detection of cancerous tissue in later studies [22][23][24][25][26][27][28][29][30].For example, Tan et al. [30] compared the performance of T2-weighted MRI, DCE-MRI, and DWI for the detection of cancer within the prostate gland.In their study they reported that DWI alone showed better specificity than DCE-MRI alone.It is also showed better overall specificity than combined DWI and T2-weighted imaging.
To the best of our knowledge, there are a very limited number of image-based approaches for automated computer-aided diagnosis of prostate cancer using DWI.These related works are discussed in the following section.

Image-Based Computer-Aided Diagnostic (CAD) Systems for Prostate Cancer Detection
In literature, a limited number of CAD systems for prostate cancer diagnosis have been proposed.For example, Chan et al. [31] proposed an in-vivo CAD system using multimodal MRI to estimate malignancy likelihood in the peripheral zone.They constructed statistical maps from T2-weighted MRI, DWI, and Proton Density (PD) images.These statistic maps were combined with textural and anatomical features of prostate cancer areas in order to detect the cancerous regions.However, this study doesn't include benign regions.Huisman et al. [32] developed a CAD system for prostate lesion classification using a Hessian-based blob detection algorithm [33].
Results showed an accuracy of 92% in classification within the peripheral region and an accuracy of 83% in classification within transitional zones of the prostate.However, their study focused on the peripheral and transitional zones of the prostate gland and excluded central zones in which up to 30% of prostate cancers can occur.Viswanath et al. [34] generated similar likelihood maps by combining information from multimodal MR images using mathematical descriptors.Their study showed, on a voxel basis, that the discrimination between benign and malignant tissue is feasible with good performances.The unsupervised classification by kmeans clustering achieved an accuracy of 77%.Unfortunately, the corresponding slice still needs to be selected between different modality.A study by Langer et al. [35] focused on the peripheral zone of the prostate gland and excluded the central and transitional zones.However, detailed anatomic studies have suggested that 70% of cancers arise in the peripheral zone of the prostate, but up to 30% of prostate cancers occur between transition zones and the central zone of the prostate [36].
To increase the sensitivity of diagnosis, accurate delineation of the prostate region is mandatory.Basically, manual outlining of the prostate borders is the most accurate segmentation that enables precise determination of the prostate volume.However, it is prohibitively time consuming and is prone to intra-and inter-observer variability.Traditional edge detection methods (e.g., [39]) are unable to extract the correct boundaries of the prostate since the gray-level distributions of the prostate and the surrounding organs are hardly distinguishable.Therefore, other automated segmentation methods are desirable.However, multiple challenges stemming from 1) the large variations of prostate shape within a specific time series as well as across subjects; 2) lack of strong edges and diffused prostate boundaries; and 3) the similar signal-intensity profile of the prostate and surrounding tissues, complicates the segmentation process.
The most successful known approaches (e.g., [37][38][39][40][41][42][43][44][45]) have addressed the segmentation challenges of the prostate by modeling object appearance and shape.In particular, Zhu et al. [40] used a combination of an Active Shape Model (ASM) and 3D statistical shape modeling to segment the prostate.Toth et al. [41] presented an algorithm for the automatic segmentation of the prostate in multi-modal MRI.Their algorithm starts by isolating the Region-Of-Interest (ROI) from MRS data.Then, an ASM within the ROI is used to obtain the final segmentation.Klein et al. [42] presented an atlas-based segmentation approach to extract the prostate from MR images.The segmentation of the prostate is obtained as the average of the best-matched registered atlas set to the test image (image to be segmented).Recently, Vikal et al. [43] used a priori knowledge of prostate shape to detect the contour in each slice and then refined them to form a 3D prostate surface.Martin et al. [44] developed an atlas-based approach for segmenting the prostate from 3D MR images by mapping probabilistic anatomical atlas to the test image.The resulting map is used to constrain a deformable model-based segmentation framework.

Current Limitations and Motivation for Our Proposal
The above-mentioned CAD systems for analyzing DWI are not sufficiently accurate and reliable for several reasons: 1) The majority of CAD systems used multimodal MRI which is cost inefficient [45].
2) The majority of these studies require user interaction to select a ROI (a small window) around the prostate.Unfortunately, such approaches not only prone to interobserver variability, but also ROI selection biases the final decision by over-or under-estimating the problem in the entire graft, just as with biopsy.
3) Automated prostate segmentation methods have one of the following limitations:  Deformable model-based methods without adequate appearance and shape priors fail under excessive noise, poor resolution, diffused boundaries, or oc-cluded shapes in the images;  Segmentation based only on the shape prior still results in large errors caused by discontinuities in object boundaries, large image noise, and other inhomogeneities;  Parametric shape-based models are unsuitable for discontinuous prostate objects due to a very small number of distinct landmarks.4) The majority of CADs assumes that the prostates (prostate contours) remain exactly the same from scan to scan.However, prostate contours may not always exactly match due to patient movement or breathing effects; therefore, image registration schemes should be applied first before ROI selection/segmentation.
To overcome these limitations, we propose an automatic framework for analyzing DWI images building on our previous work in [46][47][48].The proposed approach consists of the following steps as shown in Figure 1: 1) Segmentation of the prostate from DWI (Section 2.1) based on a Maximum a Posteriori (MAP) estimate of a new likelihood function that accounts for both appearance features of the prostate (Section 2.1.1)and their 3D spatial voxel interactions (Section 2.1.2),as well as a 3D shape prior (Section 2.1.3).
2) A non-rigid registration approach is employed to account for any local deformation that could occur in the prostate during the scanning process based on the solution of the Laplace equation (Section 2.3).
3) KNN classifier to classify the prostate into benign or malignant based on three appearance features extracted from registered images (Section 3.2).

MATERIALS AND METHODS
In this paper we introduce a new, automated, and noninvasive framework for early diagnosis of prostate cancer from DWI. Figure 1 demonstrates the steps of the proposed CAD system.Below, we will illustrate each of these steps.

Segmentation of the Prostate Using a Joint MGRF Model
The segmentation of the prostate is a challenge, since the gray-level distribution of the prostate and surrounding organs is not highly distinguishable and because of the anatomical complexity of prostate.This stage proposes a powerful framework for prostate segmentation based on a learned shape model and an identifiable joint Markov-Gibbs Random Field (MGRF) model of DWI and "object-background" region maps.
The joint-MGRF model is fundamentally a model that relates the joint probability of an image and its objectbackground region map, to geometric structure and to the nergy of repeated patterns within the image.The basic e theory behind such models is that they assume that the signals associated with each pixel depend on the signals of the neighboring pixels, and thus explicitly take into account their spatial interactions, and other features, such as the shape.
Let , , and integer gray-level, a set of object ("ob") and background ("bg") labels, and a unit interval, respectively.Let a 3D arithmetic grid . The LCDG modeling restores transitions between these components more accurately than conventional mixtures of only positive Gaussians, thus yielding a better initial region map formed by voxel-wise classification of the image gray values, the similar intensity profile of the prostate and surrounding tissues.tions between the region labels are also taken into account using the popular Potts model, i.e., the MGRF with the nearest voxel 26-neighborhood (see Figure 3).
A generic MGRF of region maps accounts only for pairwise interaction between each region label and its characteristic neighbors.Generally, the interaction structure and the Gibbs potentials can be arbitrary and are identified from the training data.
By symmetry considerations, we assume that the potentials are independent of relative orientation of each voxel pair and depend only on intra-or inter-region position (i.e.whether the labels are equal or not).Under these restrictions, it is the 3D extension of the conventional auto-binomial, or Potts model differing only in that the potentials are estimated analytically.
The 26-neighborhood has three types of symmetric pairwise interactions specified by the absolute distance a between two voxels in the same and adjacent MRI slices ( , 1 a  2 , and 3 , respectively): 1) the closest pairs with the inter-voxel coordinate offsets; 2) the diagonal pairs with the offsets ; and 3) the farthest diagonal pairs with the offsets    

3
. The Gibbs potentials of each type are bi-valued because only label coincidence is accounted for: Then the MGRF model of region maps is as follows [52,53]: ,  , , , , 1 exp , where Z is the normalizing factor (partition function).
To identify the MGRF in Eq.1, approximate analytical maximum likelihood estimate of the 3D Gibbs potentials, and are derived in line with [52]: a eq a ne a eq where,

 
, a eq f m denotes the relative frequency of the equal labels in the equivalent voxel pairs of a region map of a given DWI aligned in accord with the prior shape model.m

Probabilistic Shape Model
To enhance the segmentation accuracy, the expected shape of the goal object is constrained with a soft probabilistic 3D prostate shape model.Initially, a training database collected from different subjects are co-aligned by rigid, affine 3-D transformations.The shape prior is a spatially variant independent random field of region labels: where is the empirical probability that the voxel (x, y, z) belongs to the prostate (L = "ob") or the background (L = "bg") given the map.To enhance the segmentation of the current prostate volume, the prior probabilistic shape model is updated by adding the previous segmented 3D prostate data to the prior calculated shape model.The proposed prostate segmentation process can be summarize as follows:  Perform an affine alignment of a given 3D MRI to an arbitrary prototype prostate from the training set using mutual information [54] as a similarity measure to obtain the learned probabilistic shape model

Performance Evaluation of the Proposed Segmentation Algorithm
The proposed segmentation is evaluated based on char-acterizing the agreement (Figure 4

Nonrigid Registration
Due to patient breathing and local movement, accurate registration is a main issue in DWI.In this paper, the nonrigid motion of the DWI data at different b-values is compensated for by using our developed registration approach that is based on the solution of the second-order partial differential Laplace equation [56]: for a scalar function   , x y  B between the target and the reference prostate objects.The solution of a planar Laplace equation between two boundaries results in intermediate equipotential surfaces (dashed lines in Figure 5) and streamlines that establish natural point-to-point correspondences and are everywhere orthogonal to all the equipotential surfaces (see e.g., the line connecting the points ai and ai in Figure 5).Based on solving the Laplace equation, we perform the non-rigid registration as follows: 2) Use these distance maps to generate equispaced isocontours as shown in Figures 6(c) and (d).
3) Solve the Laplace equation between respective reference and target iso-contours to find the point-to-point correspondence.

Color Map Generation and Tumor Boundary Determination
To characterize the physiological data, color-coded maps that illustrate the propagation of diffusion in the prostate tissues are constructed.To construct the initial color maps, we have to estimate the changes in image signals , , x y z  due to the Brownian motion.These changes are estimated from the constructed normalized diffusion as the difference between the signals of image sequences at two different b-values.DWI is performed with at least two b values, including a b value of 0 sec/mm 2 and a higher b value of 500 -1000 s/mm 2 depending on the body region or organ being imaged [57].At b = 0 s/mm 2 , there is no diffusion sensitizing gradient with free water molecules have high signal intensity.We used b = 800 s/mm 2 because imaging of solid organs requires high b value specially in prostate and using high b values allows differentiation of areas of restricted from the normal high signal at the peripheral zone.During our trials we found the b = 800 s/mm 2 allows lesions differentiation with least degradation of image quality as the image quality decrease with the high b values.
To preserve continuity (remove inconsistencies), the initial estimated , , x y z  values are considered as samples from a Generalized Gauss-Markov Random Field (GGMRF) image model [58] of measurements with the 26-voxel neighborhood (Figure 3).Continuity of the constructed 3-D volume (Figure 7) is amplified by using their MAP estimates [51]: Finally, to allocate the boundary of the detected tumors, which is important to determine the cancer stage in case of malignancy, we used a level set-based deformable model controlled by a stochastic speed function [59].The latter accounts for the perfusion information and spatial interactions between the prostate voxels.

EXPERIMENTAL RESULTS
The performance of the proposed framework has been evaluated by applying it on DWI prostate images collected from 28 patients.These patients had biopsyproven prostate cancer and underwent DWI at 1.5-T (SIGNA Horizon, General Electric Medical Systems, Milwaukee, WIS).DWI were then obtained using monodirectional gradients and a multi-section Fast Spin Echo type (FSE) echo-planar sequence in the axial plane using a body coil with the following imaging parameters: TE: 84:6 ms; TR: 8.000 ms; Band Width 142 kHz; FOV 34 cm; slice thickness 3 mm; inter-slice gap 0 mm; seven excitations, water excitation with b value of 0 s/mm 2 and 800 s/mm 2 .Fifty four slices were obtained in 120 second.to cover the prostate in each patient.Note all the subjects were diagnosed using a biopsy (ground truth).

Segmentation Results
The proposed segmentation approach has been tested on 28 independent data sets of DWI images.Figure 9 shows sample examples of prostate segmentation from different data sets, with respect to the ground truth segmentation.The ground truths were obtained by manual delineation of the prostate borders by an MR imaging expert.To highlight the advantages of our segmentation technique, we compare it to the shape-based segmentation approach proposed by Tsai et al. [60].We re-implemented the method described in [60] and tested it on our locallyacquired data.Figure 9 compares qualitatively the accuracy of our approach and the shape-based approach [60] with respect to the ground truth.The segmentation accuracy for all data sets has been evaluated using the average segmentation error, given by Eq.7.Differences between the mean errors for our segmentation and the shape-based approach [60] in Table 1 are statistically significant by the unpaired t-test and thus highlight the advantages of the proposed integration of the shape prior, prostate/background marginal intensity distributions, and spatial interaction characteristics into MAP-based segmentation.
Moreover, the accuracy of our segmentation approach has been evaluated, with respect to the expert tracing, using the PPV, Sens, DSC [61], and the APD between the borders of ground truth G and automatic segmentation C (see Figure 10).Table 2 compares the segmentation over all the test data sets with the truth obtained by manual tracing by an expert.

Diagnostic Results
The ultimate goal of the proposed framework is to distinguish between benign and malignant detected tumors.The malignant tissues show higher signal intensity with a b-value of 800 s/mm 2 , and a lower Apparent Diffusion Coefficient (ADC) compared with benign and normal tissue due to the replacement of normal tissue.To distinguish between the benign and malignant cases, we used a KNN classifier learning statistical characteristics   of the DWI.The characteristics are obtained from the training sets containing both benign and malignant cases.
After training, three features namely are the mean intensity value of the DWI at 0 s/mm 2 , the mean intensity value of the DWI at 800 s/mm 2 , and the mean value of ADC maps [62], were chosen to classify the test cases.
To build the KNN classifier that characterizes the prostate tissue, we used 13 subjects for training, and the other 15 subjects for testing.The diagnostic accuracy based on the combined three features resulted in correct classifications of all 28 data sets (i.e., 100% accuracy).Tumor's contour determination (green) using the level set approach for multiple image sections for benign (B) and malignant (M) subjects.
For regional display we explore pixel-by-pixel maps of the registered diffusion data.The diffusion was computed for each pixel and superimposed on an image slice to form a parametric image.Also, for visual assessment of the prostate tumor the tumor contours are determined.Figure 11 shows the tumor contours determination for selected image sections for two subjects involved in our study.

CONCLUSIONS AND DISCUSSION
In this paper, we present a novel fully automatic framework for detecting prostate cancer using DWI.The framework includes prostate segmentation, nonrigid registration, and KNN-based classification.For prostate segmentation, we introduced a new 3D approach that is based on a MAP estimate of a new log-likelihood function that accounts for the shape priori, the spatial interaction, and the current appearance of the prostate and its background which increases the accuracy of automatic segmentation, evidenced by the error and the DSC analysis (Tables 1 and 2).Following segmentation, we used a nonrigid registration approach that deforms the prostate object on iso-contours instead of a square lattice, which provides higher degrees of freedom to obtain accurate deformation.In the classification step, the segmented prostate regions are classified into malignant or benign using the KNN classifier.Applications of the proposed framework can assist the radiologist in detecting all prostate cancer locations and could replace the use of current technologies to determine the type of prostate cancer.
Although we have obtained promising results in this initial study using DWI data in 28 patients, potential widespread adoption would require confirmation by other groups, and investigation in a larger number of subjects.Our future work will focus on comparing the diagnostic accuracy of prostate cancer detection using other imaging modalities, such as DCE-MRI.

Figure 1 .
Figure 1.Flowchart of the proposed CAD system for automatic detection of cancer from 3D DWI.

Figure 2 .Pm
Figure 2. Aligning a 3-D joint Markov-Gibbs random field model with shape prior of DWI.

2 . 1 . 2 .
In this work we focus on accurate identification of the spatial interaction between the prostate voxels   P m , and the intensity distribution for the prostate tissues,  P  g m , and the prior distribution of the prostate shape, as shown in Figure 2. Spatial Voxel Interaction Model In order to overcome noise effect and to ensure the homogeneity of the segmentation, spatially voxel interac-To perform the initial prostate segmentation, a given 3D DWI is aligned to one of the training 3D DWI.The OPEN ACCESS

Figure 3 .
Figure 3. Pairwise voxel interaction for 26 neighborhood system in a 3D GGMRF.The reference voxel is shown in red.


Estimate the conditional intensity model   P g m by identifying the bimodal LCDG;  Use the intensity model found and the learned probabilistic shape model to perform an initial segmentation of the prostate, i.e., to form an initial region map;  Use the initial region map to estimate the potential for the Potts model and to identify the MGRF model   P m of region maps;  Improve the region map using voxel-wise stochastic relaxation (Iterative Conditional Mode (ICM) [55]) through successive iterations to maximize the log likelihood function of Eq.1 until the log likelihood remains almost the same for two successive iterations;  Output: The 3D prostate segmentation is the final estimate m .
(a)) and the Average Perpendicular Distance (APD) between the segmented and ground truth contours (Figure 4(b)).To evaluate the performance, we measured True Positive (TP), True Negative (TN), False Positive (FP), and False Negative (FN) segmentation (Figure 4(a)).Let C and G denote the segmented region and its "ground truth" counterpart, respectively.Let z denote the volume (in the number of voxels) of a region z.Then, TP C G   is the overlap between C and G, FP C C G    is the difference between C and TP; and FN G C G    is the difference between G and TP; and TN R C G    .The Positive Predictive Value (PPV), Sensitivity (Sens), Dice Similarity Coefficient (DSC), and the average segmentation error are defined as:

B 1 )
Generate the distance maps inside the prostate regions as shown in Figures 6(a) and (b).

Figure 4 .
Figure 4. 2-D schematic illustration of measuring segmentation errors (a) and (b) perpendicular distances (black lines) (b) between the ground truth G and automatic segmentation C.

Figure 5 .Figure 6 .
Figure 5. 2-D illustration of co-allocation of point-to-point correspondences between two borders by a potential field.

 2 
is the GGMRF potential, and  and  are scaling factors.The paramecontrols the level of smoothing (e.g., smooth,   , vs. relatively abrupt edges,  , prior distribution of the estimator.Then, the color maps are generated based on the final estimated   (see Figure8).

Figure 9 .
Figure 9. 3-D prostate segmentation projected onto 2-D.(a) Our segmentation (red) in comparison with the ground truth (white); (b) The segmentation with the algorithm in [60] (red) comparison with the ground truth; and (c) 2-D visualization for our segmented prostates for three of the test subjects.

Figure 10 .
Figure 10.Prostate image with ground truth (blue) and automatic segmentation (green) contours, and their associated streamlines (red) obtained by the solution of the Laplace equation yielding the estimation of the APD.

Figure 11 .
Figure 11.Tumor's contour determination (green) using the level set approach for multiple image sections for benign (B) and malignant (M) subjects.

Table 1 .
[60]mparative segmentation accuracy over all test data sets for our approach and[60].Note that "STD" stands for standard deviation.

Table 2 .
Error statistics over all test data sets.Note that "STD" stands for standard deviation and "APD" values are in mm.