A semiautomatic tool for prostate segmentation in radiotherapy treatment planning

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 http://www.biomedcentral.com/1471-2342/14 /4 approaches are available in the literature, such as the medial or skeleton representation of the prostate [4][5][6][7][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 slicewise 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][17][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.
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.

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 1 , i 2 ) ∈ {1, . . . , I 1 } × {1, . . . , I 2 } are the pixel indices in a slice, where typically I 1 = I 2 , 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 1 by I n1 and I 2 by I n2 , but for simplicity we use H, I 1 and I 2 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 , . . . , X nK n ) with X nk = (x n 1k , x n 2k , x n 3k ) T ∈ R 3 , k = 1, . . . , 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 http://www.biomedcentral.com/1471-2342/14/4 in a sequentially sorted number L n of equidistant slices whereas each contour slice containsK nl contour points, l = 1, . . . , L n . Hence it follows K n = lK nl and X n = (X n1 , . . . ,X nL n ). The image information in slice l is denoted by S nl and S n = {S n1 , . . . , S nL n } ⊆ V n andX nl defines the configuration matrix in slice S nl .
In summary, the training population is given by the set {V, X}, with a set of volume information V = {V 1 , . . . , V N } and configuration matrices X = {X 1 , . . . , 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 1 , i 2 , 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 derotated PCS where volume and contour information are aligned to each other.

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, 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 n dCP maps the de-rotated data {dBFE n , dCP n } to {BFE n , CP n } in the sample space. http://www.biomedcentral.com/1471-2342/14/4 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 The rotation parameter φ nl is defined corresponding to the ICS with origin θ nl in slice S nl . The boundary of an ellipse ρ nl centered at θ nl ∈ R 2 in slice S nl is defined by The shape model described in this section requires the best fit of an ellipse C(ρ nl ) to the contour informationX nl in each slice, i.e., we modelX nl = C(ρ nl ) + 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 forX nl . The nonlinear estimate of parameters ρ nl = (θ nl 1 , θ nl 2 , α nl 1 , α nl 2 , φ nl ) T givenX nl must minimize the error whereC(ρ nl ) is a set of nearest orthogonal points ofX nl to C(ρ nl ).

Definition 1 (Best fitting ellipse (BFE)). A best fitting ellipse for slice S nl is defined by the set BFE nl
. . , L n , n = 1, . . . , N and minimizes the error function g, i.e., BFE nl =ρ nl with The first and second principal axes must be reordered after calculation of BFE n = {BFE n1 , . . . , BFE 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(|φ i − φ M |) = π, i ∈ {1, . . . , L} after reordering. Therefore, the set BFE n of reordered bestfitting ellipses is an element of (R 2 × R 2 + × (−π, π] ) 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 φ nl in case of circularity. If both principal axes have the same length, the orientation of an ellipse is undefined. Therefore we penalize φ nl in the case of high circularity by taking φ nl from the neighboring slices into account. Second, smoothing is performed between neighboring slices to avoid large forward and backwards rotations between φ n(l−1) , φ nl and φ 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 BFE n ∈ (R 2 × R 2 + × (−π, π] ) 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 n dCP . The transformation matrix n dCP maps the de-rotated prior data {dBFE n , dCP n } to {BFE n , 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., CP n = A n 1 , . . . , A n 6 , P n 1 , . . . , P n 6 , B n 1 , . . . , B n 6 as visualized in Figure 3.
In addition, we have tested alternative control point configurations. They are described together with the construction of n dCP in Additional file 1.
After transformation we have obtained a reordered and comparable set of best fitting ellipses The statistical analysis of the training data requires an equal number L 1 = . . . = L N to establish correspondence between the parameters of the best fitting ellipses. Therefore, we interpolate the set BFE 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 θ 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 ) T for 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 (R 2 × R 2 + × (−π, π] ) L .

Statistical analyses
After reconstruction of our shape space we estimate the expectation and variance of the parameters of a mean shape model μ BFE = {μ 1 BFE , . . . , μ L BFE } with μ l BFE = (μ l θ , μ l α , μ 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 θ nl = (θ nl 1 , θ nl 2 , θ nl 3 ) T in terms of a distance vector η nl of θ nl to a center curve defined by the control points. We model θ l = ξ l + η l , where ξ l is 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 where μ l θ = (μ l θ 1 , μ l θ 2 , μ l θ 3 ) T , l = 1, . . . , L. The variance and covariance are estimated by The length parameter is modeled by a log-normal distribution because α ∈ R 2 + . Thus we estimate the mean and variance of a = log(α) ∈ 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 If θ l is defined according to the center curve given by the control points as described above, we model 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 http://www.biomedcentral.com/1471-2342/14/4 is obtained through the posterior π(ρ | S) where s il ∈ S ⊆ V is the volume information and i = (i 1 , i 2 ) ∈ 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 defines the posterior density of the deformed template π(ρ | CP) 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 π(ρ | CP) models realistic variations from our mean shape μ BFE ∈ (R 2 × R 2 + × (−π, π]) 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 T1weighted 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 μ L = 1 N N n=1 L n and the standard deviation by 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 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 where | · | is the number of voxels contained in a region. Finally, accuracy is defined as with TP = 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). http://www.biomedcentral.com/1471-2342/14/4  (unit: Dice 3D and accuracy in percentage, HD mean in mm, number of slices with HD ≥ 3mm compared to the total number of slices in brackets). 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 Table 1 with a volume overlap of the manual delineation line and the deformed mean shape of 0.90 and accuracy 0.81. 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.  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/intraobserver 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.

Additional file
Additional file 1: Supplementary materials. Article containing i.) a detailed description of the relative coordinates systems (e.g. ICS, PCS) on the basis of the DICOM file structure, ii.) post-processing procedures as for example reordering and introduction of constraints, iii.) a discussion of different control point method with construction of the transformation matrix n dCP and the parameter η nl , iv.) elaboration of the posterior distribution, and v.) a section with additional data analysis. http://www.biomedcentral.com/1471-2342/14/4 data set. Furthermore, JS drafted and revised the manuscript. SOS and FG contributed to the conception and design of the study, and helped to draft the manuscript. VKT and KM coordinated and carried out the acquisition of patient data. All authors read and approved the final manuscript.