Prostate Volume Segmentation in TRUS Using Hybrid Edge-Bhattacharyya Active Surfaces

Objective: We present a new hybrid edge and region-based parametric deformable model, or active surface, for prostate volume segmentation in transrectal ultrasound (TRUS) images. Methods: Our contribution is threefold. First, we develop a new edge detector derived from the radial bas-relief approach, allowing for better scalar prostate edge detection in low contrast configurations. Second, we combine an edge-based force derived from the proposed edge detector with a new region-based force driven by the Bhattacharyya gradient flow and adapted to the case of parametric active surfaces. Finally, we develop a quasi-automatic initialization technique for deformable models by analyzing the profiles of the proposed edge detector response radially to obtain initial landmark points toward which an initial surface model is warped. Results: We validate our method on a set of 36 TRUS images for which manual delineations were performed by two expert radiation oncologists, using a wide variety of quantitative metrics. The proposed hybrid model achieved state-of-the-art segmentation accuracy. Conclusion: Results demonstrate the interest of the proposed hybrid framework for accurate prostate volume segmentation. Significance: This paper presents a modular framework for accurate prostate volume segmentation in TRUS, broadening the range of available strategies to tackle this open problem.


A. Context
Prostate cancer is the most frequent cancer and the third cause of cancer-related mortality in men of 50 years of age or above worldwide.In France, Germany and Switzerland, incident rates have been consistently increasing since 1990 by 4-5% each year [1].Among available treatments, low-dose rate permanent prostate brachytherapy (LDR-B) has emerged as one of the most effective methods for treating localized cancers.LDR-B consists in delivering radiation directly within the gland, by permanently inserting grain-sized radioactive seeds inside the prostate.A typical procedure requires dose planning prior to the operation or alternatively just before Submitted on February 28 th , 2018.This work was partly supported by the French ANR within the Investissements d'Avenir program (Labex CAMI) under reference ANR-11-LABX-0004 (Integrated project CAPRI).
* Asterisk indicates corresponding author (correspondence e-mail: vjaouen@univ-brest.fr)implantation to determine the accurate number and positioning of brachytherapy seeds adapted to the patient.Batch of seeds are inserted one after another using needles across the perineum according to the planning.In a clinical setup, the insertion procedure is supervised under continuous transrectal ultrasound (TRUS) guidance (Fig. 1) [2], [3].
Prostate volume segmentation (PVS) is a critical task of the brachytherapy procedure.It is for example necessary during dose planning before the operation, alongside the delineation of neighboring Organs at Risk (OAR) urethra and rectum, to calculate a dose plan in order to optimize the dose delivered to the tumor while minimizing adverse effects to the OARs.
Manual delineation from TRUS images is the most common practice in clinical routine during surgery.It is however particularly tedious in 3D, especially when repeated several times along the operation, increasing the total duration of the surgery in unacceptable proportions.Due to the low quality of TRUS images, alternative image modalities such as CT or MRI can be considered to perform PVS outside the operating room, at the expense of superior time and cost.The American Society for Brachytherapy also objects that "TRUS potentially offers the only practical option for performing on-line dosimetric analysis during the procedure" [4].There is therefore a clear need for automatic or semi-automatic TRUS-based PVS methods that could perform this task with satisfying accuracy, repeatability and speed

B. TRUS-based prostate volume segmentation challenges
The prostate gland globally appears in B-mode TRUS images as a hypoechoic (dark) walnut-shaped region at the center of the field of view (FOV), surrounded by more echoic (bright) medium.However, its accurate segmentation is made difficult by the low quality of the images obtained.Two major difficulties can be identified.The first main difficulty lies in accurate prostate edge detection.Edge signal in TRUS is of low to very low quality, depending on the subject and the image region considered.Fig. 2a shows mid-plane transverse slices for 12 different patients who underwent LDR-B.In these images, variations in image intensity are not only caused by the prostate interface, but also by various undesired artifacts such as hyperechoic, patient dependent, calcifications.Echoic shadows along the beamline, located behind hyperechoic regions, create spurious strong edges, the most consistent of which being caused by the urethra that crosses the gland transversally along the z axis, following the conventions of Fig. 1.This creates a strong edge response from the center towards the anterior region of the prostate.Also, whereas edges are usually well defined around the mid-gland region, they are often very difficult to identify around base and apex, where the echoicity of the prostate becomes similar to surrounding tissues (Fig. 2bc).A less reported issue is the uneven reliability of prostate edges depending on their relative orientation to the beam.Edges orthogonal to the beam axis (y direction) tend to be less visible than those parallel to it (x direction).In a 2D transverse view, this translates by well identified edges in the anterior and posterior regions of the FOV, and by weak to invisible edges in the lateral regions.Another major difficulty of TRUS lies in the intersubject variability of the prostate characteristics.The shape and volume of the gland can vary substantially from one patient to another.Fig. 3a shows manual delineations of 2D midgland contours in our image data set (36 patients).In these patients, the prostate volume varied from 17 cm 3 to 54 cm 3 (mean: 41±10 cm 3 ).More importantly, distributions of image intensity show important variations depending on acquisition parameters (gain, probe quality, probe frequency) and patientspecific properties (tissue echoicity, calcifications).Fig. 3b shows nonparametric estimations of intensity ditributions inside the gland, limited to a representative subset of 14 patients for clarity.Mean intensity, standard deviation and global shape of the distribution strongly depend on the patient considered.

C. Computer-aided segmentation
Various segmentation methods have been proposed to perform TRUS-based PVS, which are summarized in comprehensive reviews of the literature [5]- [7].The lack of publicly available dataset with reliable ground truth and standardized validation metrics, and the diversity of acquisition setups found in the clinic make comparisons of results issued from different groups difficult.Up to this date, no method can be considered as a gold standard, due mainly to the variety of the difficulties encountered in this modality.Because of the lack of consensus segmentation strategy, the majority of PVS methods make use of a combination of several paradigms such as shape or regularization priors, image statistics, classification, clustering, or edge detection, alongside with variable amount of user interaction [6], [7].The objective of these "hybrid" approaches is to address the many challenges met in TRUS in a complementary fashion.Because of the different kinds of information used, these methods are often formalized under modular deformable model (DM) frameworks like active contours [8]- [15] or level-sets [16]- [22].In DM approaches, an initial geometric model is subjected to regularity (internal) constraints and a combination of image-derived (external) terms in order to converge towards the boundary of the prostate volume.The challenges of such models then lie in appropriately selecting and weighting ad-hoc external and internal terms while reducing the amount of patient-dependent parameter tuning and user interaction.DMs also typically evolve in non-convex energy landscapes, which makes their accurate initialization critical.

D. Outline of Present Work
In this article, we propose a new hybrid edge and regionbased parametric deformable model, or active surface, for PVS.Our contribution is threefold.First, we develop a new prostate edge detector derived from the radial bas-relief (RBR) approach [23].We generalize this concept to define analogues of image derivatives by decomposing the RBR in orthogonal directions, that we embed in the structure tensor formalism [24] to yield improved prostate edge detection.
Second, we combine this edge term with a region-based force driven by the Bhattacharyya distance [25] in a narrow band at the vicinity of the active surface.The idea of using the Bhattacharyya gradient flow for prostate segmentation was first proposed in [18] under the level-set framework.However, it was applied to match template prostate intensity distributions, which cannot be applied in presence of highly heterogeneous datasets such as ours.We derive a simplified expression for the Bhattacharyya flow that depends only on the distance between normalized kernel densities of inner and outer prostate regions.We note that a narrow band Bhattacharyya hybrid level-set model was proposed very recently for hip joint segmentation in MR images [26].However, our approach differs in that we substantially modified the formulation of this region term in the case of discrete parametric surfaces.
Finally, we develop a new initialization technique for prostate deformable models.Taking advantage of the watertightness of the prostate and its relatively low shape complexity, we analyze its radial edge profiles in polar coordinates to define initial landmark points towards which an initial surface model is warped.By doing so, we enable fast and reproducible initialization, requiring only that the prostate is approximately located at the middle of the field of view (FOV) and that the base and apex slice are selected.
We validate our method on a set of 36 3D TRUS images using manual delineations from two expert radiation oncologists, using a wide variety of quantitative metrics.
This paper is organized as follows.In the next section, we introduce related image processing concepts.In section III, we describe our proposed active surface model.The validation setup used in our experiments is presented in section IV.In section V, we present the results obtained on our clinical dataset.We finally discuss the results and conclude this article in section VI.

II. RELATED WORKS
In this section, we briefly discuss the general characteristics of deformable models and introduce related TRUS-based methods, motivating our implementation choices.
Deformable models (DMs) such as active contours [27] have been extensively used for TRUS-based prostate segmentation [19], [21].Their ability to combine shape and image information like edge or image statistics under a common framework make them appealing in a context where difficulties of different nature must be addressed.DMs are dynamic geometric models superimposed onto the image domain and evolving under the simultaneous effects of regularizing and image-derived terms.Their evolution can be considered from a variational point of view, in which an energy functional is iteratively minimized, or equivalently from a force-balance point of view, in which the net force acting on the model reaches zero at equilibrium.The model forces are then designed to equilibrate at the boundary of the region(s) of interest.In this article, we consider the force-balance point of view, and focus on the case of 3D object segmentation in volumetric images.
A DM can be represented either parametrically by a set of connected vertices defined in the subpixel 3D domain [28], [29], or implicitly as the zero-level set of a 3D hyper-surface defined in the volumetric image domain [30].In general, two advantages can be attributed to level-set formulations: ease of implementation, as calculations of derivatives are performed in the image domain, and topological flexibility of the model, allowing the surface boundaries to split or merge.This last property is not relevant to 3D PVS as the topology of the prostate is known to be a 1−connected (watertight) object, and may constitute a drawback if spurious changes in topology occur owing to e.g.noise.For this reason, in this work, we use parametric models where the 1−connectedness of the segmentation is guaranteed, following several works [9], [15].
Regardless of its implementation, the DM evolution can be formulated by an Euler-Lagrange evolution partial differential equation of a model S subjected to forces [31]: where S 0 is the initial model, subscript • t denotes partial derivative with respect to Euler time t, F int i are the internal forces controlling the smoothness of the model, F ext j are the image-derived external forces, and α i and β j are weights that control the relative influence of these terms.The design of deformable models thus consists in selecting the appropriate combination and weighting of model forces alongside with providing accurate initialization.We recall in the following paragraphs some of the choices made by previous TRUS-based PVS approaches.

A. Internal forces
Internal forces can be based on various smoothness measures.For example one can define analogues of respectively tension and rigidity using the first and second spatial derivatives of a surface model with respect to arclength parameters [27].Mean curvature motion [32] is also a popular regularizer of level-set approaches, where a force proportional to the local mean curvature is applied in the direction normal to the model.Tension and mean curvature force are analogous [31] and often found in practical implementations of deformable models for TRUS-based PVS [11].However, models subjected to these forces shrink and collapse to spherical points with increasing Euler time in absence of image forces [33].Special care must therefore be taken in order to obtain the desired smoothing effect while maintaining attachment to data.
Higher order forces provide softer regularization effect, being less prone to provoking the collapse of the model while preventing the apparition of irrealistic sharp discontinuities such as corners [34].A few authors proposed to implement higher order internal terms in TRUS-based PVS based e.g. on differences in local curvature within the framework of discrete explicit dynamic models [8] [14], [15], [35].

B. Statistical region-based forces
Complementary to internal forces, external forces constitute the attachment to image data part of the model.They can be divided into two categories: region-based forces and edgebased forces.Region-based forces (RBF) like the Chan and Vese (CV) model [36] take advantage of variations in image statistics between the prostate and the background by maximizing a statistical distance between the two regions, e.g.differences in mean intensity values or between-region variance.RBF are typically robust to noise, as they do not rely on local edge information.However, due to the global nature of the analysis, they may fail in handling intensity heterogeneities.Localizing statistical analysis within a neighborhood of each model node [37] can partly alleviate this problem.RBF-based models were applied to prostate imaging using global [22] and localized versions [19], [21] of the CV approach.Alternatively, other statistical distances were also considered [18], [38].In particular, Xu et al. maximized the Bhattacharyya distance between the background and a template prostate intensity distribution using the level set framework [18], yielding higher discriminative power compared to CV models [25].The Bhattacharyya gradient flow is applicable to a large category of probability laws, which makes it relevant for TRUS-based PVS where distributions of intensity do not follow a known pattern between patients.Our proposed approach includes a region term based on the Bhattacharyya distance.However, due to the intersubject variability of the prostate characteristics, as illustrated in Fig. 3b, we do not follow the strategy proposed in [18] where a template distribution is defined from training data and describe in section III-C an alternative approach.

C. Local edge-based forces
RBF may fail due to their sensitivity to image heterogeneities.For this reason, they are generally used in combination with local edge-based forces (EBF) in coupled hybrid edge-region frameworks.The rationale is to take advantage of robustness to noise through RBF while increasing accuracy around edges with EBF.
Classical gradient-based detectors are often used for prostate edge detection [6], [10], [13], [15], [19].However, due to the low quality of the gradient signal in ultrasound images, other edge detectors based on higher level information are also considered, such as coefficients of the dyadic wavelet transform [12], [22], gradient vector flow [39], instantaneous coefficient of variation [40], interacting multiple model (IMM) [41], or contrast between line segments along the direction normal to the model [14], [21], [42].We are particularly interested in this work in the radial-bas relief (RBR) approach [23], upon which our model is based.RBR exploits the fact that the prostate is located at the center of the FOV and appears darker than the surrounding tissues.By superimposing to the original image its radially enlarged and negated version, a binary edge map can be obtained after thresholding, showing prostate edge response where classical gradient-based detectors would fail (Fig. 4).However, RBR maps cannot be directly embedded in deformable models as they generally contain spurious edges and provide only rough estimates of the boundary.In section III-B, we propose a new edge detector based on the RBR approach and applicable to deformable models.

D. Initialization
Another critical aspect of deformable models is initialization, as they are prone to converge towards local minima along the course of their evolution.TRUS-based PVS methods range from manual interaction independent of image data [43] to fully automatic methods [22].A common intermediate practice is found in the use of semi automatic initialization, where landmark points are selected on the prostate boundary usually around the mid-gland region and at base and apex [11], [15].In section III-E, we propose a new quasi-automatic initialization method based on our modified RBR edge detector, where only base and apex slices need to be identified.

E. Alternative PVS strategies
Other strategies for PVS were proposed based on multiple 2D segmentations propagated from one slice to the other across the volume [10], [19], [21], [22], [44], [45].Wang et al. propagated an initial segmentation result obtained at mid-gland, where the prostate is easily detected, towards basal and apical slices [10].Due to poor performance in peripheral slices, the authors also proposed to reslice the prostate volume in a rotational way to take advantage of the relative symmetry of the gland, lowering variations between subsequent slices.This idea was followed and adapted by other authors [19], [44], [45].Reslicing approaches are however prone to accumulation of error along the propagation.Also, rotational reslicing is applied to isotropic 3D volumes with equal pixel size.In many clinical setups including ours, volumes are anisotropic as they are reconstructed from multiple subsequent 2D images acquired with a mechanical stepper [1], which lead to much larger spacing along the stepper axis than across the 2D planes, thereby reducing the practicality of rotational reslicing.

A. General formulation
We propose a new hybrid deformable model driven by three competing forces : rigidity internal forces, region-based forces based on the Bhattacharyya distance, and edge-based forces based on a modified radial bas-relief edge detector.
Let S(m, n, t) be an active surface represented as a mapping of a bivariate parameter (m, n) on a grid superimposed onto the image domain, where t is the evolution Euler time.The general formulation of our model is: where •, • is the dot product, N is the normal direction to the contour, F rigid is the internal rigidity force, F region is the regionbased force, F edge is the edge-based force field, and α i are the corresponding weighting parameters of these forces.The three force terms are described in more detail in the remainder of this section.

B. Edge Detection with Bas-Relief Structure Tensor
We introduce a new prostate edge detector for TRUS inspired from the radial bas-relief (RBR) approach.We first recall the different steps of the original RBR method in 2D images.
Let I 2D (x, y) be a 2D TRUS image with intensity values scaled between 0 and 1 and showing the prostate at the center of the FOV of coordinates O 2 = (0, 0).A scalar (grayscale) radial bas-relief image R 2D is obtained by summing the original image with its negated, radially enlarged version: where r is a scale factor, usually set to around 110% [23].
To obtain a binary edge map, the radial bas-relief image is processed in several steps summarized in the left part of Fig. 5, which can be formulated by: First, pixels with R 2D (x, y) > 1 are set to 0 through operator T 1 .The resulting image is then binarized using a 2 class Otsu threshold operator O, the output of which is in turn smoothed by successive applications of dilations and erosions using a morphological operator M. Finally, a connected components operator C n suppresses image islands with a number of voxels below a specified value n, producing binary RBR edge maps as illustrated in Fig. 4. In this work, rather than seeking to obtain a binary edge map, we derive a scalar edge detector for 3D TRUS images from the bas-relief principles.The different steps of the procedure are described in the right part of Fig 5.
Let I(x, y, z) be a 3D TRUS image.We first extract orthogonal bas-relief components by enlarging the image along the x, y and z directions: Bas-relief maps R x , R y and R z show maximum intensity near the prostate boundary along each axis.Figs.6(a-d In order to derive a scalar edge detector for TRUS, we define a structure tensor built upon the orthogonal bas-relief components.The structure tensor, or second moment matrix, is a well studied object in image processing that characterizes local image orientations [46], [47].It is expressed as: where * denotes convolution, • T is the transposition operator and G ρ is a Gaussian kernel of scale ρ.Gaussian regularization yields more homogeneous orientation estimation across space and helps partially filling gaps of scale ρ in regions where gradient information is missing. Considering that bas-relief maps encode amplitudes of image variations in the three directions of space, we propose an analogous tensor expression based on bas-relief information, the Bas-Relief Structure Tensor (BReST): The largest (or principal) eigenvalue λ 1 of B ρ provides a measure of the amplitude of maximum variations in the image [48], that we use as our edge detector.Fig. 6e shows example edges obtained from the BReST.The edge map can be refined through additional processing.For example, boundary gaps can be further filled using a 2D stick filter [49] applied on a slice by slice basis (Fig. 6f).
For this edge information to be used by a deformable model, we derive a 3D edge-based force field from the principal eigenvalue λ 1 .To this end, we use the Vector Field Convolution (VFC) method [50].VFC is a computationally efficient generalization of the popular Gradient Vector Flow [31], which produces smooth vector fields oriented towards image edges.A VFC field is obtained by convolving a scalar edge map with a Vector Field Kernel (VFK), a vector kernel of size s whose vectors point towards its center with decreasing magnitude d −γ , where d is the distance to center and γ is the attenuation parameter.
The proposed external force field is expressed as: where K is the VFK.Normalization is performed so that force vectors have unit Euclidean norm, the strength of the force applied to the model being modulated by parameter α 3 of eq. (2).

C. Region-based force using the Bhattacharyya distance
We now consider the statistical region-based force acting on the deformable model in combination with the edge-based force.
The oriented active surface S divides the image I(x) = I(x, y, z) into two mutually exclusive sub-regions Ω in and Ω out .The sub-region Ω in corresponds to voxels inside the surface that approximate the prostate volume, whereas Ω out approximates the background.Following an idea originally proposed in the level-set framework [25], our statistical region force is based on the maximization of the Bhattacharyya distance between distributions of intensity in these two image regions.
The Bhattacharyya distance is computed between two probability density functions (PDF) inside and outside of the surface model p in (k|S) and p out (k|S), k ∈ R N , where R N is a N -dimensional feature space embedding N image features such as intensity, gradient or texture information.It is defined as − log B , where: is the Bhattacharyya coefficient.In this work, we take N = 1 and consider image intensity as the only feature.A kernel density estimation of the PDF of I can be written as: where K σ is a Gaussian kernel of bandwidth σ.A Bhattacharyya region force is derived by taking the first variation of B with respect to S [25], [31]: where A in and A out denotes the volumes (number of voxels) of Ω in and Ω out respectively and is the weighted difference between the square roots of the likelihood ratios p in /p out and p out /p in .
The Bhattacharyya gradient flow of eq. ( 12) was originally proposed in the level-set framework [25], in which inner and outer regions are defined straightforwardly, using the signed distance to the boundary of the deformable model as the levelset function.Our surface model is implemented as a parametric surface, constituted of a set of vertices V i (x, y, z) defined in the subpixel domain.The equivalent implementation of eq. ( 12) in an explicit framework would require computationnally expensive surface rasterization such as found in [29] in order to determine inner and outer image voxels.In this work, we consider an alternative strategy in which we simplify this statistical region term using several approximations.
We first assume that the initialization of the DM is sufficiently accurate so that we can restrict statistical analysis in regions Ω in and Ω out to two approximately parallel narrow bands Ω b in and Ω b out of thickness b mm around S [51].This has the effect of eliminating the influence of prostate and background heterogeneities beyond the narrow band.
To further reduce the sensitivity of the approach to heterogeneities, we increase the local intensity contrast by performing local image normalization.The image I is normalized by applying the following linear transformation: where local means µ I and standard deviations σ I are estimated using an averaging window of size L x × L y × L z mm 3 .We also suppose that the inner and outer narrow bands are of approximately equal volume, i.e.A in ≈ A out = A.This simplification makes the region force independent on the differences between A in and A out and eq.( 12) boils down to: To compute this region force in a parametric framework, we approximate the PDFs in the narrow band p in (k|S) and p out (k|S) by sampling intensity values at regular intervals along the normal to the surface N i at each vertex V i using bicubic interpolation.We then normalize F region between −1 and 1 to obtain the final region force applied to the model, weighted by parameter α 2 of eq. ( 2).Fig. 7 shows a closeup representation of the narrow band near the boundary of the deformable model, where sub-regions Ω in and Ω out are reduced to banded regions Ω b in and Ω b out .The white dots along the inner and outer normal at each vertex represent the different sampling points at which intensity is interpolated to sample p in (k|S) and p out (k|S).Fig. 8 shows the effect of local normalization on the sampled density estimations p in (k|S) and p out (k|S) in the 14 patients shown in Fig. 3.As previously, in this example the surface model S considered is the manually delineated ground truth surface.The local normalization reduces the variability of the distributions among subjects while preserving shape differences between inner and outer regions.

D. Internal force
We consider in this work high order rigidity forces acting on the deformable surface model to increase resistance to twist while reducing its tendency to collapse to a point.The internal force F rigid is expressed as: where α 1 is the weighting parameter of eq. ( 2) and ∆S is the surface Laplacian: The surface Laplacian is approximated at each vertex V i of the triangulated surface using the umbrella operator L(V i ) [52]: where n is the valence (number of neighbors) of the neighborhood N i of V i .We implement the model evolution using a classical semiimplicit finite different scheme, where internal force motion is computed implicitly and external force motion explicitly [27].

E. Model initialization based on polar edge analysis
The initialization of deformable models is a critical part of their design, as they are prone to fall into local minima along the course of their evolution.In this section, we propose a new quasi automatic initialization technique for 3D TRUS images based on the modified bas-relief edge detector described in section III-B.We analyze radial edge profiles in polar coordinates to identify 2D landmark points in the (x, y) midgland plane towards which an initial prostate model is warped.
The proposed approach requires only that the base and apex slices z base and z apex are determined by the user and that the prostate center is approximately located at the center of the FOV, of coordinates (x 0 , y 0 ) in the transaxial slices.
First, using the edge detection procedure defined in section III-B, we obtain 3D edge maps corresponding to the principal eigenvalue λ 1 (x, y, z) of the BReST B ρ .
Given knowledge of the base and apex slices z base and z apex , we identify the mid-gland slice as z 0 = (z base + z apex )/2.As the prostate shape varies slowly along the transverse axis near the mid-gland, where edges are usually well-defined, transverse slices in the vicinity of z 0 also contain relevant edge signal.For this reason, we strengthen the mid-gland 2D edge map λ 1 (x, y, z 0 ) by integrating λ 1 along the z axis in the neighborhood of the central slice: where the value δ z was empirically chosen to two consecutive slices of 1 mm in our experiments.
We then look for maximum edge responses of λmid (x, y) along radial line profiles using equal angle sampling around the center of the FOV (x 0 , y 0 ) to determine locations of candidate landmark points (Fig. 9a).To this end, we transform λmid into polar coordinates using bicubic interpolation.The angular step was set in our experiments to α = 1 • , yielding 360 candidate landmark points.Fig. 9b shows an example mid-gland edge map λmid where landmark candidates are superimposed.Maximum radial edge response is obtained mostly around prostate edges.As shown in this example, some outlier points can appear for angles where the edge signal is too weak.In order to remove these outliers, we apply a point cloud denoising algorithm based on neighborhood processing using the mean µ p and standard deviation σ p of the K nearest neighbour distances of each candidate point [53].An edge point p is considered an outlier if the average distance to its K nearest neighbour is outside µ p ± mσ p , where K = 30 and m = 0.25 are two parameters determined empirically.Fig. 9c shows the effect of outlier rejection, where remaining points are located on the prostate contour.The points are then connected to form a 1D curve C(x, y).To obtain a smooth estimate of the midgland contour without shrinking the curve, we resample C to provide even spacing between its nodes and regularize it using the one-dimensional analogue of the internal force described in section III-D for two-dimensional surfaces [27].Fig. 9d shows the final curve obtained.To enforce the robustness of the proposed approach to outliers, the mid-gland landmark candidates are reduced to 24 points by taking the average coordinates of the points belonging to C in sub-quadrants of 15 • along the contour (Fig. 9d).Two additional landmark points are finally selected at base (x 0 , y 0 , z base ) and apex (x 0 , y 0 , z apex ).
We obtain an initial surface model S 0 defined in the volumetric image domain by warping an oblate ellipsoid towards the landmarks using thin plate spline warping with point-to-point correspondence [54].We superimpose to the volumetric image a triangulated surface model E 0 of an oblate ellipsoid centered at (x 0 , y 0 , z 0 ) whose semi-principal axes (a x , a y , a z ) are defined so that E 0 is enclosed within the 26 image landmark points.In our experiments, we set (a x , a y , a z ) to 2/3 of the maximum image landmark coordinate in the corresponding dimension.Source vertices are identified on the ellipsoid model that correspond to the 26 destination landmarks.Fig. 10 shows tridimensional views of the initial ellipsoid and the result of thin plate spline deformation for the image shown in Fig. 9.The manually delineated surface is also shown for comparison.

A. Image acquisition
Transverse TRUS images were collected at the Brest University hospital, Brest, France during the examination of 36 patients undergoing preimplant brachytherapy planning (2231 transaxial 2D images in total).For each patient, series of transverse 2D TRUS images were acquired at 1-mm intervals using a BK Medical Flex Focus 500 ultrasound station using a Biplane Endocavity Z848 TRUS probe.The pixel size of each transverse slice was 0.152 mm×0.156mm.The TRUS probe was placed on a stepper, allowing to move along the transverse axis to acquire parallel transverse slices that covered the entire prostate volume.
Manual delineations of the prostate boundary were performed in every transverse slice by two radiation oncologists from the Brest University Hospital with full expertise on the brachytherapy protocol, hereafter referred to as as E1 and E2.These reference contours were used to reconstruct 3D prostate volumes that we considered as surrogates of ground truth to assess the accuracy of the proposed approach.

B. Pre-processing
Due to the protocol followed at the Brest University hospital, a point grid was overlaid to B-mode acquisitions in all the transverse images of our data set, visible for example in Fig. 2c.To suppress its influence, we pre-filtered all the image slices with a 2D median filter of dimension 5 × 5 pixels.

C. Quantitative metrics
We performed quantitative reference-based comparative studies to evaluate the segmentation accuracy of the proposed hybrid deformable model.Given the lack of consensus validation metrics [6], we considered various volume-based and distance-based metrics to ease comparison with other works.For volume-based metrics, we considered the Dice Similarity Coefficient (DSC) [6], Percent Volume Error (PVE) [41], Percent Volume Difference (PVD) [41] and Sensitivity (SE) [6].For distance-based metrics, we considered the Hausdorff distance (HD) [6] and the mean absolute distance (MAD) [6], expressed in milimeters.Descriptions of each metric can be found in the corresponding referenced article.
We divided the gland into three sectors: base, mid-gland and apex to provide more detailed analysis of the performance of the method in the different regions of the prostate [41], [45].We used the following proportions along the transverse axis: 30% for base, 40% for mid-gland, and 30% for apex.

D. Parameters
As for any deformable model, the proposed hybrid active surface is associated with parameters that are summarized and briefly described in Table I.In our experiments, these parameters were set equally for all images and determined in order to maximize the DSC score on average among the 36 patients of the data set through empirical parameter tuning, where Expert E1 was used as a reference.To provide better insight on how these parameters may be adapted to specific cases, we also give indicative ranges of values for some influential parameters.The ranges shown are minimum and maximum values which led to a maximum DSC score, i.e. superior to the average score obtained with the values set for our experiments, in at least one image of the data set.

V. RESULTS
We implemented the method using MATLAB.The average processing time, including initialization and excluding base and apex slice selection, was 40 ± 4 seconds on average using a Intel Core i7 CPU.This time is compatible with the clinical requirements of the LDR-B procedure.
Fig. 13 shows segmentations results obtained in two different images of the data set.The manual reference volume obtained through successive delineations of the prostate in every transverse slice is superimposed for visual comparison.
Table II summarizes the quantitative results obtained for the entire data set using the proposed active surface with respect to the two expert delineations E1 and E2.Results showing the agreement between the two experts are also shown, referred to as "E2 vs E2".Expert E1 is arbitrarily used as the reference for non-symmetric metrics (PVD and SE).
Segmentation results are presented in more detail in Fig. 12 for each patient and each metric, where Expert E1 is taken as  the reference.To compare the relative improvement achieved by the active surface with respect to the proposed initialization stage, average segmentation results obtained with the initial model before its evolution under eq.( 2) are shown in Table III, where E1 is used as a reference.After surface evolution, an average global DSC with respect to expert delineation E1 (resp.E2) of 0.92 ± 0.02 (resp.0.91 ± 0.02) was obtained across the entire prostate volume.In the most central slices representing 40% of the gland, the approach led to a DSC of 0.96 ± 0.01 (resp.0.96 ± 0.02) on average.Scores were lower at the gland extremities, with comparable average DSC of 0.87±0.04around the base (resp.0.84 ± 0.07) and of 0.86 ± 0.07 (resp.0.85 ± 0.08) around apex.By comparison, the initial model achieved a global DSC accuracy of 0.87 ± 0.03 (base: 0.79 ± 0.06, mid-gland: 0.93 ± 0.03, apex: 0.81 ± 0.06) with respect to expert E1.
The general ranking between base, mid-gland and apex sectors is confirmed by the other metrics.The two contourbased metrics (Hausdorff distance and MAD) were consistent with volumetric metrics, with higher MAD and Hausdorff in the peripheral slices.PVD scores assess the relative volumetric deviation of the segmentation to the reference.No significant deviation of PVD from 0% was obtained at mid-gland.However, a general overestimation of the basal region led to negative PVD values.Sensitivity was high around mid-gland and base, whereas it was lower around apex, further suggesting a systematic underestimation of this sector.On the other hand, the method tended to underestimate the apical sector of the prostate.This observation is further supported by individual results of Fig. 12e, where it can be better appreciated.The underestimation of the apex is particularly visible in the example results shown in Fig. 13, where reference E1 was not fully enclosed by the active surface around the top of the image.
The comparison between the two expert delineations E1 and E2 followed a similar pattern.For instance, higher overlap was found in the central slices (DSC of 0.96 ± 0.02) than in the peripheral slices (DSC of 0.87 ± 0.04 for base and of 0.89±0.05for apex).Sensitivity of E2 with respect to E1 was especially low around the base, suggesting a general tendency of expert E1 to label more voxels as prostate in these regions.Interestingly, quantitative scores obtained by comparing E1 to E2 were not significantly better than those reported with the proposed approach.
The only step of the workflow that is not automated is the selection of the basal and apical slices, on which depends the construction of the initial model.We therefore studied the robustness of initialization to potential errors made in selecting these two parameters.To this end, we perturbed the choice of base and apex slice and studied the corresponding segmentation accuracy achieved by the initial surface.Fig. 11 shows the average DSC score of the initial surface as a function of the distance ∆ Base , ∆ Apex , in number of slices, between the actual base and apex obtained from ground truth and the slices selected during the initialization stage.For clarity, we restrict the analysis to outward expansion from the center, i.e: more background slices are selected with increasing values of ∆ Base or ∆ Apex .As expected, segmentation accuracy

VI. DISCUSSION
Comparing the above results with other PVS approaches is a difficult task, as there is still today a clear lack of public labeled data set allowing to perform fair comparative evaluations between available methods.We followed the strategy adopted by the majority of recent works in the field by reporting segmentation performances based on a wide variety of quantitative metrics.In the current state of PVS validation, one must be careful when concluding on the absolute superiority of one approach based on these values, as there exists a large discrepancy in validation setups between centers.
Nevertheless, considering values reported in other recent works [6], [7], [15], [45] we confidently situate the proposed approach among the state of the art in prostate TRUS segmentation.In particular, the quantitative scores obtained in the central sector of the prostate, constituted of 16 ± 3 slices on average, are comparable to results reported in single 2D midgland slices in a recent work on 2D prostate segmentation [20].High segmentation performances in this sector of the prostate can be stressed by considering only the mid-gland slice, where the method obtained area DSC scores of 0.98±0.01 on average (not shown in the results table).
Reference delineation was performed manually on a sliceby-slice basis, which is the most common clinical practice.It is therefore a very tedious task which cannot be achieved in reasonable time without bringing uncertainty on the resulting surrogate of ground truth.This is particularly visible in Fig. 13, where a spurious "staircase" effect can be observed in the reference volume due to rough shape variations from one slice to the other.On the other hand, our model being a parametric active surface, it inherently produces closed and smooth prostate contours, which approximate better the natural shape of the gland.
Due to the invisibility of prostate tissue in ultrasound images around base and apex, lower performances of our approach in these regions relative to the central slices were expected.The two expert delineations E1 and E2 also showed poor agreement  in these sectors.Importantly, the global agreement between E1 and E2 was comparable to results obtained using the proposed quasi-automatic approach.
The general overestimation of the basal region led to higher sensitivity and lower DSC.This is partly explained by the fact that, depending on user expertise and patientdependent requirements, the manually segmented basal region can include more or less of the adjacent seminal vesicles.A common clinical practice is to retain only the part closest to the prostate, to an extent that is dependent on user experience.In our experiments, expert E1 included more of the seminal vesicles than expert E2, which translated in higher overlap with segmentation results.As the seminal vesicles are indistinguishable from the prostate tissue before they separate from the gland further along the z axis, it is unlikely that any current automatic method could accurately discriminate them from the prostate.Nevertheless, to increase performance in peripheral regions, the proposed model can be made more accurate by e.g.setting more landmarks or attractive anchor points around the base and apex [19].However, we believe that such additional complexity would have shadowed the main contributions of this article.
We proposed to combine a new prostate edge detector based on the Bas-Relief Structure Tensor (BReST) with a Bhattacharyya gradient flow using parametric active surfaces.The external forces obtained can be applied to any type of deformable model to facilitate its evolution towards the boundaries of the prostate [31].We also developed a robust initialization strategy based on the BReST edge detector, which only requires the selection of the basal and apical slice.The proposed initialization stage allows for high segmentation accuracy in the central sector of the prostate (table III).This was expected as our initialization method exploits the detection of prostate edges at mid-gland.Naturally, initialization accuracy is lower in peripheral slices where no image data is used.The hybrid model then evolves from this initial stage to improve segmentation accuracy in all three sectors Fig. 15.Global DSC score achieved per patient in extreme cases where only the edge force 1  1 , 0 1 or the region force 0 1 , 1 1 is used to control model evolution.Expert E1 was used as a reference.
of the gland.The proposed method achieved state of the art segmentation results despite the variablity of our dataset in terms of prostate shape, edge strength or global statistics.
The major burden of deformable models is parameter tuning.Due the highly heterogeneous nature of the data set at our disposal, where images were acquired with variable probe parameters, it was difficult to obtain standard values leading to nearly optimal segmentation in all subjects.Ad-hoc parameter values were given in table I which led to the average results presented in this paper.The values reported are likely to allow for similarly accurate segmentation in other images.However, parameters may naturally be tuned to each subject for increased performance.The parameters most influential on model evolution are α 2 and α 3 , which control the relative strength of the edge-based force (EBF) and of the region-based force (RBF).In order to obtain the best weighting scenario between these forces, we studied their respective influence on segmentation accuracy.To this end, we kept the total external force weight α † = α 2 + α 3 constant and studied the average DSC results obtained by varying the ratios α 2 /α † and α 3 /α † between 0 and 1. Fig. 14 shows the corresponding results by steps of 1/6.Poor performance was observed when the model relied only on the RBF (Fig. 14, leftmost point).This is mainly caused by the short attractive range of the RBF due to the limited spatial extent of the narrow band.When the statistical distance between the inner and outer regions of the narrow band is due to within region heterogeneities, the model may evolve in arbitrary directions, which can lead to erroneous segmentation.The BReST-based EBF, of greater capture range, allowed for accurate segmentation results without the help of the region term (Fig. 14, rightmost point).The best weighting was nevertheless achieved when both forces were considered (α 2 /α † = 1/3 and α 3 /α † = 2/3), with notable improvement brought by the RBF around the basal and apical regions, wherein edge information is of lower quality.To further support the use of a hybrid model, DSC results corresponding to the two extreme scenarios of full RBF without EBF and full EBF without RBF are shown for every patient in Fig. 15.Cases for which the surface collapsed in the full RBF scenario (DSC scores of 0.6 or less) penalized the average DSC scores shown in Fig. 14, shadowing the fact that it achieved comparable or better performance than the full EBF configuration in approximately half of the subjects.Both these terms are therefore in general useful and complementary, despite cases where, for example, spurious or missing edge information would lead to erroneous edge-driven forces, or where statistically indistinguishable background and prostate would cause their merging by the Bhattacharyya gradient force.
The new BReST-based edge detector enables more robust prostate edge detection compared to classical gradient-based approaches.Its only requirement is that the prostate must be located approximately at the center of the FOV, which is always the case in clinical practice.Additional refinement for the calculation of the BReST can be considered.For example, the morphological processing step used in the original radial bas-relief approach can help provide cleaner edge maps.Also, the edge influence on model evolution can be weighted in inverse proportion with the distance to mid-gland, as edges are usually better defined at mid-gland than in peripheral slices.Preliminary explorations indeed showed some improvement by considering such refinements.The initialization step can also be improved by using more elaborated initial surfaces than the proposed ellipsoid.For example, tapered ellipsoids [41] or template models built from average co-registered segmentations can be considered.We assumed that the midgland image is located exactly at equal distance between base and apex.This may actually not be true, and the best contrasted slice was often found closer to base.However, we consider that the results shown with the proposed approximations are sufficient to highlight the advantages of the proposed hybrid model without additional and more complex refinements.More generally, there is always a compromise to be found between the amount of complexity specific to the problem at stake and the generality of the methods developed.

VII. CONCLUSION
We have proposed a parametric deformable model for accurate prostate volume segmentation in TRUS.The proposed approach is based on a combination of new edge and region image forces tailored for prostate imaging.To reduce the variability of results induced by initialization, we have also proposed a quasi automatic initialization scheme based on the analysis of prostate edges at mid-gland.Results using fixed parameters on a large dataset of 36 patients demonstrated the relevance of the proposed approach for achieving state of the art segmentation performance with very limited user interaction.We hope in future works to automatize parameter tuning based on image data to optimize performance on a case-by-case basis and avoid cumbersome and time-consuming exploration of the parameter space.

Fig. 1 .
Fig. 1.Schematic representation of a typical brachytherapy procedure.Radioactive implants (seeds) are inserted under TRUS guidance.Multiple parallel transverse slices in the xy plane are acquired with a motorized stepper to cover the entire prostate volume.Modified image by Cancer Research UK, distributed under a CC BY-SA 4.0 license.

Fig. 2 .
Fig. 2. Limitations of TRUS prostate volume segmentation.(a) General overview.(b) Example basal slice.Manual delineation by a radiation oncologist is superimposed as a red contour, showing the common practice of including a part of the seminal vesicles (light gray region near the bottom).(c) Example apical slice, showing poor to invisible prostate contours.

Fig. 4 .
Fig. 4. Radial bas relief edge maps obtained by superimposing to the TRUS image its enlarged, negated version.Top: Original TRUS image.Middle: Edge detection with gradient magnitude.Bottom: RBR edge map.

Fig. 5 .
Fig. 5. Left: Workflow of the original RBR approach.Right: Workflow of our proposed method for scalar edge detection in TRUS images.
Fig. 6. (a-d) Comparison of edge detection using gradient amplitude and bas-relief.(e) Principal eigenvalue of the Bas-Relief Structure Tensor (f) After stick filtering (stick length: 101, stick width: 1).Gray levels were scaled for better clarity.

Fig. 7 .Fig. 8 .
Fig. 7. Estimation of inner and outer probability density estimations of image intensity in a narrow band around the surface model (solid white).Image values are sampled at regular intervals (white dots) along the normal to each vertex (yellow nodes) to obtain the sampled intensity distribution.
Fig. 10.3D illustration of the warping of an ellipsoid surface model towards image landmarks points.The initial oblate ellipsoid is shown in wireframe, the warped surface in yellow and the manually delineated ground truth surface in red.

Fig. 11 .
Fig.11.Sensitivity of initial segmentation results to perturbation of base and apex slice selection.Results shown are expressed in terms of mean DSC as a function of the distance between manually selected slices and the actual base and apex, as determined by expert E1 .A surface is interpolated between the points to ease visualization.

Fig. 12 .
Fig. 12. Quantitative segmentation results obtained on the 36 patients data set, where expert E1 was used as a reference.Results are divided into base, mid-gland and apex.

Fig. 13 .
Fig. 13.Example segmentation results.Active surface is shown in yellow, surrogate of ground truth based on slice-by-slice delineation (by expert E1) is shown in red.

TABLE I LIST
OF PARAMETERS ASSOCIATED WITH THE METHOD.

TABLE II SEGMENTATION
ACCURACY FOR 36 PATIENT 3D TRUS IMAGES (MEAN±STANDARD DEVIATION).