Autonomous Prostate Segmentation in 2D B-Mode Ultrasound Images

Abstract

Prostate brachytherapy is a treatment for prostate cancer; during the planning of the procedure, ultrasound images of the prostate are taken. The prostate must be segmented out in each of the ultrasound images, and to assist with the procedure, an autonomous prostate segmentation algorithm is proposed. The prostate contouring system presented here is based on a novel superpixel algorithm, whereby pixels in the ultrasound image are grouped into superpixel regions that are optimized based on statistical similarity measures, so that the various structures within the ultrasound image can be differentiated. An active shape prostate contour model is developed and then used to delineate the prostate within the image based on the superpixel regions. Before segmentation, this contour model was fit to a series of point-based clinician-segmented prostate contours exported from conventional prostate brachytherapy planning software to develop a statistical model of the shape of the prostate. The algorithm was evaluated on nine sets of in vivo prostate ultrasound images and compared with manually segmented contours from a clinician, where the algorithm had an average volume difference of 4.49 mL or 10.89%.

1. Introduction

Prostate brachytherapy is a procedure to treat prostate cancer where long flexible needles, containing radioactive seeds, are inserted into the prostate. This procedure is an attractive treatment option used in early-stage locally-confined prostate cancers [1,2]. In current practice, prostate brachytherapy requires a planning stage where the entire volume of the prostate is imaged and the 3D contour of the prostate is created. The 3D prostate contour is used for dosimetric calculation and, traditionally, is manually segmented by a clinician in software. While either MR or CT imaging modalities can be used for the dosimetric planning, the typical imaging modality used is ultrasound, both for cost-effectiveness and to match the imaging modality used by the clinician during the insertion procedure. This manual segmentation of the prostate from axial image slices is time-consuming, and thus algorithmic prostate segmentation poses an attractive option to reduce the amount of time that a clinician is required during prostate brachytherapy preplanning. In addition to therapeutic intervention planning, prostate gland contours can also be used for medical diagnostics, through measuring parameters such as gland asymmetry [3] or gland image appearance [4].

To assist a clinician when segmenting the prostate, we propose a fully-autonomous segmentation algorithm, shown in Figure 1, that consists of three main parts: a statistics-based active shape model, superpixel image clustering, and graph-cut-based segmentation. First, a 3D prostate contour model is defined that will be used to model, and later segment, the shape of the prostate in TRUS images. Utilizing a dataset of previously gathered manually segmented prostate contours, a Gaussian-distribution-based statistical model of prostate shapes will be created. This Gaussian distribution model thus forms the active shape prostate model and is used to estimate the range of values (and corresponding probability) for each parameter in the shape model. With the active shape model found, the next step in the prostate segmentation algorithm involves the image processing of a full set of axial images for a patient (capturing the 3D shape of the prostate). A novel superpixel algorithm is used on each TRUS image in the set to cluster contiguous pixels together, based on statistical similarity, into large regions of superpixels. By utilizing the statistical dissimilarity between neighbouring superpixel regions after processing, the edges of the prostate in each axial image can be found.

The final 3D prostate contour segmentation is carried out through a graph-cut based method, to find the most likely prostate contour based on the active shape model and informed by edge features detected by the superpixel algorithm. This algorithm is thus a hybrid method which incorporates information from an active shape model into a graph-cut based segmentation system. The work of [5] proposed a similar style of hybrid system, where an active contour model is progressively moved and deformed through graph-cut-based optimization. This hybrid method was shown to be robust to image noise and large edge discontinuities. This robustness to edge discontinuities is of particular interest in medical imaging, where [6] applied a similar hybrid active contour graph-cut model, and used a series of graph-cut-based expansion and contraction steps to move the contour and segment out a kidney in CT image slices. In a similar vein to our proposed method, ref. [7] proposed augmenting an active shape–active appearance model (incorporating both contour and image statistics for the organ of interest) with a graph-cut-based optimization to delineate 3D contours of livers, kidneys, and spleens in CT image volumes. This hybrid segmentation methodology has the significant advantage that it can match the “contouring style” of a particular clinician and is significantly robust to the poor contrast and prostate edge definition seen in TRUS prostate images.

Contributions and Outline

The proposed prostate segmentation algorithm within this paper specifically provides the three following contributions:

• It incorporates a novel superpixel-based pixel clustering algorithm that includes a secondary graph-based hierarchy to allow for graph-cut image segmentation at the superpixel level.
• A probabilistic prostate model is introduced, and through the use of this probabilistic modeling, the algorithm is highly flexible and can produce contours which statistically resemble the input manual contours from a particular clinician or group of clinicians.
• The prostate contouring algorithm is fully autonomous, and does not require any user input for the segmentation process, nor any manual post-algorithm contour correction.

This paper is laid out as follows. A brief overview of the literature is provided in Section 2. The prostate contour model and statistical prostate shape model are presented in Section 3. The proposed superpixel algorithm for the ultrasound image processing is covered in Section 4 and the graph cutting segmentation is discussed in Section 5. The experimental setup and evaluation of the proposed algorithm as applied to in vivo patient images are given in Section 6. The conclusions and results are summarized in Section 7.

2. Background

There are several artifacts and limitations inherent in ultrasound images that must be compensated for in computer-aided tissue segmentation [8]. In the specific case of prostate segmentation, the well-known shadowing and edge enhancement artifacts, which obscure the edge of the prostate in US images, have motivated the use of active shape, active contour, texture-based, and other methods for segmentation. Regardless of segmentation modality, algorithmic prostate segmentation can be grouped into two general methods, semi-autonomous and fully-autonomous systems. Semi-autonomous systems require some user input for each patient in order to perform prostate segmentation. Fully-autonomous systems, in contrast, require no such user input. The field of prostate segmentation in multiple imaging modalities, including ultrasound, CT, and MRI, has been well explored in the literature, with [9,10] providing comprehensive surveys. A more recent survey focusing on the automated diagnosis of prostate cancers based on algorithmic prostate contouring in ultrasound images is given in [11].

In active contour model segmentation, a simple contour model is deformed through some forcing function such that the contour will come to rest on the edge of the prostate. Active contour segmentation is typically semi-autonomous, requiring user input for the initial contour placement, and can be particularly sensitive to areas of low contrast or shadowing/enhancement in US images. One of the first active contour segmentation systems was presented by [12], wherein an initial user-defined 2D contour model (later extended to 3D by [13,14]) is fit to the edge of the prostate, with the user editing the segmented contour in areas affected by US artifacts. To compensate for incomplete edge information algorithmically, ref. [15] used an active contour model which incorporates a Kalman filter to probabilistically detect the edge of the prostate around user input seed points. To avoid the instability caused by the forcing functions in standard active contour model methods, ref. [16] enforced constraints on the contour utilizing a 1D wavelet transform, and [17] used nonlinear warping of the axial image to an ellipse-based contour model. The work of [18] represented the contour of the prostate with 2D superellipses. This concept was further developed in [19], where information from a full series of these 2D axial superellipses was refined into a single 3D ellipsoid contour.

To enhance the robustness of contour-based segmentation algorithms in the presence of shadowing/enhancement and low organ-to-background contrast, a statistics-based approach known as active shape modelling has been employed. In active shape segmentation, a large number of manually segmented contours (or feature points) are input into a system to derive the probability that the contour will take a particular shape (or has particular feature point locations) [20]. By incorporating both contour and image information from manually segmented contours in a hybrid method, refs. [21,22] used non-rigid registration to warp an input US image to all of the manually segmented US images in a database. A series of warped manually-segmented images which strongly matched the input image is then used to calculate the contour statistics, with the mean of this statistic being used as the contour for segmentation. The work of [23] presented a 2D prostate contouring algorithm using weak membrane that is deformed to segment the prostate based on edge information, enhanced by an edge-preserving filter, in TRUS images. A spline-based discrete dynamic contour model was presented in [24] and was used by [25] to segment the prostate after filtering with a dyadic wavelet transform.

Using the statistical inference from a series of prostate contours, ref. [26] created a point distribution model and [27] used a spherical harmonics model, which were fit to a series of user input points and then updated based on the likelihood that nearby prostate edges in the US image correspond to the statistical model. A Gabor filter method was proposed in [28] to create a smoothed and rotation-invariant shape statistical model of the prostate. Machine learning methods have also been used to create discrete statistical models of the prostate, which can then be actively deformed based on edges as found through multiresolution image-pyramid analysis [29]. In a similar manner, an adaptive learning method was used to refine patient-specific shape statistics in an online manner during segmentation using deep learning methods on a large TRUS image dataset [30].

3. Active Shape Prostate Contour Model

For planning in a clinical prostate brachytherapy procedure, a series of ultrasound images are taken as a transrectal ultrasound (TRUS) probe is translated forward along the z-axis (towards the head of the patient); see Figure 2. The prostate images are taken at regular distances apart to obtain a 3D volume and shape of the prostate. The goal of this work is to automate the 3D prostate contouring process through segmentation of the prostate in all of the TRUS image slices, ${}^P\mathbf{I}$, for a particular patient. Each TRUS image, denoted by ${}^P{\mathcal{I}}_i$, within the image set captures an axial slice of the prostate (in the ${}^{P}x$-${}^{P}y$ plane). The image set is described in terms of its placement and size in the patient attached frame, $\{P\}$, with distances measured using millimetres. Each axial image slice is taken at a depth ${}^P{z}_i$, where, for the patient image sets used in this work, the slices are spaced at 5 mm intervals. The placement of the patient-attached coordinate frame $\{P\}$ is chosen such that the first image in the set ${}^P{\mathcal{I}}_0$ is located at ${}^P{z}_0=0$. When using axial US images to create the surgical plan, planning software allows the clinician to select points in each axial TRUS image which correspond to the edge of the prostate, thereby contouring the prostate. This prostate contour can, in general, be represented by a polygon that comprises a number of vertices, $p_{(v,i)}={[{}^P{x}_{(v,i)},{}^P{y}_{(v,i)},{}^P{z}_{(v,i)}]}^{\top }$, with v being an index for the vertices, in the axial image slice i (from the index of the TRUS image ${}^P{\mathcal{I}}_i$). For each axial slice, the polygon can be represented by the vectors

$$\begin{array}{cc}\hfill {}^P{\overrightarrow{x}}_i& =[{}^P{x}_{(0,i)},{}^P{x}_{(1,i)},...,{}^P{x}_{(n_v,i)}]\hfill \\ \hfill {}^P{\overrightarrow{y}}_i& =[{}^P{y}_{(0,i)},{}^P{y}_{(1,i)},...,{}^P{y}_{(n_v,i)}]\hfill \\ \hfill {}^P{\overrightarrow{z}}_i& =[{}^P{z}_{(0,i)},{}^P{z}_{(1,i)},...,{}^P{z}_{(n_v,i)}],\hfill \end{array}$$

(1)

where $n_v$ is the number of the vertices in the polygon. The position point vector is defined in terms of the patient-attached frame, $\{P\}$, and can be converted, for a particular image slice i, to the pixel locations in the US image frame, denoted by $\{U\}$, where the full set of image slices taken for a particular patient is given by ${}^U\mathbf{I}=\left\{{}^U{\mathcal{I}}_0,{}^U{\mathcal{I}}_1,...,{}^U{\mathcal{I}}_{nImg}\right\}$, and consists of a number, $nImg$, of image slices taken at depths ${}^{P}d=\left\{0\mathrm{mm},5\mathrm{mm},...,nImg·5\mathrm{mm}\right\}$, where ${}^P{z}_i$ is the i-th element of ${}^{P}d$.

To smooth out the point-based contour, we propose a model which describes the shape of the prostate in an axial TRUS image. For each image slice i, the prostate contour model, $R_i({\theta }_i)$, describes the radius of a curve, within the cylindrical coordinate system, as a weighted sum-of-sines (a finite-length Fourier series). The general form of the model is given by

$$R_i(\theta )={\alpha }_0+\sum _{j=1}^{n_{\omega }}{\alpha }_{j}sin(j\theta +{\varphi }_j),$$

(2)

where $n_{\omega }$ is the number of sine components in the model, j controls the frequency of sine components in the model ($j=[1,2,...,n_{\omega }]$), ${\alpha }_0$ is a positive constant (${\alpha }_0>0$), the ${\alpha }_j$ terms weight the magnitude of each of the sine functions, and the ${\varphi }_j$ terms are phase offsets for the sine functions. The model consists of a number of modes, where the j-th mode refers either to the constant, ${\alpha }_0$ for mode 0, or weighted sine component of the model (${\alpha }_{j}sin(j\theta +{\varphi }_j)$). In the patient frame, the model describes the curve

$$\begin{array}{cc}\hfill {}^P{x}_i& =R_i(\theta )·cos(\theta )+{}^P{x}_c\hfill \\ \hfill {}^P{y}_i& =R_i(\theta )·sin(\theta )+{}^P{y}_c\hfill \\ \hfill {}^P{z}_i& ={}^P{z}_i,\hfill \end{array}$$

(3)

for a particular image slice i, and the model can thus be seen as imposing a $\theta$ dependence onto the radius of a planar curve. The output value of $R_i(\theta )$ must be strictly positive over the domain of theta, where $R_i(\theta )>0,\forall -\pi \le \theta <\pi$. This sliced-based contour model is designed to incorporate physical and geometric attributes of the prostate in axial images, by enforcing the constraint that $R_i(\theta )$ is strictly positive ensures that the contour model is both smooth, continuous, and has no self-intersections in the patient-attached frame.

3.1. 3D Prostate Contour Model Fitting

To evaluate the proposed contour model, a constrained least-squares optimization will be used to fit manually segmented prostate contour points to the model. While the model is defined with respect to the prostate in a single image slice, ${}^P{\mathcal{I}}_i$, the model fitting optimization incorporates information from the entire patient dataset, where all $nImg$ slice-based contours in ${}^P\mathbf{I}$ will be used. Each slice’s contour will be represented through a model, $R_i(\theta )$, with the parameters of the model being estimated for that particular slice. By evaluating the goodness-of-fit for all of the slice models across the patient image set, an optimal, patient-specific, centre point and rotation offset for the model will found.

Contour points of the prostate, segmented by the clinician, are given in the point vector form described by (1). For contour fitting, a rotation offset ${}^P{\varphi }_c$ is considered. This rotation offset aims to reduce any fitting artefacts caused by the chosen constant ${\varphi }_j$ values in the model and allows multiple prostate contours to be easily compared, which will be important when creating the statistical shape model. As the model is formulated in cylindrical coordinates, the prostate contour points, for image slice i, are defined in terms of the patient-attached frame and must be transformed by

$$\begin{array}{cc}\hfill {\overrightarrow{\theta }}_i& =\mathrm{atan}2({}^P{\overrightarrow{y}}_i-{}^P{y}_c,{}^P{\overrightarrow{x}}_i-{}^P{x}_c)-{}^P{\varphi }_c\hfill \\ \hfill {\overrightarrow{r}}_i& =\sqrt{{({}^P{\overrightarrow{y}}_i-{}^P{y}_c)}^2+{({}^P{\overrightarrow{x}}_i-{}^P{x}_c)}^2}\hfill \\ \hfill {\overrightarrow{z}}_i& ={}^P{\overrightarrow{z}}_i,\hfill \end{array}$$

(4)

with ${}^P{\overrightarrow{x}}_i$,${}^P{\overrightarrow{x}}_i$, and ${}^P{\overrightarrow{z}}_i$ being the clinician-segmented points; the values of ${}^P{x}_c,{}^P{y}_c$, and ${}^P{\varphi }_c$ come from the augmented centre point vector, and ${\overrightarrow{\theta }}_i$, ${\overrightarrow{r}}_i$, and ${\overrightarrow{z}}_i$ are the transformed contour points. The individual elements of the transformed contour point vectors are given by

$$\begin{array}{cc}\hfill {\overrightarrow{\theta }}_i& =[{\theta }_{(0,i)},{\theta }_{(1,i)},...,{\theta }_{(n_v,i)}]\hfill \\ \hfill {\overrightarrow{r}}_i& =[r_{(0,i)},r_{(1,i)},...,r_{(n_v,i)}]\hfill \\ \hfill {\overrightarrow{z}}_i& =[z_{(0,i)},z_{(1,i)},...,z_{(n_v,i)}],\hfill \end{array}$$

(5)

with each $\left\{{\theta }_{(v,i)},r_{(v,i)},z_{(v,i)}\right\}$ set of points being mapped from the set $\left\{{}^P{x}_{(v,i)},{}^P{y}_{(v,i)},{}^P{z}_{(v,i)}\right\}$ of contour points.

The parameters of the model are denoted by ${\mathit{\beta }}_i=[{\alpha }_0,{\alpha }_1,...,{\alpha }_{n_{\omega }}]$ and are found through the least squares optimization

$$\begin{array}{c}\hfill \underset{{\mathit{\beta }}_i}{argmin}∥{\overrightarrow{r}}_i-R_i({\overrightarrow{\theta }}_i|{\mathit{\beta }}_i)∥\end{array}$$

(6)

where $R_i({\overrightarrow{\theta }}_i|{\mathit{\beta }}_i)$ is the model, incorporating parameters ${\mathit{\beta }}_i$, and the optimization is constrained, such that

$$R_i({\overrightarrow{\theta }}_i|{\mathit{\beta }}_i)>0,\forall -\pi \le \theta <\pi ,$$

(7)

holds for the fit ${\mathit{\beta }}_i$ parameters. A second least-squares optimization will be performed in order to find the optimal augmented centre point, ${}^{P}c={[{}^P{x}_c,{}^P{y}_c,{}^P{\varphi }_c]}^{\top }$, for fitting the 3D model. The augmented centre point corresponds to the centre point of the prostate for each patient. Having an optimal ${}^{P}c$ allows for the prostate shape model to be found directly, as finding ${}^{P}c$ removes the relative translation and rotation between two prostate contours in different image sets. Here, the parameters to be optimized are ${\mathit{\beta }}_c=[{}^P{x}_c,{}^P{y}_c,{}^P{\varphi }_c]$. The least-squares optimization, under the constraint given in (7), is then

$$\begin{array}{c}\hfill \underset{{\mathit{\beta }}_c}{argmin}\sum _{i=0}^{nImg}\underset{{\mathit{\beta }}_i}{argmin}∥\left\{{\overrightarrow{r}}_i|{\mathit{\beta }}_c\right\}-R_i(\left\{{\overrightarrow{\theta }}_i|{\mathit{\beta }}_c\right\}|{\mathit{\beta }}_i)∥\end{array}$$

(8)

where, $\left\{{\overrightarrow{r}}_i|{\mathit{\beta }}_c\right\}$ and $\left\{{\overrightarrow{\theta }}_i|{\mathit{\beta }}_c\right\}$ are the transformed points calculated from (4). Thus, the optimal augmented centre point is found, which minimizes the slice contour model error across all slices, i, in the image set ($i=[0,1,...,nImg]$). This fitted contour model, incorporating the prostate contour information from all 2D slices, thus describes the 3D contour of the prostate. The resulting 3D contour model, after fitting to all of the patient contour slices, is shown in Figure 3, where input contour points are the red dots and the blue curves are the resulting prostate contour model.

To evaluate the accuracy of the prostate contour model, the mean-squared difference and maximum absolute error between the output model and the input contour points will be used. The mean-squared error is given as

$$\mathrm{MSE}=\frac{1}{n_{pt}}\sum _{i=0}^{nImg}\left\{{\overrightarrow{r}}_i|{\mathit{\beta }}_c\right\}-R_i(\left\{{\overrightarrow{\theta }}_i|{\mathit{\beta }}_c\right\}|{\mathit{\beta }}_i),$$

(9)

where $n_{pt}$ is the total number of clinician-segmented contour points that the model was fit to. The maximum absolute error is given as

$$\mathrm{MAE}=\underset{i}{max}\underset{v}{max}\left|\left\{{\overrightarrow{r}}_i|{\mathit{\beta }}_c\right\}-R_i(\left\{{\overrightarrow{\theta }}_i|{\mathit{\beta }}_c\right\}|{\mathit{\beta }}_i)\right|,$$

(10)

with v being the index of a vertex point within the input slice contour points; hence the MAE is the maximum absolute error between the model and all of the input contour points. These two metrics were evaluated for a set of 100 3D clinician-segmented contours (consisting of 7 to 13 2D slice contours each) based on a TRUS images of the prostate and a set of 100 3D clinician-segmented contours (consisting of 7 to 16 2D contours each) based on MR images of the prostate. While this work is focused on US prostate segmentation, the second MR dataset was evaluated to see if the model had any bias with respect to input imaging modality. In these results, and for the prostate contour after graph cutting segmentation, only four modes were used in the model ($n_{\omega }=4$), and the constant phase offset values, ${\varphi }_j$, were chosen to be ${\varphi }_1=0,{\varphi }_2=\pi /2,{\varphi }_3=0,$ and ${\varphi }_4=\pi /2$. Therefore, for each contour slice, this four mode model is given by

$$R_i(\theta )={\alpha }_0+{\alpha }_{1}sin(\theta )+{\alpha }_{2}cos(2\theta )+{\alpha }_{3}sin(3\theta )+{\alpha }_{4}cos(4\theta )$$

(11)

where the parameters to be fit for slice i are ${\mathit{\beta }}_i=[{\alpha }_0,{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4]$. The mean value and standard deviation of the MSE and MAE across the 100 contours and the range (minimum-to-maximum) of the MAE for the two datasets is presented in

Table 1. Prostate model fitting results.

Data Set	Average MSE (mm)	Average Max Abs Error (mm)	Range of Max Abs Error (mm)
US Contours	0.611 (±0.165)	1.519 (±0.310)	0.679–2.031
MR Contours	0.616 (±0.167)	1.565 (±0.305)	0.794–2.124

. The results show that the contour model has an average sub-millimetre mean-square error accuracy and a mean maximum absolute error less than 1.6 mm, regardless of imaging modality. The results also indicate that, by using only four modes, the clinician-segmented contours are modelled with acceptable accuracy, where the largest discrepancy (maximum absolute error) between the input clinician contour point and the contour model across the two datasets was under 2.2 mm.

3.2. Statistical Prostate Shape Model

For the prostate segmentation algorithm, a 3D probabilistic shape model of the prostate will be created. This shape model is used to find the likelihood of the edge of the prostate being at a particular point. To generate the prostate shape model, a statistical distribution of the prostate contour locations is found by analysis of the prostate shape information contained in a number of 3D contour models from different patients (image sets). The resulting shape model characterizes the average prostate shape, and variance of the prostate shape, measured across a large number of 3D contour models.

After fitting in (8), a particular 3D contour model, given an index m, will be denoted by

$$R(\theta ,i|{\mathit{\beta }}_m)=R_i(\theta |{\mathit{\beta }}_i)$$

(12)

where ${\mathit{\beta }}_m$ contains the model parameters of all the slices, $m_{nImg}$ for a particular prostate, such that ${\mathit{\beta }}_m=[{\mathit{\beta }}_0,{\mathit{\beta }}_i,...,{\mathit{\beta }}_{m_{nImg}}]$.

The probabilistic prostate shape model will incorporate many, $n_{Model}$, 3D contour models from different patient TRUS images, and is created using linear interpolation of Gaussian distributions. The radius R, from the 3D contour models, will be evaluated at a discrete number of ${\theta }_c$ points to build the probability model for each slice. The probabilistic prostate shape model will then be in the same slice-based 3D form as the contour model, in the cylindrical coordinate system. Each Gaussian component of the shape model will be centred at $\overrightarrow{\theta }=[{\theta }_1,{\theta }_2,...,{\theta }_{n_{\theta }}]$ points with the vector

$$\overrightarrow{\theta }=\left[-\pi ,-\pi +2\pi \frac{1}{n_{\theta }-1},-\pi +2\pi \frac{2}{\theta -1},...,-\pi +2\pi \frac{n_{\theta }-1}{n_{\theta }-1}\right]$$

(13)

giving $n_{\theta }$ points sampled at regular intervals in the range $(-\pi ,\pi ]$. The Gaussian normal distribution for a single point, ${\theta }_c$, in contour slice i, is

$$p_c(r|{\theta }_c,i)=\frac{1}{w_{(c,i)}\sqrt{2\pi {\sigma }_{(c,i)}^2}}exp\left(-\frac{{(r-{\mu }_{(c,i)})}^2}{2{\sigma }_{(c,i)}^2}\right),$$

(14)

where $w_{c,i}$ is a normalization constant, ${\mu }_{(c,i)}$ is the mean radius, and ${\sigma }_{(c,i)}$ is the standard deviation for all of the shape model at the ${\theta }_c$ point. Note that the ${\theta }_c$ points are used for all slices of the shape model. The CDF of this distribution is given by

$$\begin{array}{c}{P}_c(R\le r|\theta ,i)={\int }_0^{r}p(r|\theta ,i)dr-{\int }_{-\infty }^{0}p(r|\theta ,i)dr\\ P_c(R\le r|\theta ,i)=\frac{1}{2}\left[\mathrm{erf}\left(\frac{(r-{\mu }_{(c,i)})}{{\sigma }_{(c,i)}\sqrt{2}}\right)-\mathrm{erf}\left(\frac{(0-{\mu }_{(c,i)})}{{\sigma }_{(c,i)}\sqrt{2}}\right)\right]\end{array}$$

(15)

where $\mathrm{erf}()$ is the error function, and the normalization constant is given by

$$w_{(c,i)}=\frac{1}{2}\left[1-\mathrm{erf}\left(\frac{(0-{\mu }_{(c,i)})}{{\sigma }_{(c,i)}\sqrt{2}}\right)\right]$$

(16)

such that the probability of $P_c(0\le R\le \infty |{\theta }_c,i)=1$ as desired. With this normalization constant, the mean and standard deviation of the radius at ${\theta }_c$, across all $n_{Model}$ contour models, are calculated by

$${\mu }_{(c,i)}=\frac{1}{n_{Model}}\sum _{m=1}^{n_{Model}}R({\theta }_c,i|{\mathit{\beta }}_m)$$

(17)

and

$${\sigma }_{(c,i)}=\frac{1}{n_{Model}+1}\sum _{m=1}^{n_{Model}}{\left(R({\theta }_c,i|{\mathit{\beta }}_m)-{\mu }_{(c,i)}\right)}^2$$

(18)

in the same manner as a standard Gaussian distribution.

To find the contour edge probability for any value of $\theta$ in the slice, linear interpolation between the two closest ${\theta }_c$ distributions is used. For a point $\theta$, let ${\theta }_c^{-}$ be the closest ${\theta }_c$ point that is less than $\theta$ and ${\theta }_c^{+}$ be the closest ${\theta }_c$ which is greater than $\theta$. If $\theta$ is between the largest ${\theta }_c$ and $\pi$ (${\theta }_{n_{\theta }}<\theta \le \pi$), then ${\theta }_c^{-}={\theta }_{n_{\theta }}$ and ${\theta }_c^{+}=-\pi$, so that the contour probabilities at $-\pi$ and $\pi$ are the same. The interpolated probability for all $\theta$ values is then

$$p(r|\theta ,i)=\left(\frac{\theta -{\theta }_c^{-}}{n_{\theta }-1}\right)p_c(r|{\theta }_c^{-},i)+\left(1-\frac{\theta -{\theta }_c^{-}}{n_{\theta }-1}\right)p_c(r|{\theta }_c^{+},i)$$

(19)

where normalization of (14) means that the interpolated shape model value is also a normalized (standard) probability. The CDF of the probability for the shape model for any point $\theta$ is given by

$$\begin{array}{c}P(R\le r|\theta ,i)=\left(\frac{\theta -{\theta }_c^{-}}{n_{\theta }-1}\right)P_c(R\le r|{\theta }_c^{-},i)+\\ \left(1-\frac{\theta -{\theta }_c^{-}}{n_{\theta }-1}\right)P_c(R\le r|{\theta }_c^{+},i)\end{array}$$

(20)

where $R\le r$ is the probability that the edge of the prostate is located at some point less than r for a given value of $\theta$. The resulting mean and standard deviation of the probabilistic prostate shape model, from a set of 100 prostate contours segmented by clinicians from TRUS images, is shown in Figure 4.

When creating the shape model, the number of contour slices between 3D models may vary; $nIm{g}_0$ may or may not equal $nIm{g}_1$, for instance. The number of slices, $nSlice$, in the prostate shape model will be limited to the smallest number of slices, $nImg$, across all of the 3D contour models. To line up the slices, the prostate shape model components are calculated for all models with $nSlice$ images. For any contour model with more than $nSlice$ images, $nSlice$ components of the 3D contour model are compared to the mean prostate shape components, such that

$$MD(i_z)=\sum _{i=0}^{nSlice}\sum _{c=1}^{n_{\theta }}\left|R({\theta }_c,i+i_z|{\mathit{\beta }}_m)-{\mu }_{(c,i)}\right|$$

(21)

with $i_z$ being a slice index offset. The model difference is evaluated starting at index $i_z=0$, the first $nSlice$ components are compared to the mean 3D contour in the prostate shape model, then the next $nSlice$ components, $i_z=1$, are evaluated, continuing until $i_z=nIm{g}_m-nSlice$. The set of slices with the lowest model difference, starting at $i_z$, will be incorporated into the prostate shape model by updating ${\mu }_{(c,i)}$ and ${\sigma }_{(c,i)}$. Essentially, the $nSlice$ components from each 3D contour model which best match the prostate shape model are used.

4. Superpixel Algorithm Incorporating Image Statistics

In this section, we propose a superpixel algorithm which incorporates statistical information from an image to detect edges and salient features. Superpixel algorithms, such as the SLIC algorithm [31], group individual pixels within an image into large clusters (or superpixels). An iterative method is used to update which pixels are grouped into each cluster to minimize a score or distance function. This work will follow a similar design for pixel-based clustering and proposes an additional optimization (or clustering) layer that is built on top of, and groups together, the pixel-based clusters. The regions in the input image delineated by this second layer will be considered the output superpixels of the algorithm, and will be used as the basis for the prostate contouring and segmentation. The algorithm updates and optimizes the pixel clusters and superpixel regions in an iterative manner, where in each iteration the pixel clusters will first be updated (Section 4.1) and then the superpixel regions are updated (Section 4.2).

4.1. Pixel Clustering

At the pixel level, superpixel clustering will be performed on an input image $\mathcal{I}(i)$, where ${\mathcal{I}}_{width}$ and ${\mathcal{I}}_{height}$ are the width and height of the image. An iterative method will be used to group individual pixels, i, in $\mathcal{I}(i)$ into clusters $C_j$. A pixel may only belong to a single cluster, and from the work of [31], a label image, $\mathcal{L}(i)$, is created to indicate which cluster the pixel is grouped into. The label image $\mathcal{L}(i)$ is identical in size (width and height) to the input image $\mathcal{I}(i)$, where $\mathcal{L}(i)=j$ signifies that pixel i is a member of cluster $C_j$. The pixel-based clustering will be described in two separate stages, the initialization stage, where initial pixel clusters are created from the input image, and the update stage, which is performed iteratively to optimize which cluster a particular pixel belongs to.

At initialization, the label image is divided into a regular grid of squares (or approximately a regular grid for a non-square input image). The chosen height and width of each grid, $s_{grid}$, breaks the input image into $n_s$ square sections, where $n_s=\frac{{\mathcal{I}}_{width}·{\mathcal{I}}_{height}}{s_{grid}^2}$ (or $n_s\le \frac{{\mathcal{I}}_{width}·{\mathcal{I}}_{height}}{c_{grid}^2}$ when the input image is non-square). To assign the pixels to an initial set of clusters, each grid square is given a unique index, j, where $j=\{1,2,\dots ,n_s\}$. The label for all of the pixels within a grid square is then updated, $\mathcal{L}(i)=j$, to match.

With the label image created and the initial pixel clustering performed, an array of cluster elements, $\mathbf{C}$, will be created. The array takes the form

$$\mathbf{C}={[C_0,C_1,...,C_{n_c}]}^{\top },$$

(22)

where each element is a pixel cluster, $C_j$, containing a set of parameters to be calculated. The cluster parameters are expressed as

$$C_j=[{\mu }_j,x_j,y_j],$$

(23)

with $(x_j,y_j)$ being the centroid of the cluster and ${\mu }_j$ denoting the average intensity of the pixels in the cluster. The values of ${\mu }_j$,$x_j$, and $y_j$ are computed through

$$\begin{array}{cc}\hfill {\mu }_j& =\sum _{i\in C_j}\frac{\mathcal{I}(i)}{n_{jpx}}\hfill \\ \hfill x_j& =\sum _{i\in C_j}\frac{x_i}{n_{jpx}}\hfill \\ \hfill y_j& =\sum _{i\in C_j}\frac{y_i}{n_{jpx}},\hfill \end{array}$$

(24)

where $x_i$ and $y_i$ are the x–y coordinates of the pixel i and $n_{jpx}$ is the number of pixels $i\in C_j$ (the set of all pixels i in the cluster $C_j$).

The distance of an individual pixel i with respect to the cluster to which it belongs, $C_j$, will be stored in a distance image, $\mathcal{D}(i)$, which is of the same size as the input image, $\mathcal{I}(i)$, and which will be used when updating the pixel label image $\mathcal{L}(i)$. Before clustering, preprocessing was used to restrict the area of the US image to a region of interest (ROI); see Appendix A. To ignore regions in the US image which contain invalid pixels, for instance, those pixels in the warped US image that are outside the desired ROI ($i\notin {\mathcal{I}}_{ROI}$, see Figure A1b), the distance image will be initialized to allow only valid pixels to be considered during cluster updating. The pixel value in the distance image will be set to infinity (or a large positive value) for pixels inside of the ROI and $-1$ for pixels outside of the ROI, such that

$$\mathcal{D}(i)=\left\{\begin{array}{cc}\infty \hfill & \mathrm{for}i\in {\mathcal{I}}_{ROI}\hfill \\ -1\hfill & \mathrm{for}i\notin {\mathcal{I}}_{ROI}\hfill \end{array}\right.$$

(25)

with $i\in {\mathcal{I}}_{ROI}$ and $i\notin {\mathcal{I}}_{ROI}$ being the pixels inside of and outside of the ROI, respectively. The invalid pixels are removed from their respective clusters through modifying the label image, such that $\mathcal{L}(i)=0,\forall \{i\notin {\mathcal{I}}_{ROI}\}$. If there are large areas outside of the ROI then there may be several clusters which no longer contain any pixels. In this case, the clusters without any pixels are removed from the cluster array so that $\mathbf{C}$ only contains valid pixel clusters.

4.2. Superpixel Regions

The input image pixels are now grouped into SLIC style clusters, as proposed in [31], through the creation of the label image. We now propose our extension to the SLIC superpixel algorithm with a second optimization layer. This additional layer is the mechanism through which image statistics will be incorporated into the superpixel algorithm. The underlying implementation of the second layer incorporates a graph component, which will allow for the output of the superpixel algorithm to be segmented by graph cutting in a straightforward manner. The components of this second layer will be referred to as superpixel regions, which will be stored in a list $\mathbf{R}$, given by

$$\mathbf{R}=[R_0,R_1,...,R_n]{,}^{\top }$$

(26)

with $R_k$ denoting the individual regions. In a similar manner to a pixel cluster, $C_j$, being composed of a number of individual pixels, a region $R_k$ is composed of a number of pixel clusters ($C_j$), thus building a second optimization layer on top of the pixel clusters. This second layer creates a pyramid-type hierarchy in this superpixel implementation, as diagrammed in Figure 5.

Each region consists of a series of calculated values, given by

$$R_k=[{\widehat{\mathbf{F}}}_{\mathbf{k}},x_k,y_k,j_k^{★}],$$

(27)

where $j_k^{★}$ is the list of indexes to the pixel clusters ($C_j$) which make up the superpixel region, $(x_k,y_k)$ is the centroid of the region, and ${\widehat{\mathbf{F}}}_{\mathbf{k}}$ is the empirical cumulative density function (CDF) of all of the pixels contained within the region. As a region is made up of a number of pixel clusters, all of the pixels belonging to the region’s subordinate clusters are defined to be pixels of the region. Formally, this is described as the union of all of the pixels in the component clusters of $R_k$, given as

$$i\in R_k=\bigcup _{j\in k}i\in C_j,$$

(28)

where $j\in R_k$ are the component clusters of the region. The centroid values of the region, $x_k$ and $y_k$, are calculated as follows:

$$\begin{array}{cc}\hfill n_{kpx}& =\sum _{j\in R_k}{n}_{jpx}\hfill \\ \hfill x_k& =\sum _{j\in R_k}\frac{n_{jpx}·x_j}{n_{kpx}}\hfill \\ \hfill y_k& =\sum _{j\in R_k}\frac{n_{jpx}·y_j}{n_{kpx}},\hfill \end{array}$$

(29)

with $n_{jpx}$ being the number of pixels in a cluster and $n_{kpx}$ being the total number of pixels within the region. The empirical CDF, ${\widehat{\mathbf{F}}}_{\mathbf{k}}$, allows for a metric based on pixel intensity (image) statistics within a region; the CDF is calculated through

$$\begin{array}{c}{\widehat{F}}_k(x)=\mathrm{Count}\left(\mathcal{L}(i\in R_k)\le x\right)\\ {\widehat{\mathbf{F}}}_{\mathbf{k}}=[{\widehat{F}}_k(0),{\widehat{F}}_k(1),...,{\widehat{F}}_k(255)],\end{array}$$

(30)

given a counting function, $\mathrm{Count}()$, which returns the number of pixels within a set, in this case, the set of pixels in the region $R_k$ with an intensity value, $\mathcal{L}(i)$, less than I. This is essentially the cumulative summation of all the elements in the image intensity histogram, for all of the pixels in the region $i\in R_k$.

The superpixel region optimization uses a graph-based component, $\mathcal{G}$, to keep track of which regions are adjacent to one another, called neighbours. For this work, an undirected graph of the form

$$\mathcal{G}=\{\mathcal{V},\mathcal{E}\},$$

(31)

is used. Here, $\mathcal{V}$ are the vertices of the graph which are interconnected by edges $\mathcal{E}$. The vertices of the graph will be the superpixel regions, and edges define which regions are neighbours to one another; the vertices and edges have the form

$$\begin{array}{cc}\hfill \mathcal{V}& =\{0,1,...n_R\}\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \mathcal{E}& =\{k_{v1}^{★},k_{v2}^{★}\},\hfill \end{array}$$

(33)

with $\mathcal{V}$ being a set of all of the region indexes (k). The edges in the graph $\mathcal{E}$ are described by a set of double-valued elements, $k_{v1}$ and $k_{v2}$, with an individual line edge element connecting two vertices given by ${\mathit{E}}_v=\left(k_{v1},k_{v2}\right)$, where $k_{v1}$ and $k_{v2}$ are the indexes of two neighbouring regions $R_{k1}$ and $R_{k2}$. With the graph being undirected, the ordering of the elements in an edge is ignored, such that $\left(k_{v1},k_{v2}\right)\equiv \left(k_{v2},k_{v1}\right)$. Additionally, under this independence from element ordering, the connection between two vertices $k_{v1}$ and $k_{v2}$ is described by a single edge without duplication. For example, only one of the two edge elements ${\mathit{E}}_a=\left(k_{v1},k_{v2}\right)$ and ${\mathit{E}}_b=\left(k_{v2},k_{v1}\right)$ are allowed to be a member of the set $\mathcal{E}$.

The adjacency graph is used to definite a region’s neighbours and will be used when optimizing region placement and size. The neighbourhood of a region will also be used as part of the weighting to determine the optimal prostate segmentation graph-cut. The neighbours to a region, $R_k$, are defined as a set of all of the regions, $R_m$, which have a graph edge connected to $R_k$. The set of neighbours, ${\mathcal{N}}_k$, is derived from the elements ${\mathit{E}}_v=\left(k_{v1},k_{v2}\right)$ in the edge set $\mathcal{E}$, where ${\mathcal{N}}_k$ is given by

$$\begin{array}{c}{\mathit{N}}_{k1}=\{k_{v1},\forall {\mathit{E}}_v\mathrm{where}{k}_{v2}=k\}\\ {\mathit{N}}_{k2}=\{k_{v2},\forall {\mathit{E}}_v\mathrm{where}{k}_{v1}=k\}\\ {\mathcal{N}}_k={\mathit{N}}_{k1}\cup {\mathit{N}}_{k2}\end{array}$$

(34)

such that ${\mathcal{N}}_k$ contains all of the $k_{v1}$ and $k_{v2}$ regions that are paired with (connected to) k in any edges of the graph $\mathcal{E}$.

To initialize the regions and the graph, a region $R_k$ is created for each cluster $C_j$ such that $k=j$, resulting in an initial one-to-one mapping between the clusters and regions. As pixel clusters are created in a grid pattern (see Section 4.1 and Figure 6a), the graph $\mathcal{G}$ can be easily created by finding the (four-connected) grid squares that adjoin each other; see Figure 6b,c. With the pixel clusters, superpixel regions, and graph structure created, the algorithm will then work to reduce the number of superpixel regions while simultaneously enlarging the area (number of pixels) contained within those regions. The relationship between the individual pixels, the pixel clusters, and superpixel regions, along with the graph structure, are shown in Figure 5. Explicitly, in Figure 5, $R_1$ is a superpixel region comprising three pixel clusters, $C_1,C_2,$ and $C_3$, with neighbouring regions $R_2$ and $R_3$ connected via graph edges.

4.3. Superpixel Algorithm Outline

After both the pixel clusters and superpixel regions have been initialized, the algorithm will step through a series of update steps, such that large superpixel regions will be created. These update steps work to reduce the number of regions in an image by increasing the area (number of pixels) belonging to each region. Using the region layer and its associated graph, the algorithm works to decrease the number of regions by joining together statistically similar regions, returning a smaller number of larger regions each process iteration. This statistics-based approach creates regions from which edges and salient features can be found in a straightforward manner.

Each update step can be broken into two sub-steps, the first of which works at the individual pixel level to determine which pixels belong to a particular cluster, acting to increase or decrease the size of the clusters with each step. The second sub-step controls the growth of the superpixel regions by increasing the number of clusters belonging to each region and grouping together pixel clusters to maximize region size. This region update step incorporates the region statistics to only group together regions which have comparable image statistics, thus separating regions which vary in average pixel intensity and standard deviation (or covariance) of pixel intensity. This statistical separation uses the premise developed in [32], where edge and corner features in a generic image have been shown to have distinct statistical properties from smooth features in the image. During the region update step, a hard limit is placed on how dissimilar two regions can be when joined together, with small regions being retained if they are distinct enough from neighbouring regions.

The algorithm uses a few parameters which can be tuned to control how the pixel clusters and superpixel regions are modified during updating. The use of these parameters will be described in full in the text for the cluster update and region update sub-steps. One of the parameters of the algorithm is a weighting vector, $\mathit{\lambda }$, used to determine the goodness-of-fit for a particular pixel to a cluster and to limit region-to-region statistical similarity. The values within this vector are

$$\mathit{\lambda }=\left[{\lambda }_{\mu },{\lambda }_{L2},{\lambda }_{KS}\right],$$

(35)

with ${\lambda }_{\mu }$ and ${\lambda }_{L2}$ being used in the cluster update step to weight a distance function (score) between a pixel and nearby clusters. The value ${\lambda }_{KS}$ will be described in the region update step and defines the limit on how dissimilar two regions must be before they cannot be merged.

There are a few constraints on the output superpixels of the algorithm which will either be enforced by the algorithm or be checked at the end of each iteration to determine when the optimization is complete (to be discussed below). The desired size of the output superpixel regions is defined as $r_{darea}$ with constraints on the minimum pixel cluster and region area, $c_{min}$ and $r_{min}$, respectively. The minimum acceptable size for the regions is chosen to be equivalent to or larger than the size of the initial grid squares such that $r_{min}\ge s_{grid}$ with the clusters being required to be larger than $c_{min}$, where $c_{min}={\alpha }_c{r}_{min}$ with ${\alpha }_c$ being a scaling constant taking values within the range $0<{\alpha }_c\le 1$. For this work, a value of ${\alpha }_c=0.75$ was found to be sufficient. A limit is also placed on the maximum number of pixel clusters, $c_{max}$, allowed within each region. The desired number of regions returned by the algorithm is controlled by the parameter $r_{dnum}$, where, nominally, there will be $r_{dnum}$ or fewer regions when the algorithm terminates.

4.3.1. Cluster Update

Two distance functions are used during the pixel cluster update operations which score which cluster, $C_j$, an individual pixel, i, should belong to. The calculation of this score is based on two different distance metrics, $D_{\mu }(i,j)$ and $D_{L2}(i,j)$. The metric $D_{\mu }(i,j)$ is the absolute difference between the pixel’s intensity (in the input image $\mathcal{I}(i)$) with respect to the average pixel intensity, ${\mu }_j$, of a cluster. The other metric $D_{L2}(i,j)$ is defined as the Euclidean distance from the pixel’s coordinates, $(x_i,y_i)$, to the centroid, $(x_j,y_j)$, of the cluster. A distance score, $D_{Score}(i,j)$, for a pixel is calculated as

$$\begin{array}{c}{D}_{\mu }(i,j)=\left|{\mu }_i-\mathcal{I}(i)\right|\\ D_{L2}(i,j)=\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2}\\ D_{Score}(i,j)=\left(\frac{D_{\mu }(i,j)}{{\lambda }_{\mu }}\right)+\left(\frac{D_{L2}(i,j)}{{\lambda }_{L2}}\right)\end{array}$$

(36)

with ${\lambda }_{\mu }$ and ${\lambda }_{L2}$ being two tunable parameters that weight the importance of the average intensity distance and Euclidean distance. The distance value of each pixel, with respect to the cluster to which it currently belongs, is stored in the distance image $\mathcal{D}(i)$.

The goal of the cluster update step is to minimize the total of the distance score for all (valid) pixels in the image. Here, the method for updating the clusters is similar in most respects to the SLIC superpixel algorithm [31]. Minimization of the total distance score is achieved by moving a pixel to a different cluster if the pixel’s score, with respect to that new cluster, is less than its current score value $\mathcal{D}(i)$. To determine if a pixel should be moved, the distance score for each pixel, i, is calculated with respect to each cluster $C_j$. This calculation is evaluated individually for all of the clusters within the cluster array $\mathbf{C}$ in a loop. If, during the loop, the distance score, $D_{Score}(i,m)$, for a cluster, $C_m$, is less than the current value in $\mathcal{D}(i)$, then the pixel is relabelled through the function

$$\mathcal{L}(i)=m\mathrm{if}{D}_{Score}(i,m)<\mathcal{D}(i),$$

(37)

so that the pixel is now a member of the cluster $C_m$. To speed up processing, only the distance values for pixels within a square region (of size $2r_{desarea}\times 2r_{desarea}$) around each cluster’s centroid are evaluated, based on the suggested method in [31]. After all of the clusters have been iterated through, with pixels relabelled to minimize the distance scores, the values in the cluster array, $\mathbf{C}$, are updated using (24). Note that the distance image value for invalid pixels, those outside of the ROI, was set to $-1$ when the distance image was created. As the distance score is positive semi-definite ($D_{Score}(i,j)\ge 0$), invalid pixels will never be incorporated into a cluster and do not need to be considered separately.

4.3.2. Region Update

With the cluster update sub-step completed, the region update sub-step is then performed. Region updating incorporates pixel cluster information and the calculation of image statistics to modify the individual regions $R_k$ and the region adjacency graph $\mathcal{G}$. The goal of the region update sub-step is to join neighbouring regions together so that, at the termination of the superpixel algorithm, large, statistically similar regions in the input image will be delineated. In a similar manner to the pixel cluster update, a distance metric, $D_{KS}(k_1,k_2)$, will be used to inform the region update process. The distance metric used is the Kolmogorov–Smirnov (KS) test, which gives a distance between two empirical CDFs [33]. The KS test is the mechanism through which changes in image statistics in the input image will be incorporated into the region update process. Using the KS test, regions with similar statistical properties (or equivalently similar empirical CDFs) will be joined together and, conversely, dissimilar regions will remain separated, thus creating regions with distinct image statistics. The region update step consists of two separate optimizations, a region-statistic optimization through performing a region union operation, and a subsequent optimization of the pixel cluster through cluster merging.

Conceptually, measuring the statistical dissimilarity between superpixel regions is a form of texture discrimination and edge feature extraction. In [34], a model incorporating pixel cluster covariance was used to create a line and edge feature extractor for generic images. This concept was incorporated in the work of [32], where a line/edge detector was formulated based only on the gradient of the pixel region covariance (change between covariance in neighbouring pixel regions). Expanding on the use of covariance for edge extraction, ref. [35] showed that a statistical distance measurement, the Bhattacharyya space, allowed for the discrimination of various textures within an image. The Bhattacharyya space was applied to large-scale region of interest (ROI) finding and feature segmentation in [36]. Refs. [37,38] demonstrate that comparable ROI and feature segmentation can be performed on ultrasound images using speckle and image statistics.

Region Union:

The first operation, the region union, will incorporate the adjacency information stored in the graph $\mathcal{G}$ to evaluate the KS distance between neighbouring regions and to choose which regions a union operation should be performed on. As each region consists of a pixel cluster (or pixel clusters after the first iteration of the algorithm), the values in the region list $\mathbf{R}$ must be updated to incorporate the changes to the pixel clusters during their update step. Values of the centroid of each region, given in (29), and the empirical CDF of each region, given in (30), need to be calculated. Using these new region list values, the distance metric is computed for each edge ${\mathit{E}}_v=\left(k_{v1},k_{v2}\right)$ in the graph $({\mathit{E}}_v\in \mathcal{E})$. The edges describe the interconnection between the region elements, a region’s neighbours, and the score measures how statistically similar neighbouring regions are. The score is stored in an array,

$$\mathbf{S}={[s_0,s_1,...,s_{\ell }]}^{\top },$$

(38)

which contains the same number of entries as the edge structure, $\ell =v_{max}$ where $v_{max}$ is the largest index of edges ${\mathit{E}}_v\in \mathcal{E}$. The distance score is then calculated along each edge via the function

$$\begin{array}{c}{D}_{KS}(k_i,k_j)=\underset{x}{sup}\left|{\widehat{F}}_i(x)-{\widehat{F}}_j(x)\right|\\ s_v=D_{KS}(k_{v1},k_{v2})\end{array}$$

(39)

which returns the largest mismatch between the values of the empirical CDFs, ${\widehat{F}}_k(x)$ (stored in ${\widehat{\mathbf{F}}}_{\mathbf{k}}$), across all of the input pixel intensities, $x=\{0,1,...,255\}$, seen in the input image $\mathcal{I}(i)$. The KS test has several useful properties, in that it is symmetric $D_{KS}({\widehat{F}}_i(x),{\widehat{F}}_j(x))\equiv D_{KS}({\widehat{F}}_j(x),{\widehat{F}}_i(x))$ and does not require the data to be modelled as belonging to some underlying distribution. The symmetry property of the KS test allows for the distance values to be evaluated only once per graph edge irrespective of the order of elements in an edge ${\mathit{E}}_v$. The second property, the independence from an underlying distribution, is critical in the context of ultrasound image processing as the pixel intensities seen in B-mode images have been shown to follow a number of different statistical distributions varying with tissue type, location in the body, ultrasound intensity, and transducer focusing [39].

A low KS score between two regions indicates that they are statistically similar. In a loop, each region $R_k$ and its neighbours $R_m\in {\mathcal{N}}_k$ will be compared to see if they are statistically similar enough, and if so, a union operation will be performed. For each $R_k$, the neighbour with the minimum KS distance,

$$D_{KSmin}(R_k)=\underset{k_m}{min}{D}_{KS}(k,k_m)\forall k_m\in {\mathcal{N}}_k,$$

(40)

giving the index m of the neighbouring region $R_m$ with the minimum statistical difference $D_{KSmin}(R_k)$ to $R_k$ is identified. A binary threshold is then used to determine if a union will be performed:

$$R_k^{\prime }=R_k\cup R_i\mathrm{if}\left(D_{KSmin}(k)\le {\lambda }_{KS}\mathrm{or}{n}_{kpx}<r_{min}\right),$$

(41)

which results in a new region $R_k^{\prime }$ when $R_k$ and $R_m$ are either similar to one another, or forcing a union, when the area of $R_k$ is smaller than the minimum allowable region size $r_{min}$. When $n_{kpx}\ge r_{min}$, the threshold ${\lambda }_{KS}$ prevents unions between statistically dissimilar regions. If neither condition is met, then $R_k$ will be ignored and the next region in the set $\mathbf{R}$ will be tested.

When either of the the union conditions are met, the union results in the two regions being combined to create one larger region; $R_k^{\prime }=R_k\cup R_m$. This new region, $R_k^{\prime }$, takes ownership over all of the pixel clusters that originally belonged to $R_k$ and $R_m$. From (27), the region is described by its parameters, $R_k^{\prime }=[{\widehat{\mathbf{F}}}_{\mathbf{k}}^{\prime },x_k^{\prime },y_k^{\prime },j_k^{\prime ★}]$. The indexes for the all of the pixel clusters in this new region, $j_k^{\prime ★}$, are given by

$$j_k^{\prime ★}=j_k^{★}\cup j_m^{★},$$

(42)

where $j_k^{★}$ are indexes for the clusters $\{C_j,j\in R_k\}$ and, likewise, $j_m^{★}$ are indexes for the clusters $\{C_j,j\in R_m\}$. The values in the rest of $R_k$ are then updated, from (29) and (30), to incorporate all of the information from the pixel clusters in the new region. A diagram of the region union operation is given in Figure 7, where a union is performed on the two regions in Figure 7a (indicated in yellow), combining them into the single region in Figure 7b.

As shown in Figure 7, after the union operation, the region list R and the adjacency graph, $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$, need to be updated. For simplicity, the newly created region $R_k^{\prime }$ will replace the parameters of $R_k$ (at row k) in the region list, $\mathbf{R}$, and the row m, corresponding to $R_m$, is removed from $\mathbf{R}$. To update the graph, the m vertex ${\mathit{V}}_m$ is deleted and the edge ${\mathit{E}}_v=\left(k,m\right)$ is removed. Any additional edges, ${\mathit{E}}_v=\left(k_{v1},k_{v2}\right)$, which originally connected to m, $\{k_{v1},k_{v2}=m\}$, will be updated to connect to k. Any duplicate edges in $\mathcal{E}$ will also be cleaned up to ensure that there is only a single edge between vertices connected to one another.

This procedure can continue until all of the regions $R_k$ have been tested. From empirical observation, however, it is advantageous to limit the number of unions performed during one algorithm iteration. This allows more time for the pixel cluster layer to converge, where the base SLIC algorithm in [31] was shown to require a number of iterations before satisfactory cluster convergence. For this work, a subset of all of the regions is evaluated for union operations in each iteration. The subset of regions to analyze was chosen to encourage regions to grow in size and to guarantee that small regions, those with an area less than $r_{min}$, had a union performed on them. To do this, the size values $n_{kpx}$ in the superpixel region list $\mathbf{R}$ were sorted such that the smallest regions are investigated first. A limit, $n_{union}$, was placed on the number of unions performed at each iteration and, empirically, the superpixels were found to converge well when $n_{union}=50$ regions per iteration. Using the list of region indexes sorted according to size, all regions smaller than $r_{min}$ are evaluated. If removing those small regions does not exceed the $n_{union}$ unions-per-iteration limit, then the algorithm continues working through the sorted list to evaluate regions until $n_{union}$ union operations have been performed. This ensures all regions have an area larger that $r_{min}$ while achieving adequate convergence of the pixel clusters and growth of the superpixel regions.

Cluster Merge:

After the cluster update and region union operations, the optimization constraints for the pixel cluster size $c_{min}$ and number of clusters per region $c_{max}$ must be enforced. After the region union procedure, all regions smaller than $r_{min}$ are placed into nearby regions and the following cluster merge operation works on the pixel cluster components inside of the updated regions. Regions may also now contain more than the maximum allowable number of clusters, which will be corrected through the use of the cluster merge operation. In a similar manner to the region update method, the cluster merge method works to join together neighbouring pixel clusters to create new, larger, and contiguous pixel clusters. The cluster merge operation respects the pixel cluster–superpixel region hierarchy, such that a pixel cluster belongs only to a single region, and the merge operation does not move clusters between regions.

The cluster merge update loops through all of the individual pixel clusters, $C_j$, with the size of the cluster, $n_{jpx}$, being calculated. Any regions smaller than $c_{min}$ are merged with their nearest neighbouring clusters inside of their respective regions. The neighbours of a cluster $C_j$ belonging to a region $R_k$ are defined to be all of the other clusters within that region, $\{C_i\in R_k\ne C_j\}$. The nearest neighbour to the cluster $C_j$ is defined as the neighbouring cluster with the minimum Euclidean distance between their respective centroids,

$$N_N(C_j)=\underset{i}{min}\sqrt{{(x_j-x_i)}^2+{(x_j-x_i)}^2}\forall C_i\in R_k\ne C_j,$$

(43)

where $x_i,x_j,y_i,and{y}_k$ are pixel coordinates of the centroids, and $N_N(C_j)$ returns the index i of the nearest neighbour to $C_j$. The merge operation is performed by updating the label image $\mathcal{L}(i)$ such that,

$$\mathcal{L}(j)=i,\forall \mathrm{pixels}j\in C_j,$$

(44)

which relabels all of the pixels in the cluster $j\in C_j$ to belong to the cluster $C_i$. Lastly, the row j in the pixel cluster list $\mathcal{C}$ and the pointer $j_k^{★}=j$ in the parameters of region $R_k$ are then removed. Note that the scaling constant ${\alpha }_c$ between the minimum region size and minimum cluster size $c_{min}={\alpha }_c{r}_{min}$ is within the region $0<{\alpha }_c\le 1$. This ensures that, at the output of the region union operation, any region containing a cluster of size $n_{jpx}<c_{min}$ will necessarily have a second cluster that it can be merged with. An example of the cluster merge operation is shown in Figure 8, with the two clusters in Figure 8a (indicated in yellow) being combined into the single cluster in Figure 8b.

After all the pixel clusters smaller than $c_{min}$ have been merged, a second loop is used which calculates the number of cluster components $n_{ck}$ within a region $R_k$, where $n_{ck}=\mathrm{Count}(j_k^{★})$, the number of indexes in $j_k^{★}$. If $n_{ck}$ is greater than $c_{max}$, then the $n_{rc}-c_{max}$ smallest clusters will be merged with their nearest neighbouring cluster. After evaluating the loop for all of the superpixel regions in the image, the last step is to update the parameters of each of the pixel clusters $C_j$ in the cluster list $\mathbf{C}$ using (24). The parameter update recalculates the position of the centroids and average pixel intensity of the cluster for use in the distance-image-based cluster optimization in the next iteration of the algorithm.

4.4. Algorithm Iteration and Output

With all of the cluster merging having been performed and the parameters of the pixel clusters updated, the current iteration step of the superpixel algorithm is completed. The next and subsequent algorithm iterations will start again at the cluster update step, followed by the region update step (consisting of the region union and cluster merge operations). The algorithm will continue iterating until a stopping condition occurs, thus terminating the algorithm. There are three cases under which the algorithm will be stopped, with the first case occurring when the statistical dissimilarity between all regions is greater than the threshold ${\lambda }_{KS}$, as no further region union operations can be performed. The second case arises when all of the superpixel regions within the image are larger than the desired superpixel size $r_{max}$, and the third case happens when the current number of superpixel regions (in $\mathbf{R}$) is less than or equal to $r_{ndes}$. The algorithm will be stopped whenever any of the cases occur and the condition can be given as a Boolean value:

$$E_{stop}=\left\{\begin{array}{cc}\mathrm{true}\hfill & \mathrm{if}\underset{\mathcal{E}}{min}\left(D_{KS}(k_{v1},k_{v2})\right)>{\lambda }_{KS}\forall E_v\hfill \\ \mathrm{true}\hfill & \mathrm{if}{n}_{kpx}>r_{max}\forall R_k\hfill \\ \mathrm{true}\hfill & \mathrm{if}\mathrm{Count}(R_k\in \mathbf{R})\le r_{ndes}\hfill \\ \mathrm{false}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(45)

with the first case being evaluated for all of the edges $E_v$ in the graph, and the second and third cases being evaluated for all of the regions in the image. When the algorithm terminates, the resulting output is a label image ${\mathcal{L}}_R(i)$ that indicates the corresponding superpixel region for each pixel in the input image and the superpixel adjacency graph. The region image is created by initializing a new image of the same size as the input image, initializing all pixel values to zero, and then updating the label (value) of all of the pixels based on which superpixel (region) they are contained within, with the labelling given by

$${\mathcal{L}}_R(i)=k,\forall \{i\in C_j,\mathrm{for}\mathrm{each}{C}_j\in R_k\},$$

(46)

such that all of the pixels i in the component clusters $C_j$ of the superpixel $R_k$ are labelled with an index k.

The results of running the superpixel algorithm on an input TRUS image, showing the pixel cluster, superpixel regions, and adjacency graph during multiple algorithm iterations, are given in Figure 6. The boundaries of the pixel clusters are indicated with blue lines and the centroids of the pixel clusters are shown as blue dots in Figure 6a,d,g. The superpixel region boundaries are designated with green lines and the region centroids are shown in Figure 6b,e,h. Red dots are used to indicate the centroids of each of the regions (which are also the graph vertices) in Figure 6b,c,e,f,h,i. The region adjacency graph is drawn in Figure 6c,f,i, where a red line between two region centroids represents a graph edge.

5. Graph-Cut Segmentation

The last step in the prostate segmentation algorithm, as shown in Figure 1, combines the output of the superpixel processing with the statistical prostate shape model. The prostate shape model, described in Section 3.2, can be used to find the likelihood that a particular pixel will be on the edge of the prostate. The superpixel image processing algorithm returns an image which has been divided into a large number of superpixel regions, with each superpixel region being a large area of pixels which have particular pixel intensity distributions. The superpixel algorithm also returns an associated graph describing the adjacency between the superpixel regions. A graph-cut algorithm will be used to merge the contour position information with statistical properties of the superpixel regions in order to segment the prostate.

Given an input set of TRUS ultrasound images slices, ${}^U\mathbf{I}$, containing $nImg$ images, the first step in segmentation is to preprocess each image, as covered in Appendix A, to obtain the ROI for segmentation. After preprocessing, each of the ROI images in the set, ${}^R\mathbf{I}$, will be input into the superpixel processing algorithm. Each image i in the set will then have a corresponding superpixel region image ${}^R{\mathcal{L}}_{Ri}$ and graphs describing the neighbours of each region ${\mathcal{G}}_i=\{{\mathcal{V}}_i,{\mathcal{E}}_i\}$. For final segmentation, the superpixel region image is warped back into the ultrasound image frame $\{U\}$ by inverting the transform given in (A1). To do this, the location of each pixel inside of the ROI in ${}^U\mathcal{I}{({}^{U}px,{}^{U}py)}_i$ will be used to find the corresponding pixel in ${}^R\mathcal{I}{({}^{R}px,{}^{R}py)}_i$, and the region label for that pixel, at ${}^R{\mathcal{L}}_{Ri}({}^{R}px,{}^{R}py)$, will be transferred to a new region image ${}^U{\mathcal{L}}_{Ri}({}^{U}px,{}^{U}py)$, and any pixels outside of the ROI are given a label of 0. For a region $R_k$, the pixels inside of that region are defined to be $\{{}^U{\mathcal{L}}_{Ri}({}^{U}px,{}^{U}py)=k:\forall ({}^{R}px,{}^{R}py)\in R_k\}$; the result of this warping procedure is shown in Figure 9. Warping the region images back in the $\{U\}$ frame is required in order to use the prostate shape model, as it was derived based on the prostate contour segmented by clinicians in standard TRUS images. The probabilistic prostate shape model is generated offline before segmentation and so the prostate can be segmented in all of the input images.

For segmentation, an energy-based cost function [40] will be minimized to perform the final segmentation (contouring) of the prostate in the TRUS images. The segmentation and contouring will be performed by labelling pixel regions $R_k$ as being inside of the prostate, ${\mathbf{L}}_{In}$, or outside of the prostate, ${\mathbf{L}}_{Out}$. For a single image i, the cost function is given by

$$E{(\mathbf{L})}_i={\lambda }_D\sum _{k\in {\mathcal{V}}_i}{\psi }_k({\mathbf{L}}_k)+{\lambda }_S\sum _{k,m\in {\mathcal{E}}_i}{\psi }_{km}({\mathbf{L}}_k,{\mathbf{L}}_m)$$

(47)

where $\mathbf{L}$ is a labelling for the images. The cost function contains a data term

$$\sum _{k\in {\mathcal{V}}_i}{\psi }_k({\mathbf{L}}_k)$$

(48)

for each region $R_k$ in the image corresponding to the vertices in ${\mathcal{V}}_i$, and a smoothness term

$$\sum _{k,m\in {\mathcal{E}}_i}{\psi }_{km}({\mathbf{L}}_k,{\mathbf{L}}_m),$$

(49)

which incorporates a score based on the neighbours connected by an edge ${\mathcal{E}}_i$, of the region $R_k$. The values ${\lambda }_D$ and ${\lambda }_S$ weight the data terms and smoothing terms, respectively, and this cost function will be minimized through the graph-cut procedure.

The data term will incorporate the contour location probabilities from the prostate shape model and a prostate appearance term into the graph-cut score. For this, the probability from the prostate shape model, given in (19), will be discretized to give a pixel-based probability score which will be stored in an image, ${}^U{\mathcal{P}}_i$, of the same size as the TRUS images (with a height ${}^U{\mathcal{I}}_{height}$ and width ${}^U{\mathcal{I}}_{width}$). For mapping between the ultrasound image set and prostate, a centre point, ${}^{P}c={[{}^P{x}_c,{}^P{y}_c,{}^P{\varphi }_c]}^{\top }$, will be used. Only a single contour point is used for the entire image set and, in a similar manner to the contour fitting procedure, this centre point will be optimized during the graph-cut procedure. Using a given centre point, the location of a pixel in the TRUS image, $({}^{U}px,{}^{U}py)$, is converted to the contour coordinate system through the transform

$$\begin{array}{cc}\hfill {}^{P}px& =\frac{{}^{U}px-t_{px}}{s_{px}}\hfill \\ \hfill {}^{P}py& =\frac{{}^{U}py-t_{py}}{s_{py}}\hfill \\ \hfill {\theta }_{px}& =\mathrm{atan}2({}^{P}py-{}^P{y}_c,{}^{P}px-{}^P{x}_c)-{}^P{\varphi }_c\hfill \\ \hfill r_{px}& =\sqrt{{({}^{P}py-{}^P{y}_c)}^2+{({}^{P}px-{}^P{x}_c)}^2}\hfill \end{array}$$

(50)

where ${}^{P}px$ and ${}^{P}py$ are the pixel coordinates in the patient attached frame, using the constants in the transform (A1), giving the output values ${\theta }_{px}$ and $r_{px}$ in the contour coordinate system. Using this transform, the probability image, ${}^U\mathcal{P}$, is created by evaluating the CDF of the prostate shape, given in (20), for every pixel within the ROI of in the input image, such that

$${}^U\mathcal{P}{({}^{U}px,{}^{U}py)}_i=\left\{\begin{array}{cc}P(R\le r_{px}|{\theta }_{px},i)\hfill & \mathrm{for}({}^{U}px,{}^{U}py)\in {}^U{\mathcal{I}}_{ROI}\hfill \\ 1\hfill & \mathrm{for}({}^{U}px,{}^{U}py)\notin {}^U{\mathcal{I}}_{ROI}\hfill \end{array}\right.$$

(51)

where $P(R\le r_{px}|{\theta }_{px},i)$ is the probability that the prostate contour is smaller than the pixel’s radial distance $r_{px}$, relative to the centre point, and the value of 1 will prevent the prostate contour from being outside the ROI.

With the probability image defined, the data term components ${\psi }_k({\mathbf{L}}_k)$ for the cost function (47) are given as

$$\begin{array}{c}{\psi }_k(R_k\in {\mathbf{L}}_{In})={\lambda }_{DP}\left({\displaystyle \sum _{({}^{U}px,{}^{U}pyk)\in k}}\frac{{}^U\mathcal{P}{({}^{U}px,{}^{U}py)}_i}{n_{kpx}}\right)+\\ {\lambda }_{SD}\left(\frac{1}{n_{px}}{\displaystyle \sum _{({}^{U}px,{}^{U}pyk)\in k}}\frac{\left|{}^U\mathcal{L}{({}^{U}px,{}^{U}py)}_i-{\mu }_L\right|}{n_{kpx}}\right)\end{array}$$

(52)

and

$$\begin{array}{c}{\psi }_k(R_k\in {\mathbf{L}}_{Out})={\lambda }_P\left({\displaystyle \sum _{({}^{U}px,{}^{U}pyk)\in k}}\frac{1-{}^U\mathcal{P}{({}^{U}px,{}^{U}py)}_i}{n_{kpx}}\right)+\\ {\lambda }_L\left(\frac{1}{n_{px}}{\displaystyle \sum _{({}^{U}px,{}^{U}pyk)\in k}}\frac{255/{\mu }_L-\left|{}^U\mathcal{L}{({}^{U}px,{}^{U}py)}_i-{\mu }_L\right|}{n_{kpx}}\right)\end{array}$$

(53)

with $n_{kpx}$ being the number of pixels in the region k, and ${\lambda }_{DP}$ weighting contour probabilities. The ${\mu }_L$ value is an input constant corresponding to the average pixel intensity of pixels inside the prostate, and ${\lambda }_{DL}$ weights this pixel intensity term. The two weighting components are represented in (47) as ${\lambda }_D$, such that ${\lambda }_D=[{\lambda }_{DP},{\lambda }_{DL}]$. The smoothness terms within the cost function use the KS distance to weight the statistical distances between regions with the same and different labels, given by

$${\psi }_{km}({\mathbf{L}}_k,{\mathbf{L}}_m)=\left\{\begin{array}{cc}\left(\underset{x}{sup}\left|{\widehat{F}}_k(x)-{\widehat{F}}_m(x)\right|\right)\hfill & \mathrm{where}{\mathbf{L}}_k={\mathbf{L}}_m\hfill \\ \left(1-\underset{x}{sup}\left|{\widehat{F}}_k(x)-{\widehat{F}}_m(x)\right|\right)\hfill & \mathrm{where}{\mathbf{L}}_k\ne {\mathbf{L}}_m\hfill \end{array}\right.,$$

(54)

with ${\widehat{F}}_k(x)$ and ${\widehat{F}}_m(x)$ being the empirical CDFs of the neighbouring regions $R_k$ and $R_m$.

Each TRUS image slice i is then labelled to optimize the cost function using graph cutting [40]. The cost of the particular labelling L is calculated using (47) for all of the slices. To obtain the most accurate contour probability in the images ${}^U{\mathcal{P}}_i$ from the prostate shape model, the centre point ${}^{P}c={[{}^P{x}_c,{}^P{y}_c,{}^P{\varphi }_c]}^{\top }$ is optimized to remove any translational and rotational offset between the TRUS image and the prostate shape model. The centre point is found through the optimization

$$\begin{array}{c}\hfill \underset{{}^{P}c}{argmin}\sum _{i=0}^{nImg}\underset{L}{argmin}E{(\mathbf{L})}_i\end{array}$$

(55)

where the ${argmin}_{L}E{(\mathbf{L})}_i$ corresponds to the graph-cut procedure, and the ${argmin}_{{}^{P}c}$ optimization is carried out through a pattern search algorithm. Once the optimal centre point has been found, with a minimal labelling cost across all of the slices, the prostate segmentation is complete, as shown in Figure 10. All regions labelled with ${\mathbf{L}}_{In}$ after graph cutting are inside of the prostate, which could be represented in a binary image where all of the pixels inside of the regions $R_k\in {\mathbf{L}}_{In}$ have a value of one and all other pixels are zero. To smooth out the segmented prostate contour, the contour model is fit to the pixels on the perimeter of the binary region $R_k\in {\mathbf{L}}_{In}$. The perimeter pixels are transformed to the contour model coordinate system using (50), and then the model is fit to these pixel points from the least-squares fitting given in (6), thus creating a 3D contour model for the segmented prostate; see Figure 11. The results of the segmentation on a number of patient image sets is given in Section 6.

6. Results

The input dataset used to create the probabilistic prostate contour model was a series of contours exported from Variseed (Varian Medical Systems Inc., Palo Alto, CA, USA) after preplanning for clinical prostate brachytherapy (approval for study granted from Alberta Cancer Research Ethics Committee under reference 25837). For this work, the prostate segmentation algorithm was tested on the clinical ultrasound images collected from nine patients during a prostate brachytherapy procedure. The ultrasound images were taken using a B&K Medical Pro Focus 2202 ultrasound machine with a B&K Medical Type 8808 Biplane Transducer (B&K Medical Systems, Peabody, MA, USA). The image processing, model fitting, and graph-cut algorithms were all programmed or implemented in Matlab 2018a (The Mathworks Inc., Natwick, MA, USA) and run using the Simulink Real-Time environment, on an Intel Core i7-3930K running at 3.20 GHz (Intel Corporation, Santa Clara, CA, USA).

The mean absolute difference (MAD), max absolute difference (MaxD), segmented prostate volume, Jaccard index, and prostate contour width–height are used to compare the 3D contour volumes returned by algorithm, $R(\theta ,i|{\mathit{\beta }}_A)$, to the target clinician segmented 3D contours, $R(\theta ,i|{\mathit{\beta }}_C)$, to evaluate the effectiveness and accuracy of the algorithm. The mean absolute difference (MAD) and max absolute difference (MaxD) are commonly used performance metrics in the literature which provide a measure of how closely the component 2D algorithm and clinician segmented contours match across all slices in the 3D contour. The difference between corresponding points, with the same i and $\theta$ values, on the algorithm segmented contour and clinician segmented contours are given by

$$D_{CA}(i,\theta )=\sqrt{{\left({}^P{x}_C(i,\theta )-{}^P{x}_A(i,\theta )\right)}^2+{\left({}^P{y}_C(i,\theta )-{}^P{y}_A(i,\theta )\right)}^2}$$

(56)

where $\{{}^P{x}_C(i,\theta ),{}^P{y}_C(i,\theta )\}$ and $\{{}^P{x}_A(i,\theta ){}^P{y}_A(i,\theta )\}$ are the corresponding points between the algorithm and clinician contours in Cartesian coordinates, respectively, where the Cartesian coordinates are found from (3). A set of ${\theta }_j$ points sampled at regular intervals in the range $(-\pi ,\pi ]$, where $thet{a}_j=-\pi +2\pi \frac{j}{n_{\theta }-1}$, is used for calculating the MAD and MaxD values; the value $n_{\theta }$ is the number of points to be compared ($n_{\theta }=90$ for the presented results). The mean absolute difference and max absolute difference metrics are given by

$$\begin{array}{c}\mathrm{MAD}=\frac{1}{nSlice\times n_{\theta }}\sum _i\sum _j\left|D(i,{\theta }_j)\right|\\ \mathrm{MaxD}=max\left|D(i,{\theta }_j)\right|\forall \{i,j\}\end{array}$$

(57)

where $i=\{0,1,...,nSlice\}$ and $j=\{0,1,...,n_{\theta }\}$. The height and width of the prostate contour, measured in mm in the ultrasound image frame ${}^{U}y$ and ${}^{U}x$ axes, are found from the height and width of a bounding box containing the prostate contour; see Figure 12. The average width and height of the algorithmic and clinician segmented prostate contours, with respect to image depth, are presented in Figure 13. For measuring the area and volume of the prostate, a polygon is fit to points on the prostate contour model $R(\theta ,i|{\mathit{\beta }}_A)$ or clinician-segmented contour $R(\theta ,i|{\mathit{\beta }}_C)$ for each slice i, as shown in Figure 12, with the resulting polygon areas, $Are{a}_A(i)$ and $Are{a}_C(i)$ for the algorithmic and clinician contours respectively, calculated in metric units (mm²). To calculate the volume of the prostate for the algorithm $Vo{l}_A$ and clinician contours $Vo{l}_C$, the area of each slice was multiplied by the 5 mm TRUS image slice stepping distance and summed for all of the slices in the image set, such that

$$\begin{array}{c}Vo{l}_A={\displaystyle \sum _{i=0}^{nSlice}}\frac{5\mathrm{mm}\ast Are{a}_A(i)}{1000}\\ Vo{l}_C={\displaystyle \sum _{i=0}^{nSlice}}\frac{5\mathrm{mm}\ast Are{a}_C(i)}{1000}\end{array}$$

(58)

with the factor of 1000 used to scale the output volumes to mL. The calculated prostate volumes, for each of the nine prostates incorporated in the experimental results, are shown in Figure 14. For the results presented in

Table 2. Prostate segmentation contour results.

Dataset	Mean Absolute Difference (mm)	Max Absolute Difference (mm)	Volume Difference (mL)	Volume Difference
Prostate 1	2.01 (±1.30)	5.89	4.65	8.85
Prostate 2	1.87 (±1.29)	6.91	2.27	4.37
Prostate 3	3.03 (±2.21)	8.50	6.75	13.77
Prostate 4	2.54 (±1.48)	6.00	5.77	11.88
Prostate 5	2.46 (±1.60)	6.82	2.80	5.76
Prostate 6	2.24 (±1.43)	6.70	6.20	12.68
Prostate 7	3.08 (±2.06)	8.82	7.86	17.56
Prostate 8	2.73 (±1.71)	8.90	-3.39	6.07
Prostate 9	2.76 (±1.65)	6.20	7.47	17.11
Average	2.52 (±1.66)	7.19 (±1.22)	4.49 (±3.53)	10.89 (±4.90)

, the volume difference in mL and % were calculated by

$$\begin{array}{l}\mathrm{Volume}\mathrm{Difference}(\mathrm{mL})=Vo{l}_C-Vo{l}_A\\ \mathrm{Volume}\mathrm{Difference}(\%)=\frac{\left|Vo{l}_A-Vo{l}_C\right|}{Vo{l}_C}\end{array}$$

(59)

where a positive value for the prostate contour volume difference indicates that the algorithm segmented prostate contour was larger than the clinician-segmented contour, and a negative difference indicates that the algorithmic segmented contour was smaller.

As an additional metric to evaluate the performance of the prostate contouring algorithm, the Jaccard index was also analyzed. The Jaccard index gives a measure between 0 and 1 that indicates how well the algorithmic output of the prostate contour matches with the clinician-segmented contour, with a value of 1 indicating a perfect match. The Jaccard index is given by

$$JI=\frac{\left|R(\theta ,i|{\mathit{\beta }}_A)\cap R(\theta ,i|{\mathit{\beta }}_C)\right|}{\left|R(\theta ,i|{\mathit{\beta }}_A)\cup R(\theta ,i|{\mathit{\beta }}_C)\right|}$$

(60)

where $R(\theta ,i|{\mathit{\beta }}_A)$ is the output contour of our proposed prostate segmentation algorithm, for a particular TRUS image set, and $R(\theta ,i|{\mathit{\beta }}_C)$ is the clinician-segmented contour for the same image set. Figure 15 shows the average Jaccard index score values with respect to the z-axis depth of the TRUS image, where the first image slice at a depth of 0 mm corresponds to the base of the prostate, the midgland of the prostate is located at 15 mm, and the apex of the prostate is at a depth of 30 mm, as indicated in Figure 11.

Table 2 and

Table 3. Prostate segmentation Jaccard index results.

Dataset	Average Jaccard Index	Min Jaccard Index	Max Jaccard Index
Prostate 1	0.84 (±0.04)	0.79	0.92
Prostate 2	0.85 (±0.06)	0.75	0.91
Prostate 3	0.77 (±0.02)	0.73	0.79
Prostate 4	0.79 (±0.06)	0.73	0.87
Prostate 5	0.79 (±0.04)	0.72	0.83
Prostate 6	0.82 (±0.06)	0.71	0.90
Prostate 7	0.74 (±0.07)	0.66	0.85
Prostate 8	0.77 (±0.07)	0.66	0.84
Prostate 9	0.77 (±0.09)	0.61	0.89
Average	0.79 (±0.07)	0.71 (±0.05)	0.87 (±0.04)

show the resulting values of the metrics, comparing the output contours from our proposed prostate segmentation algorithm (the input parameters used for the superpixel algorithm for the results given in Table 2 were ${\lambda }_{\mu }=100,{\lambda }_{L2}=125,{\lambda }_{KS}=0.2,$ $s_{grid}=256\mathrm{px},r_{darea}=512\mathrm{px},r_{min}=256\mathrm{px},c_{min}=192\mathrm{px},c_{max}=4$, and the parameters for the graph-cut segmentation were ${\lambda }_{DP}=3.1,{\lambda }_{DL}=1.8,{\mu }_l=40,{\lambda }_S=1.5$) for nine patient TRUS image sets, where each image set consisted of seven images slices. Table 2 presents the mean absolute difference, max absolute difference, and prostate volume difference, and Table 3 presents the contour difference, Jaccard index, and prostate volume.

The average mean absolute difference (average MAD) for all of the trials was 2.52 (±1.66) mm and the average maximum absolute difference (average MaxD) was 7.19 (±1.22) mm. For context, using values reported by [9], we will discuss these results with respect to three semi-autonomous 3D prostate TRUS segmentation algorithms within the literature [14,26,27] that are the most comparable to our proposed method. The active contour model algorithm proposed by [14] reported an MAD of $7.01\pm 1.04$ mm and an MaxD of $1.09\pm 0.49$ mm. The active-shape-model-based contouring algorithm of [26] provided an MAD of $7.27\pm 2.32$ mm and an MaxD of $1.26\pm 0.41$ mm. The hybrid active contour/active shape model system developed in [27] had an MaxD of $4.06\pm 1.25$. While the mean absolute difference of our algorithm is slightly larger than comparable results presented in the literature, the maximum absolute difference is in line with those reported in the literature. A major factor which may account for this discrepancy in mean absolute difference performance between our approach and those in the literature is that many methods, including [14,26], require a clinician to manually correct the algorithmic prostate contours after segmentation [9,10,11], and report only the post-corrected segmentation values. In many cases in the literature, this manual clinician correction is required in up to 30% of the input set TRUS images for segmentation.

As a deep comparison, in the literature it is well known that there is a high amount of variation between clinicians when manually segmenting the prostate in TRUS images [41,42]. Therefore, the performance of the algorithm can also be compared to clinician segmentation of TRUS images, where [41,42] analyzed the interobserver clinician segmentation accuracy between observers (clinicians) and intraobserver repeatability from a single observer (a clinician segmenting the same TRUS image set multiple times). Interestingly, the difference between the prostate volume measured from the algorithm and clinician-segmented contours is directly comparable to the interobserver variance reported in [41,42], indicating that the segmented contours from the proposed algorithm are within the accuracy range accepted by clinicians in practice.

The average Jaccard index values across all images slices are slightly lower than those reported for clinicians [42], indicating that some additional work to increase contouring accuracy may be required. The per-slice Jaccard index values presented in Figure 15 indicate that the contouring has roughly the same contouring performance across all prostate image slices, with the exceptions of the base image slice (0 mm), which has slightly lower than average accuracy, and the apex image slice (30 mm), which has higher than average accuracy. These per-slice results compare favourably with the variances reported in [42], with the base slice segmentation variance being higher for both the algorithm and between clinicians. The accuracy of the prostate midgland contouring is the primary area where the algorithm performs slightly worse than clinicians [42] and indicates the primary area where the algorithm accuracy can be improved.

7. Conclusions

In this work, we have outlined a method for prostate segmentation which incorporates a prostate contour model, active shape probabilistic model of the prostate, a superpixel image processing algorithm, and graph-cut-based contour optimization. The proposed contour model was validated independently of the prostate segmentation algorithm by comparing the mean-squared error and maximum absolute error between input, clinician-segmented, prostate contours and the resulting contour models after parameter fitting. Two sets of 100 3D clinician-segmented prostate contours were used for validation, with one set consisting of 3D contours containing 7 to 13 2D TRUS image-derived contours of the prostate and the second set having 3D contours with each containing 7 to 16 2D contours segmented using MR images of the prostate. The results of this validation, given in Table 1, show that the clinician-segmented contours, from either imaging modality, can be represented in the contour model with sub-millimetre accuracy on average. The average mean-squared error across the two input contour sets was 0.61 (±0.16) mm and the average maximum absolute error was 1.54 (±0.31) mm. The largest absolute error observed in parameter fitting of all 200 contours was 2.12 mm. With these results, the prostate contour model was shown to be unaffected by imaging modality and capable of modelling a variety of prostate shapes with an acceptable accuracy.

The prostate contour model was then incorporated, using interpolation between Gaussian normal distribution models, into an active shape statistical contour model of the prostate. This statistical shape model describes the probability of the edge of the prostate being at a particular location in a TRUS image slice. The CDF for the statistical shape model was incorporated as a parameter within the graph-cutting cost function to utilize prior prostate contour information for algorithmic prostate segmentation. The graph to be segmented was output by a novel superpixel image processing algorithm which uses a hierarchical structure to simultaneously optimize pixel clusters, with respect to pixel intensity and Euclidean-based distance metrics, and superpixel regions, which are composed of a number of pixel clusters. These superpixel regions separate an input image into areas with similar image statistics, in a type of edge and texture discrimination, with the statistical similarity between two superpixel regions being another parameter within the graph-cutting cost function.

The superpixel image processing algorithm created and updated a graph structure, with the superpixel regions corresponding to the vertices of the graph and neighbouring superpixel regions being interconnected by a graph edge. Graph-cuts are performed on the resulting superpixel adjacency graph, after image processing, in order to segment out superpixel regions which are inside of the prostate in the TRUS image. This prostate segmentation algorithm was tested on nine prostate TRUS image sets, consisting of six or seven prostate image slices, to assess the algorithmic segmentation accuracy. The output of the prostate segmentation algorithm was compared to prostate segmentation performed by a clinician, in those TRUS image sets, with an average mean absolute difference of 2.52 (±1.66) mm and an average maximum absolute difference of 7.19 (±1.22) mm. The average Jaccard index of the algorithm was 0.79 (±0.07). The average Jaccard index result is satisfactory for the limited number of TRUS image sets that were tested; however, the continuation of this work will require a significantly larger number of the TRUS prostate image sets to fine-tune the optimization parameters used during superpixel image processing and segmentation.