An Efficient Convex Optimization Approach to 3D Prostate MRI Segmentation with Generic Star Shape Prior Jing Yuan, Wu Qiu, Eranga Ukwatta, Martin Rajchl, Yue Sun, and Aaron Fenster Medical Imaging Laboratory, Robarts Research Institute, Western University, London, ON, N6A 5K8 Canada

Abstract. In this work, we propose a novel global optimization-based contour evolution approach for the segmentation of 3D prostate magnetic resonance (MR) images, which incorporates histogram-matching and a novel variational formulation of a generic star shape prior. Our method overcomes the existing challenges of segmenting 3D prostate MRIs: heterogeneous intensity distributions and a wide variety of prostate shape appearances. We introduce a novel convex relaxation-based method to evolve a contour to its globally optimal position during each discrete time frame. The proposed generic star shape prior provides robustness to the segmentation when the image suffer from poor quality, noise, and artifacts. Our approach provides a fully time implicit scheme to contour evolution, which allows a large time step-size to accelerate the speed of convergence. Moreover, a new continuous max-flow formulation is proposed which is dual to the convex relaxation formulation, obtains global optimality of contour evolution, and is implemented in a GPU.

1

Background

Prostate cancer is the most common non-skin cancer in men of the developed countries. It is also the third leading cause of death due to cancer [1]. MR guided focal therapies such as radio-therapy, laser ablation are increasingly used for prostate cancer treatment due to the excellent soft-tissue contrast provided by MR imaging [2]. Accurate and precise segmentation of the prostate boundary from 3D MR images is an essential and crucial step of radio-therapy planning to spare the neighbouring tissues such as bladder, rectal wall and seminal vesicles, from exposing to any radiation. In addition, the segmented prostate boundary is also required for the surface-based MR to 3DUS registration in 3D ultrasound (3DUS) guided targeted biopsies, where a surface based registration is preferred over an image-based registration due to the deformations caused by the endorectal coil in MR. Even though there are intensive studies in prostate segmentation, specifically in the segmentation of transrectal ultrasound (TRUS), the segmentation of the prostate from in vivo endorectal T2 weighted 3D MR images needs specifically

82

developed techniques due to its existing challenges: first, sharp imaging edges of sub-regions are widely distributed in scanned MRIs, i.e. the existing high inhomogeneity of intensities, which seriously bias finding the correct segmentation boundaries; second, the wide variety of appearing prostate shapes makes it often a computationally intensive task to learn the precise shape information, especially for 3D prostate shapes, so as to successfully guide a shape-based segmentation task [3,4]; moreover, dealing with variations of prostate shape and size across different patient studies also needs additional computation cost for segmenting the correct prostate boundaries. In this work, we propose a new global optimization based contour evolution approach to segmenting 3D prostate MRIs efficiently and robustly, which addresses the associated challenges by matching histograms in- and outside the object [5,6] and incorporating the generic star-shape prior [7,8] (see Fig. 1(a) for illustration). By means of convex relaxation, we show the given contour is gradually moved to its globally optimal position during each discrete time frame subject to the introduced star shape constraint, which is .essentially distinct from previous local optimization based methods, e.g. active contours or levelsets [9,10,11] etc, due to its new theoretical perspective of global optimization. In this regard, our approach introduces a fully time implicit contour evolution scheme, which allows a large time step-size to significantly speed up contour evolution and improves numerical reliabilities. The new definition of star shapes makes it possible to be implicitly embedded in the proposed optimization algorithm and solved globally with a low extra computation cost. In numerics, a new continuous max-flow based algorithm is proposed and implemented on GPU to achieve high computational performance.

2

Methodology

(a)

(b)

(c)

Fig. 1. 1(a) gives examples of the star shapes w.r.t. the central point. 1(b) shows the variational definition of the star shape prior to handle the variety of prostate shape appearances. 1(c) illustrates the evolution of contours.

83

In this paper, we introduce a global optimization based contour evolution approach to 3D prostate MRI segmentation and use histogram matching and the star-shape prior to address the two challenges of the segmentation task. 2.1

Optimization Formulation

Histogram Matching The probability density functions (PDFs) of the intensities inside and outside the object region of interests provide a global and practical description of the object to be segmented [5,6], which helps to significantly reduce the bias of segmented boundaries by the highly inhomogeneous intensity distributions. The model PDFs can be learned from either the training datasets or the sampled pixels. Let I(x) be the given 3D prostate MRI volume and u(x) ∈ {0, 1} be the indicator function of the estimated prostate region C. Given the intensity distribution models inside and outside the prostate region, i.e. qin (z) and qout (z) where z ∈ Z gives the photometric value of intensities, the Bhattacharyya distance [6] is used for PDFs’ matching of both inside and outside regions, which results in the following model-matching term:     Ematching (u) = − pin (u, z) qin (z) − pout (1 − u, z) qout (z) (1) z∈Z

z∈Z

where pin (u, z) and pout (u, z) are the respective PDFs for the estimated inside and outside regions of prostate, and computed by the Parzen method [12]:   K(z − I(x)) u dx K(z − I(x)) (1 − u) dx Ω   , pout (1 − u, z) = Ω pin (u, z) = u dx (1 − u) dx 1 exp(−x2 /2σ 2 ). where K(·) is the Gaussian kernel function K(x) = √2πσ 2 Star Shape Prior A star shape was first proposed in [7] by graph-cuts, which is defined with respect to the center point o (see Fig. 1(b)): an object has a star shape if for any point x inside the object, all points on the straight line between the center o and x also lie inside the object; in another word, the object boundary can only pass any radial line starting from the origin o one single time. In this work, we introduce a new variational formulation of the star-shape prior in the spatially continuous setting: for the center point o, let do (x) be the distance map with respect to o and e(x) = ∇dO (x); the star shape prior can be defined as follows: ∇u(x) · e(x) ≥ 0 , ∀x ∈ Ω . (2)

The star shape prior (2) models a subset of the simply connected regions, where the holes and disjoint regions are exactly ruled out (illustrated by Fig. 1(b)). Optimization Model In view of the histogram matching energy function (1) and the star shape constraint (2), we propose to segment 3D prostate MRI by minimizing the following energy function  Ematching (u) + g(x) |∇u(x)| dx , s.t. ∇u(x) · e(x) ≥ 0 . (3) min u(x)∈{0,1}

Ω

where the anisotropic total-variation term gives the boundary smoothness term.

84

2.2

Global Optimization Based Contour Evolution

We solve the challenging non-convex optimization function (3) by a new global optimization based contour evolution approach, which can efficiently and reliably evolve the contour to its globally optimal position at each discrete time frame: for the given contour Ct at time t, we propagate Ct to its globally optimal position Ct+1 at the next time frame, in terms of minimizing (3). The associated theory is also introduced in [13]. For the global optimization based contour evolution, we introduce the costs c+ (x) and c− (x) w.r.t. two distinct regions C + and C − (see Fig. 1(c)): 1. C + indicates the expanded region w.r.t. Ct : for ∀x ∈ C + , it is initially outside Ct at time t, and ’jumps’ to be inside C at t + 1; for such a ’jump’, it pays the cost c+ (x). 2. C − indicates the shrunk region w.r.t. Ct : for ∀x ∈ C − , it is initially inside Ct at t, and ‘jumps’ to be outside C at t + 1; for such a ’jump’, it pays the cost c− (x). We propose to propagate C t to C by solving the following optimization problem globally and exactly:    + − c (x) dx + c (x) dx + g(s) ds . (4) min C

C+

C−

∂C

In fact, (4) addresses the contour evolution by achieving the minimum total costs paid by the two regions C + and C − , where the label of pixels change: either from 1 to 0, i.e. x ∈ C − , or from 0 to 1, i.e. x ∈ C + . We can apply (4) to solve the proposed optimization problem (3) such that the two cost functions c+ (x) and c− (x) are set up to be the first-order derivatives w.r.t. the two terms of (1) respectively (see [6] for details), along with the distance functions w.r.t. Ct . Moreover, we define the cost functions Cs (x) and Ct (x):  −  + c (x) , where x ∈ C c (x) , where x ∈ /C Cs (x) := , Ct (x) := . (5) 0, otherwise 0, otherwise By (5), the optimization problem (4) can be equally reformulated as  min

u(x)∈{0,1}

1 − u, Cs  + u, Ct  +

Ω

g(x) |∇u| dx , s.t. ∇u(x) · e(x) ≥ 0 . (6)

Convex Relaxation and New Continuous Max-Flow Approach Clearly, (6) gives rise to a star-shape constrained continuous min-cut model. We show it can be solved globally and exactly by its convex relaxation, i.e.  g(x) |∇u| dx , s.t. ∇u(x) · e(x) ≥ 0 ; (7) min 1 − u, Cs  + u, Ct  + u(x)∈[0,1]

Ω

85

under its convex duality perspective. In this regard, we propose the novel continuous max-flow formulation based on a similar flow configuration as [14], for which an additional spatial flow q(x) = λ(x)e(x) and λ(x) ≥ 0 is introduced. The new continuous max-flow model is formulated as follows:  ps (x) dx (8) min ps ,pt ,p,λ

Ω

subject to |p(x)| ≤ g(x) , ps (x) ≤ Cs (x) , pt (x) ≤ Ct (x) ,  

div p(x) + λ(x)e(x) − ps (x) + pt (x) = 0,

∀x ∈ Ω ;

(9)

∀x ∈ Ω ;

(10)

along with λ(x) ≥ 0. We can prove the introduced continuous max-flow model (8) is dual to the convex relaxation problem (7), by means of similar analytical steps in [14] together with the variational analysis over λ(x) ≥ 0. With helps of the proposed continuous max-flow model (8), we can also prove Proposition 1 The non-convex combinatorial optimization problem (6) can be solved by its convex relaxation (7) exactly and globally, such that the optimum u∗ (x) ∈ [0, 1] of (7) can be thresholded by any given constant α ∈ [0, 1) to be u# (x) ∈ {0, 1} which solves (6) globally. Its proof follows the similar steps in [14] and omitted here due to limit space. A big advantage to utilize the proposed continuous max-flow model (8) is that it results in a numerically simple and efficient algorithm to its dual formulation (7), where the nonlinear total-variation function and the associated star shape constraints are encoded by simple projections to the corresponding convex sets. The continuous max-flow algorithm is similar to the one proposed in [14], beside an extra flow adjustment step at each iteration for λ which can be implemented by an additional projected-gradient step, also see [13] for details. Clearly, the contour evolution approach proposed in this paper is essentially different from previous methods, e.g. level-sets, which explicitly solve the evolution PDE with local optimization based schemes. It does introduce a fully time implicit method to the contour evolution, where the time step-size is implicitly adapted into the regularization weight function g(x) and allows a large value to speed up contour evolution in practice.

3

Experimental Design

We summarize our experiment settings as follows: – The approach presented in this study is initialized by a closed surface, which is constructed by the thin-plate spline with positioning eight to ten initial points on the boundary of the prostate [six at the transverse view (red points in Fig. 2], two or four at the sagittal view (green points in Fig. 2)). Such smooth closed surface approximates the prostate shape and provides a relatively better initialization condition (see [15] for more details); such that the inside and outside

86

voxels of the estimated surface are used to generate prior PDFs for the foreground and background and the closed surface also defines the initial value of the indicator function u(x) for the proposed algorithm. – The segmentation results are compared to the given manual segmentation using volume-based metrics: sensitivity(Se), DSC (dice similarity coefficient), and distance-based metrics: the mean absolute surface distance(MAD) and maximum absolute surface distance(MAXD)[16]. – Five images from the training datasets are used to optimize the parameters of the algorithm, and the rest are used to validate the proposed method. The parameter values are empirically chosen. Afterwards, the parameters are optimized sequentially by changing a single parameter at a time while holding other parameters fixed. Then the optimized values of parameter used in our experiments are fixed during the validation experiments.

Fig. 2. The scheme of initialization (six red points on the transverse view and two green points on the sagittal view).

4 4.1

Results and Discussion Qualitative results

Table 1. Overall performance results DSC(%) 86.2 ± 3.0

Se(%) 85.4 ± 5.3

MAD(vx) 6.38 ± 4.2

MAXD(vx) 15.7 ± 6.3

The validation results in Table.1 show that our method can obtain a DSC of 86.2 ± 3.0%, a sensitivity of 85.4 ± 5.3%, a M AD of 6.38 ± 4.2 voxels, and M AXD of 15.7±6.3 voxels. To validate the variability introduced by the manual initialization, an intra-observer variability test was performed. All images were segmented three times. The coefficient of variation(CV ) [17] of DSC was used

87

to evaluate the intra-observer variability of our method. The result shows that our proposed method yielded a CV of 3.8%. 4.2 Efficiency The proposed convex max-flow algorithm was implemented using parallel computing architecture (CUDA, NVIDIA Corp., Santa Clara, CA) and the user interface in Matlab (Natick, MA). The experiments were conducted on a Windows desktop with a Intel Core i7-2600 CPU of 3.4 GHz and a GPU of NVDIA GTX580. More detailed information about the runtime environment is listed in Table 2. The initialization for our algorithm takes about 30 seconds using the software developed by our laboratory. The computational time for convergence is about 15 seconds in average, which is achieved within only 5-10 outer iterations for a single 3D MR image. It significantly outperforms over the level-set based methods for which hundreds of iterations will be applied to achieve convergence! Table 2. Detailed information about the runtime environment.

Algorithm

Parameter Value Language: Matlab 7.12.0 Libraries/Packages: CUDA Thrust GPU Optimizations: Yes Multi-Threaded: No

Machine

CPU Clock Speed: 3.4 GHz

5

Machine CPU Count: 1 Machine Memory: 8 GB Memory Used During Segmentation: 2 GB

Concluding Remarks

We describe and evaluate a novel global optimization based contour evolution approach to 3D prostate MRI segmentation with the integration of the generic star shape prior, which properly addresses the existing challenges of prostate MRI segmentation and achieves numerical efficiency and reliabilities. In practice, the high anisotropy of the sampled 3D prostate MRI data does interfere achieving high segmentation accuracy. A slice wise segmentation may be an option to improve our segmentation results in the future, for which the proposed algorithm can be directly adapted to such 2D cases. In addition, the learned image texture information provides another robust imaging clue to prostate segmentation.

Acknowledgments The authors are grateful for the funding support from the Canadian Institutes of Health Research (CIHR), the Ontario Institute of Cancer Research (OICR), and the Canada Research Chairs (CRC) Program for this work.

88

References 1. Jemal, A., Bray, F., Center, M.M., Ferlay, J., Ward, E., Forman, D.: Global cancer statistics. CA: A Cancer Journal for Clinicians 61(2) (March 2011) 69–90 2. Xu, S., Kruecker, J., Turkbey, B., Glossop, N., Singh, A.K., Choyke, P., Pinto, P., Wood, B.J.: Real-time MRI-TRUS fusion for guidance of targeted prostate biopsies. Computer aided surgery : official journal of the International Society for Computer Aided Surgery 13(5) (September 2008) 255–264 PMID: 18821344 PMCID: 2664902. 3. Heimann, T., Meinzer, H.P.: Statistical shape models for 3d medical image segmentation: A review. Medical Image Analysis 13(4) (2009) 543 – 563 4. Toth, R., Tiwari, P., Rosen, M., Reed, G., Kurhanewicz, J., Kalyanpur, A., Pungavkar, S., Madabhushi, A.: A magnetic resonance spectroscopy driven initialization scheme for active shape model based prostate segmentation. Medical Image Analysis 15(2) (2011) 214225 5. Aubert, G., Barlaud, M., Faugeras, O., Jehan-Besson, S.: Image segmentation using active contours: calculus of variations or shape gradients? SIAM J. Appl. Math. 63(6) (2003) 2128–2154 6. Michailovich, O., Rathi, Y., Tannenbaum, A.: Image segmentation using active contours driven by the bhattacharyya gradient flow. IEEE Transactions on Image Processing 16(11) (2007) 2787–2801 7. Veksler, O.: Star shape prior for graph-cut image segmentation. In Forsyth, D., Torr, P., Zisserman, A., eds.: ECCV 2008. Volume 5304 of Lecture Notes in Computer Science. (2008) 454–467 8. Strekalovskiy, E., Cremers, D.: Generalized ordering constraints for multilabel optimization. In: Computer Vision (ICCV), 2011 IEEE International Conference on. (nov. 2011) 2619 –2626 9. Caselles, V., Kimmel, R., Sapiro, G.: Geodesic active contours. International Journal of Computer Vision 22(1) (1997) 61–79 10. Sundaramoorthi, G., Yezzi, A., Mennucci, A., Sapiro, G.: New possibilities with sobolev active contours. In: SSVM. (2007) 153–164 11. Chan, T., Vese, L.A.: Active contours without edges. IEEE Image Proc., 10, pp. 266-277 (2001) 12. Parzen, E.: On estimation of a probability density function and mode. The Annals of Mathematical Statistics 33(3) (1962) pp. 1065–1076 13. Yuan, J., Ukwatta, E., Tai, X.C., Fenster, A., Schnoerr, C.: A fast global optimization-based approach to evolving contours with generic shape prior. Technical report CAM-12-38, in submission, UCLA (2012) 14. Yuan, J., Bae, E., Tai, X.: A study on continuous max-flow and min-cut approaches. In: CVPR 2010 15. Hu, N., Downey, D.B., Fenster, A., Ladak, H.M.: Prostate boundary segmentation from 3d ultrasound images. Med. Phys. 30(7) (Jul 2003) 1648–1659 16. Garnier, C., Bellanger, J.J., Wu, K., Shu, H., Costet, N., Mathieu, R., de Crevoisier, R., Coatrieux, J.L.: Prostate segmentation in hifu therapy. IEEE Trans. Med. Imag. 30(3) (2011) 792–803 17. Zou, K.H., Mcdermott, M.P.: Higher-moment approaches to approximate interval estimation for a certain intraclass correlation coefficient. Statistics in Medicine 18(15) (1999) 2051–2061

89

An Efficient Convex Optimization Approach to 3D ...

evolution, which allows a large time step-size to accelerate the speed of .... (8) is dual to the convex relaxation problem (7), by means of similar analytical .... the high anisotropy of the sampled 3D prostate MRI data does interfere achieving.

137KB Sizes 2 Downloads 223 Views

Recommend Documents

LNCS 7510 - A Fast Convex Optimization Approach to ...
scar tissue solely from 3D cardiac delayed-enhancement MR images (DE-. MRI). ..... Foundation for Innovation (CFI) grant #20994, Canadian Institutes of Health.

Convex Optimization
Mar 15, 1999 - 5.1 Unconstrained minimization and extensions . ..... It is (as the name implies) a convex cone. Example. ..... and lies in the domain of f (i.e., c. T.

Efficient Global Optimization Based 3D Carotid AB-LIB ...
London, ON, Canada ... black blood MR images, by simultaneously evolving two coupled surfaces ... flow, image segmentation, GPGPU, coupled level sets.

an efficient signal-matching approach to melody ...
audio (average 3.6 minutes per song). 155 of the songs .... An integrated pitch tracking algorithm for speech systems. In IEEE Int. Conf. Acoust., Speech.

Guaranteed Non-convex Optimization: Submodular ...
Submodular continuous functions naturally find applications in various real-world settings, including influence and revenue maximization with continuous assign- ments, sensor energy management, multi-resolution data summarization, facility location,

A Study on Convex Optimization Approaches to Image ...
ful applications include image denoising [22,8,30], image decomposition [2,17] ..... + ui/c)}) . (30). Interestingly, our experiments show that just one single step of ... experiments are computed by a Ubuntu desktop with AMD Athalon 64 X2 5600.

boyd convex optimization pdf
Download. Connect more apps... Try one of the apps below to open or edit this item. boyd convex optimization pdf. boyd convex optimization pdf. Open. Extract.

Particle Swarm Optimization: An Efficient Method for Tracing Periodic ...
[email protected] e [email protected] ..... http://www.adaptiveview.com/articles/ipsop1.html, 2003. [10] J. F. Schutte ... email:[email protected]

3D Object Retrieval using an Efficient and Compact ...
[Information Storage and Retrieval]: Information Search and Retrieval. 1. ... shape descriptor that provides top discriminative power and .... plexity which is prohibitive for online interaction as well as .... degree of significance, we can intuitiv

Conscience online learning: an efficient approach for ... - Springer Link
May 24, 2011 - as computer science, medical science, social science, and economics ...... ics in 2008 and M.Sc. degree in computer science in 2010 from Sun.