Segmentation Evaluation with Sparse Ground Truth Data: Simulating True Segmentations as Perfect/Imperfect as Those Generated by Humans

1.1. Background

Numerous 3D anatomy segmentation methods have emerged since the advent of tomographic imaging modalities in the 1970s. Early methods were based purely on the information available in the image to be segmented (Herman et al., 1979; Liu 1977). Since they did not harvest information available via anatomic priors, they needed ground truth (or reference) segmentations only for segmentation evaluation. Although such approaches continue to seek new frontiers, methods that exploit priors in various forms have emerged during the past 2-3 decades and have shown significant gain in segmentation robustness and accuracy. These later methods may be generically referred to as model-based since they employ some form of model to encode prior information. Such models include shape and geographic models (Cootes et al., 1995; Pizer et al., 2003; Shen et al., 2011; Staib et al., 1992; Udupa et al., 2014), atlases (Ashburner et al., 2009; Christensen et al., 1994; Chu et al., 2013; Gee et al., 1993; Shi et al., 2017), and deep neural network models (Agn et al., 2019; Cerrolaza et al., 2019; Drozdzal et al., 2018; Moeskops et al., 2016; Wang et al., 2019a). However, for these methods, fully annotated ground truth (GT) segmentations that capture the very variability over a human population of focus they purport to encode is of fundamental importance, some of them requiring a large number of such data sets for robust model building alone, not to mention for evaluation as well.

There are two main issues with generating GT segmentations: expense and imprecision. GT reference is most typically provided by manual (human expert) contouring of anatomical objects in medical imaging. Thus, generating large fully manual GT sets becomes impractical and expensive (Schipaanboord et al., 2018). Unsupervised and semi-supervised methods may be utilized that do not require fully annotated data sets which can ease somewhat the expense issue. However, accuracy and convergence in learning are superior with supervised learning than with semi-supervised or unsupervised methods (Park et al., 2019). Crowd-sourced non-expert annotation can be another solution to the expense issue. However, the second issue – imprecision – that naturally and commonly exists among expert annotations becomes more pronounced when non-expert annotations are employed. Several factors (Joskowicz et al., 2019) lead to imprecision in manually annotated GT segmentations including image quality issues, lack of standard ways of defining objects or variations in the interpretation of (pseudo) standards (when they exist), human subjectivity in interpreting boundaries in images, and institution- and application-specific vagaries in clinical contouring culture. The magnitude of these imprecisions is object-specific and application-specific. Small, non-compact, and sparse objects entail larger degrees of imprecision inversely proportionate to their size compared to large, well-defined, and compact objects. The Expectation-Maximization-based STAPLE (Warfield et al., 2004) framework and its extensions are a series of methods commonly used to generate consensus GT from multiple manual segmentations. However, unsupervised/semi-supervised methods and generating consensus segmentations deal with only one of the two key issues and not both simultaneously.

1.2. Related work

Numerous works have investigated solutions to the issues of expense and imprecision in generating GT segmentations. Bounding box (Rajchl et al., 2016) or partial annotation via partial slices or scribble strategics (Can et al., 2018; Çiçek et al., 2016; Koch et al., 2017) are used in semi-supervised learning based algorithms in medical image segmentation, while it has been pointed out that unsupervised and semi-supervised methods cannot surpass the accuracy of fully supervised methods (Rajchl et al., 2016; Papandreou et al., 2015; Wang et al., 2019b). Cost-effective annotation methods deal with the expense issue by active learning or interactive segmentation to iteratively improve performance of segmentation models (Tajbakhsh et al., 2020), where active learning methods select most representative and uncertain samples to generate as few GT segmentations as possible (Yang et al, 2017), and interactive segmentation methods modify the auto-segmentation manually and fine-tune the model in a sample-specific manner (Wang et al., 2018). They are not suitable to generate a large GT set which is needed in evaluation.

Machine learning and deep learning methods have been utilized to conduct segmentation evaluation and estimate evaluation metrics without ground truth. A SVM regressor (Kohlberger et al., 2012) trained by a space of shape and appearance features or regression neural networks trained directly from binary and intensity images (Robinson et al., 2018) or the difference images of intensity images reconstructed from binary segmentations (Zhou et al., 2020) are utilized to yield metrics without knowing GT delineation. Reverse testing (Bhaskaruni et al., 2018) also contributed to segmentation evaluation without GT in medical images. (Valindria et al., 2017) introduce a reverse classification accuracy (RCA) framework which took predicted segmentations on test samples as pseudo ground truth and train a classifier to reversely test on labeled training samples. The predicted segmentations are considered as of good quality if the trained classifier works well on at least part of the training samples. Extensions of RCA are also utilized to estimate Dice evaluation metric values in cardiac magnetic resonance (MR) images in (Robinson et al., 2017; Robinson et al., 2018; Robinson et al., 2019). Instead of directly generating metrics for a single method or a single segmentation, there are also unsupervised methods to compare effectiveness of different segmentation algorithms by computational statistical measures (Chabrier et al., 2006) or region-correlation matrix (Sikka and Deserno, 2010).

Besides the expense issue, the imprecision issue also naturally and commonly exists in the human-drawn ground truth segmentations (Joskowicz et al., 2019; Park et al., 2019; Sharp et al., 2014; Yang et al., 2017) and will lead to different segmentation results and metric values when the samples are evaluated on ground truth generated by different algorithms (Lampert et al., 2016) or annotated by segmenters different from those employed for training samples (Shwartzman et al., 2019). (Heller et al., 2018) investigated the influence the errors in labels may have on the segmentation quality. (Bø et al., 2017) demonstrated that even for radiologists, intra-segmenter differences still exist depending on their familiarity with the segmentation tool. To minimize inter-segmenter disagreements and generate more precise consensus ground truth, segmentations from multiple human segmenters are utilized by averaging (Cheplygina and Pluim, 2018), majority voting (Nowak et al., 2010; O’Neil et al., 2017), maximizing topological agreements (Yang and Choe, 2011), and Expectation-Maximization-based STAPLE (Warfield et al., 2004) method and its extensions (Gordon et al., 2009; Li et al., 2011; Schlesinger et al., 2017; Shwartzman et al., 2019) . Since generating expert-level high quality annotations is a time-consuming and expensive task, methods have been investigated to evaluate the quality of crowd-sourced non-expert annotations and fuse crowd-sourced labels (Gurari et al., 2015). (Cheplygina and Pluim, 2018) observed that although the agreement from crowd-sourced annotations are best when utilized in medical image analysis, the disagreement of segmenters is also informative. Inter-segmenter differences provide good reference for algorithm evaluation (Joskowicz et al., 2019; Park et al., 2019; Popović et al., 2017) and also can be utilized to estimate regions with variability (Zhou et al., 2020) or uncertainty (Jungo et al., 2018) for helping in making clinical decisions.

Expense is purely a labor/cost issue. Imprecision, however, raises several conceptual issues. Although great strides have been made in examining these dual issues in the literature separately as delineated above, they have not been examined jointly or one as a function of the other. This area calls for a lot more attention in view of the promises suggested by deep neural network models. Most importantly, the practical question of the savings that ensue in cost as a function of the imprecision in GT data as a result of its “sparsilication” has not been examined so far. In other words, is it feasible to simulate full GT segmentations from sparse GT data such that the simulated GT is as perfect/ imperfect as, but not worse than, the GT generated by human experts? The cost saving then will directly depend on the degree of sparsity affordable for the sparse GT data and will be tied with the second issue of imprecision. In this work, keeping the expense issue in mind and recognizing the fact that perfect ground truth does not exist, we take inter-segmenter variability as reference to create segmentations as perfect/imperfect as those produced by expert segmenters with the attendant reduction of manual workload via proper strategies for sparse slice selection and segmentation filling. The created segmentations should be able to replace full manual GT in segmentation evaluation. There is no work in the published literature to address this important and very practical issue.

1.3. Outline of approach

In this paper, we propose a method, named SparseGT, to deal with the expense and imprecision issues jointly by exploiting natural variability existing among human segmenters, generating pseudo ground truth (p-GT for short) from sparse manual segmentations, and utilizing p-GT in actual segmentation evaluation. The method is fully described in Section 2 and is schematically presented in Figure 1. p-GT segmentation is generated in two sequential steps: selecting slices sparsely to conduct manual annotation and filling segmentations for other slices.

Each of the two steps is investigated by two strategies, including equal interval based uniform strategy and shape change function based non-uniform strategy for slice selection, and shape-based interpolation strategy (Raya and Udupa, 1990) and object-specific 2D U-net (Ronneberger et al., 2015) based deep learning strategy for segmentation filling. The largest sparseness levels where p-GT segmentations for the object under consideration are statistically indistinguishable from the full manual GT, taking into account the variability existing among a group of expert segmenters, are considered as the object-specific optimal sparseness to determine the cost saving that is feasible for the SparseGT method. The created p-GT can generate accurate metric values as full manual GT did in segmentation evaluation.

In Section 3, we describe the experimental set up, anatomic objects and data sets utilized, results, and their analysis. Experiments are conducted on the four possible combinations of strategies for sparse slice selection and segmentation filling. The optimal sparseness factors are object-, strategy- and application-specific, and the optimal combinations of strategics are determined by examining the yielded actual workload reduction and their actual availability in practice. Our conclusions, gaps remaining in this work, and avenues for potential improvements are discussed in Section 4.

Notation:

$\mathcal{B}$: Body region of interest.

O: Anatomic object of interest in $\mathcal{B}$.

$\mathcal{I}=\{I_1,\dots ,I_N\}$: Image data sets of $\mathcal{B}$ for which complete GT delineations for O are available.

${\mathcal{I}}_b=\{I_{1,b},\dots ,I_{N,b}\}$: Binary images representing complete GT delineations of O in $\mathcal{I}$.

${\mathcal{J}}_b=\{J_{1,b},\dots ,J_{N,b}\}$: Binary images representing complete segmentations of O in $\mathcal{I}$ by a segmentation algorithm A.

${\mathcal{P}}_b=\{P_{1,b},\dots ,P_{N,b}\}$: Optimally generated p-GT corresponding to ${\mathcal{I}}_b$.

Abbreviations and acronyms commonly used in this paper are listed in Table 1 for ease of reference.

2.0. Determining slice range for object O

Let N_O be the number of slices covering an object O (such as the mandible) in a patient image that includes O. The coverage of O should be decided by expert segmenters to guarantee that a proper and consistent anatomic definition of O is adhered to. Without loss of generality, we assume that the slices are transaxial with respect to $\mathcal{B}$ and denote the cranio-caudal direction orthogonal to the slice plane by z.

2.1. Sparse slice selection

One of the main aims of the SparseGT method is to reduce manual workload needed in creating expert-quality p-GT for segmentation algorithm evaluation, where only a few out of all slices are selected to conduct manual annotation and segmentations on other slices are automatically created via segmentation filling strategies. Two strategies are investigated to determine sparse slice positions, as illustrated in Figure 2. Along the z axis, slices of object O are sparsely selected in a uniform or non-uniform manner.

Step 1. Estimate a shape change measure M_SC(s_i) for each slice position.

Let s_i denote a slice through object O in a GT binary image I in ${\mathcal{I}}_b$ at position z_i. To estimate the change in shape of O in I at s_i, a neighborhood is first defined around s_i to measure local shape change at s_i. A factor R is used to define this neighborhood, where N_n = max(⌊N_O/R⌋, 1) slices previous and next to s_i are considered in defining M_SC(s_i). The shape change is measured by two metrics: exclusive or (EOR) and Jaccard index (JI). M_SC(s_i) is calculated as a weighted average of shape changes between s_i and each slice in its neighborhood, where the weight factor is determined based on the distance of the neighboring slices from s_i, as shown in (1):

where m(s_i, s_j) stands for the shape change measured by one of the metrics (EOR or JI) between slices s_i and s_j at locations z_i and z_j, and the nearer neighboring slices are assigned larger weights.

Step 2. Normalize slice positions and M_SC(s_i) both to range [0, 1].

Since different object samples are likely to contain different numbers of slices and different magnitudes for M_SC(s_i), to devise a standardized shape-change function that avoids interpolation of images, the slice positions of different object samples are normalized to a fixed slice range and the magnitudes of M_SC(s_i) are also normalized. A sufficiently large parameter N_R is determined such that N_R > N_O for all considered object samples, and the normalized range [0, 1] is discretized into N_R positions with N_R − 1 intervals. Then, all slices of the object samples are assigned with proportional discrete positions, where only the slice positions are interpolated into a unified range and the actual image slices are not changed. We denote the normalized positions by the variable $\overline{\text{z}}$. Besides slice positions, M_SC(s_i) also shows different magnitude ranges in different object samples irrespective of whether the shape change is measured based on EOR or JI. To facilitate deriving the mean shape-change function ${\overline{\text{δ}}}_{\text{S}}(\overline{\text{z}})$, we also normalize M_SC(s_i) by min-max normalization:

where M_max and M_min denote the maximum and minimum shape change of each object sample.

Step 3. Estimate standardized mean shape-change function ${\overline{\text{δ}}}_{\text{S}}(\overline{\text{z}})$.

The value of the shape-change function ${\text{δ}}_{\text{S}}(\overline{\text{z}})$ at each discrete position $\overline{\text{z}}$ for each object sample in a population of object samples S is estimated as in (3). The mean shape-change function ${\overline{\text{δ}}}_{\text{S}}(\overline{\text{z}})$ is then derived by averaging ${\text{δ}}_{\text{S}}(\overline{\text{z}})$ over the population S. It is possible that some of the positions on the whole range of $\overline{\text{z}}$ are not accounted for by data from any of the subjects, so values at those positions are set as null. This will not influence the modelling of the average shape-change function ${\overline{\text{δ}}}_{\text{S}}(\overline{\text{z}})$.

(1). Shape-based interpolation (SI) strategy

Given two adjacent sparse slices l₁ and l₂ with manual annotations, the segmentation on a slice l_i at any position in between is estimated by using the shape-based interpolation approach (Raya and Udupa, 1990) by following the slice-to-slice shape of the object. Signed (2D) distance transform is first applied to l₁ and l₂ with the convention that distances from the object boundary are positive for pixels inside the object and negative for outside pixels. Then, the distance map for l_i is estimated by linearly interpolating the distance maps for l₁ and l₂, and the binary mask is finally determined by pixels with positive distance values on slice l_i.

(2). Deep-learning (DL) based strategy

The deep-learning based approach is developed using the 2D U-net architecture (Ronneberger et al., 2015), where the input and output images are designed in different manners for the two sparse slice selection strategies. For the t-based uniform selection strategy, as shown in Figure 3(a), the range of continuous four sparsely selected slices contains a total of 3t + 4 slice positions with four selected sparse slices and three blocks of non-selected slices in between. The intensity image slices and binary slices are combined into 2 × (3t + 4) = 6t + 8 channel images as the input to the network, where only the selected sparse slices contain annotated segmentations while the binary image slices on non-selected positions are all empty or with 0 value. The training target is to estimate GT segmentations for the central block of non-selected t slices, and the loss is calculated between predicted GT results and actual GT on the central block of slices. The other two blocks are for giving contextual information.

For the k-bascd non-uniform strategy, we will follow the above spirit of including information from four consecutive selected sparse slices to design a 10-channel input as illustrated in Figures 3(b) and (c). M_a, M_b, M_c, and M_d are four consecutive sparse slice positions. When filling segmentations in the skipped block between slices at positions M_b and M_c, intensity and binary images on M_a, M_b, F₁, M_c, and M_d compose the 10-channel input at first, where F₁ denotes the first skipped slice position next to M_b. Subsequently, after the slice at F₁ is filled with segmentation predicted by the network in Figure 3(b), the input to the network is updated with the composition of images at positions M_a, F₁, F₂, M_c, and M_d, where F₂ is the target slice next to F₁. The second slice position is always the slice previous to the target slice with manual or predicted binary mask. This process continues iteratively until all skipped slices in this block are filled with segmentations.

The networks are composed of convolutional layers with 3 × 3 kernels followed by Batch Normalization and ReLU nonlinear activation, and the output layers are activated by a sigmoid function. Down-sampling is achieved by convolutional layers with stride 2 while others are with stride 1. Adam optimizer is used to minimize the crossentropy loss function in training the networks.

3.1. Data sets and experiments

This study was conducted following approval from the Institutional Review Board at the Hospital of the University of Pennsylvania along with a Health Insurance Portability and Accountability Act waiver. Experiments are conducted on computed tomography (CT) images of two body regions, Head & Neck and Thorax, with three objects in each body region – cervical esophagus (CtEs), mandible (Mnd), and orohypopharynx constrictor muscle (OHPh) in the Head & Neck body region, and heart (Hrt), trachea & proximal bronchi (TB), and right lung (RLg) in the Thorax body region. The objects are chosen to represent different shape and size characteristics and different degrees of challenges for segmentation (RLg and Mnd: least challenging, Hrt and TB: moderately challenging, CtEs and OHPh: most challenging). Our set $\mathcal{I}$ of images consists of 498 3D images – 298 in Head & Neck region and 200 in Thorax region. For the above 6 objects, GT data are available to us which constitute real clinical data as contoured by several dosimetrists (as explained in the previous section) for the routine RT planning of patients with Head & Neck or Thoracic cancer. For 81 of the Head & Neck studies, each of two dosimetrists contoured the above 3 objects separately. Similarly, another two dosimetrists contoured the 3 Thoracic objects separately on 87 Thoracic studies. These data sets constitute the training set we denoted by V previously and were utilized for estimating (μ_M, σ_M).

We note that not every study in the cohort $\mathcal{I}$ necessarily has all 6 objects contoured since the actual objects contoured in clinical practice for RT planning depend on the location and size of the tumor in each study. The data sets in the two body regions are separately divided into disjoint sets of training and test samples, where about 20% samples are randomly selected to compose the test set Te, and the remaining samples compose the training set Tr to determine optimal sparseness factors t_O and k_O, which are verified on Te to check if the test set can also yield indistinguishable deviation with respect to expert variability. The voxel size in our data sets varies from 0.93 × 0.93 × 1.5 mm³ to 1.6 × 1.6 × 3 mm³. The object sizes are also variable among subjects, where object samples of different subjects occupy varying numbers of slices. The bounding box size is determined for each object based on the largest occupied range among its samples. To fit the input size to the DL network, the 2D region of interest (ROI) of object samples is set with size in multiples of 8, since there are three convolutional layers with stride 2 in both of the designed U-net based networks. The splitting of data sets and their roles in different stages of SparseGT are listed in Table 3 for clarity, and key information about data sets is summarized in Table 4 for ready reference.

3.2. Results and discussion