A semiautomatic tool for prostate segmentation in radiotherapy treatment planning

Abstract

Background

Delineation of the target volume is a time-consuming task in radiotherapy treatment planning, yet essential for a successful treatment of cancers such as prostate cancer. To facilitate the delineation procedure, the paper proposes an intuitive approach for 3D modeling of the prostate by slice-wise best fitting ellipses.

Methods

The proposed estimate is initialized by the definition of a few control points in a new patient. The method is not restricted to particular image modalities but assumes a smooth shape with elliptic cross sections of the object. A training data set of 23 patients was used to calculate a prior shape model. The mean shape model was evaluated based on the manual contour of 10 test patients. The patient records of training and test data are based on axial T1-weighted 3D fast-field echo (FFE) sequences. The manual contours were considered as the reference model. Volume overlap (Vo), accuracy (Ac) (both ratio, range 0-1, optimal value 1) and Hausdorff distance (HD) (mm, optimal value 0) were calculated as evaluation parameters.

Results

The median and median absolute deviation (MAD) between manual delineation and deformed mean best fitting ellipses (MBFE) was Vo (0.9 ± 0.02), Ac (0.81 ± 0.03) and HD (4.05 ± 1.3)mm and between manual delineation and best fitting ellipses (BFE) was Vo (0.96 ± 0.01), Ac (0.92 ± 0.01) and HD (1.6 ± 0.27)mm. Additional results show a moderate improvement of the MBFE results after Monte Carlo Markov Chain (MCMC) method.

Conclusions

The results emphasize the potential of the proposed method of modeling the prostate by best fitting ellipses. It shows the robustness and reproducibility of the model. A small sample test on 8 patients suggest possible time saving using the model.

Background

Prostate cancer is the second most diagnosed cancer accounting for 14 percent of all cancers diagnosed worldwide [1]. It is most common in males over the age of 50, and has the highest incidence rate in the developed countries. Aggressive tumors are usually treated with extern radiotherapy or brachytherapy which requires a precise treatment plan for the target volume. In any type of radiotherapy treatment, radiation of healthy tissue should be minimized while maintaining the desired dose to the target volume. Therefore, a successful treatment of prostate cancer relies on an accurate segmentation of the prostate from the surrounding tissue, by image-based description of the shape and location of the target volume. The volume of interest is characterized by a smooth shape, and for this reason an algorithmic description of the volume is feasible.

Transrectal ultrasound (TRUS), magnetic resonance (MR) and computed tomography (CT) images are the three main imaging techniques used in diagnosis, treatment planning and follow-up examination of prostate cancer. Smith et al. [2] investigated the effects of these imaging techniques on the properties of the prostate volume. A collection of methods available for prostate segmentation is reviewed by Ghose et al. [3]. In addition to the methods presented by Ghose et al., alternative approaches are available in the literature, such as the medial or skeleton representation of the prostate [4-8]. The present work proposes a segmentation method which falls into the category of deformable meshes in Ghose et al. [3], but refers to the term geometrical parametrization as described in Dryden and Mardia [9]. The main focus of this paper is the development of a statistical shape model for the prostate. An overview about this type of models in 3D medical image segmentation is presented for example by Davies et al. [10] and Heimann and Meinzer [11].

The works of Saroul et al. [12] and Mahdavi et al. [13] are related to the stacked ellipses parametrization method used in this paper. Mahdavi et al. [13] proposes a 3D ellipsoid shape of the prostate in warped transrectal 3D ultrasound images based on control points. This method extends the warping idea proposed in Badiei et al. [14] from 2D to 3D ultrasound images. On the contrary, we focus on slice-wise best fitting ellipses which will introduce more flexibility into the model, e.g., between the positions and lengths of the first and second axes of the ellipses between neighbor slices. The approach of slice-wise best fitting ellipses has similarities to a tubular medial representation [15].

Beside the single segmentation of the prostate, several attempts have been tried out for a joint segmentation of neighbor organ and structure to gain improved segmentation results [16-18].

To our knowledge, despite the substantial effort in this area, no widely implemented algorithm exists. In oncology departments this means that the physician has to delineate the prostate slice by slice. This is time-consuming and inefficient. We propose a less ambitious approach compared to more sophisticated models, such as skeletal models as discussed above, in that we use a method that gives a useful starting point for the physician after the definition of few control points. Given the initial estimate of the volume of interest, the physician can adjust the estimate according to their evaluation of the image rather than starting from scratch. By this approach, we obtain the same accuracy with less effort. The main points in our approach are as follows: First, we accept that the algorithm cannot give a fully precise description of the volume. Our main aim is therefore to give a good estimate which can be used as a starting point for the physician. Second, we use a simple ellipse model that is easy to interpret and understand. Our hypothesis is that a more efficient use of physicians in Radiotherapy Treatment Planning (RTP) of patients with prostate cancer can be obtained by an easy-to-interpret semiautomatic tool.Figure 1 shows an example of the initial estimate we typically obtain for a single image slice. The dashed line in (a) to (e) describes the manual contour while the solid line shows the best fitting ellipse including the two principal axes for the observed data of this slice. Note that the fitted model is very much in agreement with the manual line, indicating that the stacked ellipses model gives a good description of the object of interest. The solid lines in (f)-(j) shows the outcome from our model in this situation together with few defined control points.

Figure 1

Selected slices of MR data set 3 from the test data set. (a)-(e) Manual delineation of the prostate (dashed line) and best fitting ellipse (solid line). (f)-(j) Manual delineation of the prostate (dashed line), deformed mean shape (solid line) and defined control points in the first, center and last slice.

This result shows a typical performance of the method, and that the estimate is close to the best fit we can obtain with the ellipse model. The full processing demands little computational resources, such that the suggested delineation can be presented immediately. The example is discussed further in the Methods and Results and discussion section.

The rest of the paper is organized as follows. In the Methods section, we introduce the data sources and the proposed stacked ellipses model, and discuss the shape space and statistics along with constraints and parameters. Results are presented in the Results and discussion section using a test data set to show the potential of the mean shape model, followed by a Conclusions section. Additional file 1 with further detailed discussion is available online.

Methods

Preliminaries

Each prostate must be described by a shape model in order to calculate statistics, e.g., by stacked ellipses as a parametric shape model. The parameters of a parametric shape model can be estimated from a training set. The training set models also the geometric variability of anatomical structures by a shape probability distribution. The training set contains volume and contour information of segmented prostates from N patients. The volume information describes the image modalities (e.g., CT or MR) and the contour information the volume of interest as defined in the following.

The volume information of each training set n = 1, …,N is defined by a 3-dimensional matrix V _n where V _n (i,h) contains the observed gray level in voxel (i,h), i = (i ₁,i ₂) ∈ {1,…,I ₁} × {1,…,I ₂} are the pixel indices in a slice, where typically I ₁ = I ₂, and h ∈ {1,…,H} is the number of slices per data set. The number of slices H is not necessarily the same for all patients in the training data sets. Therefore, we indicate H by H _n and in the same manner I ₁ by I _n1 and I ₂ by I _n2, but for simplicity we use H, I ₁ and I ₂ if the meaning is clear.

In addition to the volume information, each training set n =1,…,N consists of contour information of the prostate, manually drawn by a physician. The contour information can be modeled by a (M × K _n ) configuration matrix $X_n:=(X_{n1},\dots ,X_{{\mathit{\text{nK}}}_n})$ with $X_{\mathit{\text{nk}}}={(x_{1k}^n,x_{2k}^n,x_{3k}^n)}^T\in {\mathbb{R}}^3,k=1,\dots ,K_n$, where K _n defines the total number of available contour information points in a data set and M = 3 defines the dimension. We assume the contour information for an object is defined in a sequentially sorted number L _n of equidistant slices whereas each contour slice contains ${\tilde{K}}_{\mathit{\text{nl}}}$ contour points, l = 1,…,L _n . Hence it follows $K_n=\sum _l{\tilde{K}}_{\mathit{\text{nl}}}$ and $X_n=({\tilde{X}}_{n1},\dots ,{\tilde{X}}_{{\mathit{\text{nL}}}_n})$. The image information in slice l is denoted by S _nl and $S_n=\{S_{n1},\dots ,S_{{\mathit{\text{nL}}}_n}\}\subseteq V_n$ and ${\tilde{X}}_{\mathit{\text{nl}}}$ defines the configuration matrix in slice S _nl .

In summary, the training population is given by the set $\{\mathbb{V},\mathbb{X}\}$, with a set of volume information $\mathbb{V}=\{V_1,\dots ,V_N\}$ and configuration matrices $\mathbb{X}=\{X_1,\dots ,X_N\}$. We assume X _n defines the configuration matrix for the corresponding data set V _n and matches the volume information V _n exactly.

The contour information is often defined in a Patient based Coordinate System (PCS) whereas the volume information is defined in an Image based Coordinate System (ICS). The ICS can be transformed to PCS by a transformation matrix Λ _DCM , which transform an image coordinate p ^im= (i ₁,i ₂,h)^T to patient coordinate p ^p= (x,y,z)^T. The definition of Λ _DCM and the relation between PCS and ICS (see Figure 2) is discussed in detail in the Additional file 1. In addition, we introduce a de-rotated PCS where volume and contour information are aligned to each other.

Figure 2

Visualization of different coordinate systems with example data. PCS: Patient based coordinate system (manual delineation line). ICS: Image coordinate system (volume data). de-rotated PCS: de-rotated patient based coordinate system with same scale and origin as PCS but same orientation as ICS (de-rotated best fitting ellipse dBFE _nl , de-rotated control points dCP _n ). sample space: The transformation matrix Λ dCP n maps the de-rotated data {dBFE _n ,dCP _n } to $\{{\overline{\mathit{\text{BFE}}}}_n,{\overline{\mathit{\text{CP}}}}_n\}$ in the sample space.

Modeling

The prior information inferred from the training set is incorporated into a shape model. We assume a stacked ellipse model as a shape prior for the prostate. Specifically, the prostate outline in slice S _nl ,l = 1,…,L _n , n = 1,…,N is modeled by a slicewise best-fitting ellipse, as visualized in Figure 3. An ellipse in slice S _nl can be uniquely described by ${\rho }_{\mathit{\text{nl}}}={({\theta }^{\mathit{\text{nl}}},{\alpha }^{\mathit{\text{nl}}},{\varphi }^{\mathit{\text{nl}}})}^T\in {\mathbb{R}}^2\times {\mathbb{R}}_{+}^2\times (-\frac{\pi }{2},\frac{\pi }{2}]$with

• position ${\theta }^{\mathit{\text{nl}}}={({\theta }_1^{\mathit{\text{nl}}},{\theta }_2^{\mathit{\text{nl}}})}^T\in {\mathbb{R}}^2$ defines the center in slice S _nl ,

Figure 3

The prostate model of stacked slicewise best-fitting ellipses illustrated by (a) a 3D view of the model and (b) the corresponding 2D model of the prostate contour in a slice.

• length of principle axes ${\alpha }^{\mathit{\text{nl}}}={({\alpha }_1^{\mathit{\text{nl}}},{\alpha }_2^{\mathit{\text{nl}}})}^T\in {\mathbb{R}}_{+}^2$ and

• rotation angle ${\varphi }^{\mathit{\text{nl}}}\in (-\frac{\pi }{2},\frac{\pi }{2}]$ in slice S _nl .

The rotation parameter ϕ ^nlis defined corresponding to the ICS with origin θ ^nlin slice S _nl . The boundary of an ellipse ρ _nl centered at ${\theta }^{\mathit{\text{nl}}}\in {\mathbb{R}}^2$ in slice S _nl is defined by

$$C\left({\rho }_{\mathit{\text{nl}}}\right)=\left\{\mathit{\text{Rx}}+{\theta }^{\mathit{\text{nl}}}:\frac{x_1^2}{{\left({\alpha }_1^{\mathit{\text{nl}}}\right)}^2}+\frac{x_2^2}{{\left({\alpha }_2^{\mathit{\text{nl}}}\right)}^2}=1,x\in {\mathbb{R}}^2\right\}$$

(1)

$$\text{and}R=\left(\begin{array}{cc}\text{cos}{\varphi }^{\mathit{\text{nl}}}& -\text{sin}{\varphi }^{\mathit{\text{nl}}}\\ \text{sin}{\varphi }^{\mathit{\text{nl}}}& \text{cos}{\varphi }^{\mathit{\text{nl}}}\end{array}\right)$$

is a rotation matrix in ${\mathbb{R}}^2$ with rotation angle ϕ ^nland x = (x ₁,x ₂)^T.

The shape model described in this section requires the best fit of an ellipse C (ρ _nl ) to the contour information ${\tilde{X}}_{\mathit{\text{nl}}}$ in each slice, i.e., we model ${\tilde{X}}_{\mathit{\text{nl}}}=C\left({\rho }_{\mathit{\text{nl}}}\right)+\epsilon$ where ε is an error with mean zero. The best-fitting ellipses provide us with a slice-by-slice parametrization of the prostate for all slices in each training shape.

The problem of fitting an ellipse to geometric features like the contour is discussed widely in the literature (e.g., [19,20]). This work follows Ahn et al. [19], who proposed a least-square minimizer for ${\tilde{X}}_{\mathit{\text{nl}}}$. The nonlinear estimate of parameters ${\rho }_{\mathit{\text{nl}}}={({\theta }_1^{\mathit{\text{nl}}},{\theta }_2^{\mathit{\text{nl}}},{\alpha }_1^{\mathit{\text{nl}}},{\alpha }_2^{\mathit{\text{nl}}},{\varphi }^{\mathit{\text{nl}}})}^T$ given ${\tilde{X}}_{\mathit{\text{nl}}}$ must minimize the error

$$g\left({\widehat{\rho }}_{\mathit{\text{nl}}}\right)={\left({\tilde{X}}_{\mathit{\text{nl}}}-\tilde{C}\left({\widehat{\rho }}_{\mathit{\text{nl}}}\right)\right)}^T\left({\tilde{X}}_{\mathit{\text{nl}}}-\tilde{C}\left({\widehat{\rho }}_{\mathit{\text{nl}}}\right)\right)$$

where $\tilde{C}\left({\widehat{\rho }}_{\mathit{\text{nl}}}\right)$ is a set of nearest orthogonal points of ${\tilde{X}}_{\mathit{\text{nl}}}$ to $C\left({\widehat{\rho }}_{\mathit{\text{nl}}}\right)$.

Definition 1

(Best fitting ellipse (BFE)). A best fitting ellipse for slice S _nl is defined by the set ${\mathit{\text{BFE}}}_{\mathit{\text{nl}}}:={({\theta }^{\mathit{\text{nl}}},{\alpha }^{\mathit{\text{nl}}},{\varphi }^{\mathit{\text{nl}}})}^T\in {\mathbb{R}}^2\times {\mathbb{R}}_{+}^2\times (-\frac{\pi }{2},\frac{\pi }{2}]$, l = 1,…,L _n , n = 1,…,N and minimizes the error function g, i.e., ${\mathit{\text{BFE}}}_{\mathit{\text{nl}}}={\widehat{\rho }}_{\mathit{\text{nl}}}$ with

$$g\left({\widehat{\rho }}_{\mathit{\text{nl}}}\right)=\underset{{\rho }_{\mathit{\text{nl}}}\in {\mathbb{R}}^2\times {\mathbb{R}}_{+}^2\times (-\frac{\pi }{2},\frac{\pi }{2}]}{\text{min}}g\left({\rho }_{\mathit{\text{nl}}}\right).$$

(2)

The first and second principal axes must be reordered after calculation of ${\mathit{\text{BFE}}}_n=\{{\mathit{\text{BFE}}}_{n1},\dots ,{\mathit{\text{BFE}}}_{{\mathit{\text{nL}}}_n}\}$ in order to establish correspondence between parameters of adjacent slices and across the population. Improved correspondence will support accurate statistics. The basic idea in our reordering procedure is to carry out the reordering corresponding to the lowest rotation angle of both principal axes to the first principal axis of the neighbor slice where the center slice is chosen as the basis. The rotation between the center slice M and an arbitrary slice is constrained by max(|ϕ ⁱ- ϕ ^M|) = π, i∈ {1,…,L} after reordering. Therefore, the set BFE _n of reordered best-fitting ellipses is an element of $({\mathbb{R}}^2\times {\mathbb{R}}_{+}^2\times {(-\pi ,\pi ])}^{L_n}$.

A further improvement of correspondence is achieved by the introduction of two additional constraints in the parameter model.

First, we relax the rotation parameter ϕ ^nlin case of circularity. If both principal axes have the same length, the orientation of an ellipse is undefined. Therefore we penalize ϕ ^nlin the case of high circularity by taking ${\varphi }^{n{l}^{\prime }}$ from the neighboring slices into account. Second, smoothing is performed between neighboring slices to avoid large forward and backwards rotations between ϕ ^{n (l - 1)},ϕ ^nland ϕ ^{n (l + 1)}. The reordering algorithm and implementation of constraints are described in detail in Additional file 1.

The current implementation assumes the definition of control points CP _n in the training data set {V _n ,X _n ,BFE _n }, where ${\mathit{\text{BFE}}}_n\in ({\mathbb{R}}^2\times {\mathbb{R}}_{+}^2\times {(-\pi ,\pi ])}^{L_n}$ is a reordered set of best fitting ellipses, n = 1,...,N. Furthermore, the control points have to be defined manually by a physician in a new patient data set. The control points are used to make the best fitting ellipses BFE _n comparable and to transform the parametrized ellipses model to a common position, scale and orientation by a transformation matrix Λ dCPn. The transformation matrix Λ dCPn maps the de-rotated prior data {dBFE _n ,dCP _n } to $\{{\overline{\mathit{\text{BFE}}}}_n,{\overline{\mathit{\text{CP}}}}_n\}$ in the sample space, as depicted in Figure 2. In this article, we assume 6 control points in the first, center and last slice at the boundary of the prostate, i.e.,

$${\mathit{\text{CP}}}_n=\left\{A_1^n,\dots ,A_6^n,P_1^n,\dots ,P_6^n,B_1^n,\dots ,B_6^n\right\}$$

as visualized in Figure 3. In addition, we have tested alternative control point configurations. They are described together with the construction of Λ dCPn in Additional file 1.

After transformation we have obtained a reordered and comparable set of best fitting ellipses

$${\overline{\mathit{\text{BFE}}}}_n=\left\{{\overline{\mathit{\text{BFE}}}}_{n1},\dots ,{\overline{\mathit{\text{BFE}}}}_{{\mathit{\text{nL}}}_n}\right\}$$

with ${\overline{\mathit{\text{BFE}}}}_{\mathit{\text{nl}}}={({\overline{\theta }}^{\mathit{\text{nl}}},{\overline{\alpha }}^{\mathit{\text{nl}}},{\overline{\varphi }}^{\mathit{\text{nl}}})}^T$, n = 1,…,N, l = 1,…,L _n . The statistical analysis of the training data requires an equal number L ₁ = … = L _N to establish correspondence between the parameters of the best fitting ellipses. Therefore, we interpolate the set ${\overline{\mathit{\text{BFE}}}}_{\mathit{\text{nl}}}$ to a common number L.

When L is chosen, interpolation is done by independent cubic interpolation in each dimension, i.e., we find points of a one-dimensional function that underlies the data ${\overline{\theta }}_1^{\mathit{\text{nl}}},{\overline{\theta }}_2^{\mathit{\text{nl}}},{\overline{\theta }}_3^{\mathit{\text{nl}}},{\overline{\alpha }}_1^{\mathit{\text{nl}}},{\overline{\alpha }}_1^{\mathit{\text{nl}}}$ and ${\overline{\varphi }}^{\mathit{\text{nl}}}$. The final interpolated best fitting ellipses are denoted by

$${\mathit{\text{iBFE}}}_n=\left\{{\mathit{\text{iBFE}}}_{n1},\dots ,{\mathit{\text{iBFE}}}_{\mathit{\text{nL}}}\right\}.$$

(3)

These ellipses are used for the statistical analysis and computation of a mean shape model. To keep things simple, we denote such a reordered, transformed and interpolated set of best-fitting ellipses by BFE _nl = (θ ^nl,α ^nl,ϕ ^nl)^Tfor the number L of contour slices with l = 1,…,L and n = 1,…,N. The comparable set of best fitting ellipses BFE _n  is an element of the shape space $({\mathbb{R}}^2\times {\mathbb{R}}_{+}^2\times {(-\pi ,\pi ])}^L$.

Statistical analyses

After reconstruction of our shape space we estimate the expectation and variance of the parameters of a mean shape model μ _BFE = {μ BFE 1,…,μ BFE L} with ${\mu }_{\mathit{\text{BFE}}}^l=({\mu }_{\theta }^l,{\mu }_{\alpha }^l,{\mu }_{\varphi }^l)$ from the training set BFE _nl , l = 1,…,L. We denote the mean shape mean best fitting ellipses (MBFE). In addition to the described ellipse parameters we define the position ${\theta }^{\mathit{\text{nl}}}={({\theta }_1^{\mathit{\text{nl}}},{\theta }_2^{\mathit{\text{nl}}},{\theta }_3^{\mathit{\text{nl}}})}^T$ in terms of a distance vector η ^nlof θ ^nlto a center curve defined by the control points. We model θ ^l= ξ ^l+ η ^l, where ξ ^lis analytically defined by L intersection points of the curve within each slice. Thereby, we are describing the mean shape which is closest to the control points. This approach is reasonable under the assumption that the control points are well defined. In Additional file 1 we explore various ways of describing the position parameter for different control point methods.

The mean curve of the expected location is given by

$${\mu }_{{\theta }_j}^l=\frac{1}{N}\sum _{i=1}^N{\theta }_j^{\mathit{\text{il}}},j\in \left\{1,2,3\right\},$$

(4)

where ${\mu }_{\theta }^l={({\mu }_{{\theta }_1}^l,{\mu }_{{\theta }_2}^l,{\mu }_{{\theta }_3}^l)}^T$, l = 1,…,L. The variance and covariance are estimated by

$${\left({\sigma }_{{\theta }_j}^l\right)}^2=\frac{1}{N-1}\sum _{i=1}^N{({\theta }_j^{\mathit{\text{il}}}-{\mu }_{{\theta }_j}^l)}^2,j\in \{1,2,3\},\text{and}$$

(5)

$${\Sigma }_{\theta }^l=\frac{1}{N-1}\sum _{i=1}^N({\theta }^{\mathit{\text{il}}}-{\mu }_{\theta }^l){({\theta }^{\mathit{\text{il}}}-{\mu }_{\theta }^l)}^T.$$

(6)

The length parameter is modeled by a log-normal distribution because $\alpha \in {\mathbb{R}}_{+}^2$. Thus we estimate the mean and variance of $a=\text{log}(\alpha )\in {\mathbb{R}}^2$. The estimation of means and variances of the remaining parameters a,ϕ,η is according to (4-5).

Following Dryden and Mardia [9] we suggest a prior distribution for a new data set as

$$\begin{array}{l}{\theta }_1^l\sim N\left({\mu }_{{\theta }_1}^l,{\left({\sigma }_{{\theta }_1}^l\right)}^2\right),{\theta }_2^l\sim N\left({\mu }_{{\theta }_2}^l,{\left({\sigma }_{{\theta }_2}^l\right)}^2\right),\\ a_1^l\sim N\left({\mu }_{a_1}^l,{\left({\sigma }_{a_1}^l\right)}^2\right)\iff \\ {\alpha }_1^l\sim \text{log-}N\left({\mu }_{a_1}^l,{\left({\sigma }_{a_1}^l\right)}^2\right)\text{with}{a}_1^l=\text{log}\left(\underset{1}{\overset{l}{\alpha }}\right),\\ a_2^l\sim N\left({\mu }_{a_2}^l,{\left({\sigma }_{a_2}^l\right)}^2\right)\iff \\ {\alpha }_2^l\sim \text{log-}N\left({\mu }_{a_2}^l,{\left({\sigma }_{a_2}^l\right)}^2\right)\text{with}{a}_2^l=\text{log}\left(\underset{2}{\overset{l}{\alpha }}\right),\\ {\varphi }^l\sim N\left({\mu }_{\varphi }^l,{\left({\sigma }_{\varphi }^l\right)}^2\right),\end{array}$$

l = 1,…,L. If θ ^lis defined according to the center curve given by the control points as described above, we model ${\eta }_i^l\sim N\left({\mu }_{{\eta }_i}^l,{\left({\sigma }_{{\eta }_i}^l\right)}^2\right),i=1,2$. Since the rotational parameter is expected to have small variance it is not necessary to apply a circular distribution, and we assume normality.

After constructing the shape model we estimate the best fitting ellipse BFE _l parametrized by ρ _l = (θ ^l,α ^l,ϕ ^l)^T, l = l,…,L in a new data set given the control points CP. This is obtained through the posterior π (ρ ∣ S) where s _il ∈ S ⊆ V is the volume information and i = (i ₁,i ₂) ∈ I (ρ) is a set of indices within the ellipses ρ. The control points CP  are used to deform the prior model π (ρ). Therefore we model the posterior by an empirical Bayes approach [21]. The posterior

$$\pi (\rho \mid S,\mathit{\text{CP}})\propto L(S\mid \rho )\ast \pi (\rho \mid \mathit{\text{CP}})$$

(7)

defines the posterior density of the deformed template π (ρ ∣ C P) given the the observed image. The Likelihood or image model L (S  ∣ ρ) is the joint probability density function of the gray levels given the parametrized object ρ | _CP , while ρ | _CP defines the ellipses ρ deformed by the control points CP. The prior π (ρ ∣ C P) models realistic variations from our mean shape ${\mu }_{\mathit{\text{BFE}}}\in {({\mathbb{R}}^2\times {\mathbb{R}}_{+}^2\times (-\pi ,\pi \left]\right)}^L$ given the control points. We are estimating the posterior distribution using a Markov chain Monte Carlo (MCMC) approach. The method and results are discussed in detail in Additional file 1.

Evaluation

We have evaluated the proposed method using 33 patient case studies. The training data set consists of N = 23 T1-weighted Fast Field Echo (FFE) 3D Magnetic Resonance (MR) data. The mean shape model and variance is calculated from the training data set and applied to a test data set of 10 MR FFE case studies. The splitting in test and training data is done according to the sequence of data acquisition. Each data set consists of H _n Digital Imaging and Communications in Medicine (DICOM) image files and one DICOM region structure file, while the contour information of the prostate is stored in the header of a DICOM file without any image information. The voxel size (l _x ,l _y ,l _z ) is (0.559mm,0.559mm,3mm) with a slice distance of 3.3mm of the data sets V _n . Each slice consists of 288 × 288 voxels. The average number of slices with manual prostate contour information is 10.478±2.626 (mean±standard deviation) in the training set and 10.5 ± 2.799 in the test set where the mean is given by ${\mu }_L=\frac{1}{N}\sum _{n=1}^N{L}_n$ and the standard deviation by ${\left(\frac{1}{N-1}\sum _{n=1}^N{(L_n-{\mu }_L)}^2\right)}^{\frac{1}{2}}$. Figure 1 illustrates test patient 3, whose image set consists of 24 MR FFE slices whereas 12 slices contain contour information. Slice 6 is the first slice where contour information of the prostate is available and the last slice is 17.

An ethics approval was not required for this study under Norwegian law because the aim of this study was to develop a tool and not to obtain new knowledge on medicine or diseases. The study uses solely data that were collected during routine medical treatment independent of this study at the University hospital Northern Norway (UNN). Data was provided after a full anonymization following the required guidelines at UNN.

Three metrics are used to compare the manual and the semi-automatic contours. In the axial slices, where the expert manual delineations are present, we calculate the Hausdorff distance (HD) by

$$d_{\mathit{\text{hd}}}(X,Y)=\text{max}\left\{\underset{x\in X}{\text{max}}\underset{y\in Y}{\text{min}}d(x,y),\underset{y\in Y}{\text{max}}\underset{x\in X}{\text{min}}d(x,y)\right\}.$$

(8)

The Hausdorff distance measures the maximum distance of a point in a set X to the nearest point in Y or vice versa. Generalization to 3D uses mean or median over all slices. The measure indicates how much manual corrections are required. An ideal value of HD equal to zero reflects complete agreement of the contours. In addition, the number of slices with a HD greater than 3 mm is reported and compared to the total number of slices with prostate information. A threshold of 3 mm is often seen as clinically acceptable [22]. A second criteria is the volume overlap (or Dice similarity coefficient) defined by

$$d_{\mathit{\text{vol}}}(X,Y)=2\frac{|X\cap Y|}{\left|X\right|+\left|Y\right|},$$

(9)

where |·| is the number of voxels contained in a region. Finally, accuracy is defined as

$$d_{\mathit{\text{acc}}}(X,Y)=1-\frac{\left|\mathit{\text{FP}}\right|+\left|\mathit{\text{FN}}\right|}{\left|\mathit{\text{TP}}\right|+\left|\mathit{\text{FN}}\right|}$$

(10)

with T P = X ∩ Y volume included in both X and Y (true positive), FN = X ∩ (¬Y) volume of X not included in Y (false negative) and FP = (¬X) ∩ Y volume of Y not included by X (false positive). Both values range from 0 to 1, with optimal value 1. Volume overlap indicates how much of the prostate has been detected by the approach while accuracy shows how incapable the method is to select the true prostate pixels.

In addition to the quantitative metrics, we have performed a small pilot test on 8 new patients comparing time expenditure using the proposed method and manual delineation. The time expenditure for the proposed method includes the definition of control points and the correction of the contour obtained by the method for each patient. The used mean shape model and variance was calculated from the training data as described above. Time measurements were obtained by two independent physicians for each case. Manual delineations, definition of control points and corrections were performed using the treatment planning system Eclipse^TM.

Results and discussion

The evaluation is performed on the deformed mean best fitting ellipses, i.e., on π (ρ ∣ CP) in formula (7). Additional evaluations are done in Additional file 1 for π (ρ ∣ S,CP).

Table 1 contains the distance metrics defined in (8) - (10) comparing the manual delineation and BFE for each test data set, and comparing the manual delineation and the deformed MBFE described by π (ρ ∣ CP). The high Dice similarity coefficient and accuracy values and small Hausdorff distances between manual delineations and BFE confirm the stacked ellipses model. This is also reflected by the count of slices with a HD greater than 3 mm in each test data set. Only 7.6% of the in total 105 slices have a HD greater than 3 mm. The values show the best possible description of the test cases by the proposed model. Furthermore, Table 1 presents the metrics comparing the manual delineation and the deformed MBFE for each test data set. The distance metrics reveal the fairly accurate results. The values indicate that the deformed MBFEs used as initial contours for final delineations will lower the time expenditure of the delineation procedure. However, more slices with a HD greater than 3 mm can be observed, 59% of the 105 slices. Particularly, the test sets 3, 7 and 10 have a ratio greater than 75%. Figures 1a to 1e illustrate 5 slices of the BFE evaluation of test patient 3 from Table 1. Figures 1f to 1j illustrate 5 slices of the deformed MBFE evaluation of test patient 3 from Table 1 with a volume overlap of the manual delineation line and the deformed mean shape of 0.90 and accuracy 0.81.

Table 1 Evaluation metrics comparing BFE to manual delineations, and comparing deformed MBFE to manual delineations

Table 2 shows the median and median absolute deviation (MAD) for the data groups “test data”, “training data” and “all data”. A BFE volume overlap for all data of 0.954 ± 0.010 (median ±MAD) and accuracy of 0.908 ± 0.020 confirm the model further (Table 2). Similar values in the subset of test data and training data are indicating model robustness. In addition to the BFE results, Table 2 summarizes the results by median and MAD of the distances between manual delineation and the deformed MBFE. A median volume overlap of 0.899 ± 0.021 and accuracy of 0.807 ± 0.035 of the test data show further the power of the prior. The deformation of the prior is done by the control points and can be computed directly since there is no sampling or estimation involved at this point.To evaluate the robustness of the model, we randomly split 10-times the set of 33 patients into a training set with 23 cases and a test set with 10 cases. Figure 4 shows the evaluation distances between the manual delineation and the deformed MBFE. The central mark is the median, the edges of the box are the 25th and 75th percentiles and the whiskers extend to the extreme data points. The figure shows only small variation between the different permutations, thereby demonstrating robustness of the stacked ellipses model.

Table 2 Evaluation metrics comparing BFE/MBFE and manual delineations for different data groups

Figure 4

Evaluation results between MBFE and manual delineations for 10 random permutations in training and test sets consisting of 23 and 10 patients, respectively. The Hausdorff distance is depicted in (a), the percent-volume overlap measured by the Dice coefficient in (b) and the accuracy in (c). The central mark is the median, the edges of the box are the 25th and 75th percentiles and the whiskers extend to the extreme data points.

Results from using MCMC to further optimize the delineation, as described in the Methods section, are only presented in Additional file 1 since a slight improvement comes at the cost of large computation time.

The time comparisons indicated an average of 30% time saving using the proposed method compared to manual delineation. The time measurement of the proposed method includes the definition of controls points as well as the correction of the estimated contour by the physician.

Conclusions

The presented results demonstrate the potential of the proposed method in modeling the prostate by slicewise best fitting ellipses. Deformation of the mean shape using control points gives very good results with little computational cost. Hence we believe that providing physicians with a good initial contour is beneficial in the clinical praxis of radiotherapy treatment.

The corrections of generated delineations based on few control points were not streamlined in the workflow of the physicians, and the task of correcting contours is not part of their everyday activity. Furthermore, corrections were not done directly after the definition of the control points and sometimes by different physicians, and physicians had to deal with a different orientation of the data set in the treatment planning system than in the diagnostic MRI. These issues must and can be solved for a well designed system. Therefore, a time saving of 30% likely represent a lower limit, and has to be validated in a well designed and properly powered study. Furthermore, we expect larger time savings in data sets where the prostate is imaged in a higher number of slices. In the extreme case, if the prostate is visible in only three slices, the BFE approach would not give any benefit using the current control point method. The study of inter/intra-observer variability using the proposed method compared to manual delineation was considered to be beyond the scope of this study and is left open for interesting future work.

In addition, the results show a precise description of the prostate by the BFE model with an average volume overlap of 95%. The high performance of the deformed mean shape model using the control points explains the small improvement by applying MCMC. Nevertheless, an improvement of the likelihood in the posterior distribution or by an active appearance model [23] is a field of further research as elaborated in Additional file1. A clear disadvantage of an additional method like MCMC is the extra computation time.

Further improvements can be achieved in the constraint and regularization terms, e.g., by considering the surface curvature versus changes of the ellipses parameters. We do not expect abrupt changes between neighboring slices around the central slice, but larger changes between slices towards the ends can be permitted, particularly in the length of the first and second principal axis. Also, the reduction of manual interaction in the proposed method is left for future work.