Digital phantoms for characterizing inconsistencies among radiomics extraction toolboxes

Medical imaging research labs are rapidly creating quantitative image analysis techniques that can potentially achieve expert human performance or aid clinical decision making. In particular, radiomics is emerging as one of the leading techniques for decision support, which is a high-dimensional extraction of quantitative features from medical images (i.e., CT, MRI, PET, etc) (Lambin et al 2012). These quantitative imaging features (i.e., radiomic features) can be implemented to develop predictive models for various clinical endpoints including overall survival (Parmar et al 2015a, Parmar et al 2015b), tumor staging (Liang et al 2016), tumor histology classification (Wu et al 2016, Zhu et al 2018, Lafata et al 2018), cancer recurrence (Saha et al 2018, Lafata et al 2019a), functional information (Lafata et al 2019b) etc Inherently based on standard-of-care imaging, radiomics may eventually be fully developed into a non-invasive, low-cost, and personalized decision support tool within clinical environments.

While radiomics is a promising technique, clinical translation is limited due to considerable variation in research pipelines. Typical radiomics pipelines consist of several essential steps, including image acquisition, image segmentation, feature extraction, and feature analysis and/or model building. Often, detailed descriptions of these key steps are missing (e.g., image acquisition parameters (Gilles, 2016), modeling parameters, radiomics extraction software etc), leading to a lack of standardization (Traverso et al 2018). This is a primary reason for the non-reproducibility of radiomics research. Feature extraction consistency is a particularly challenging issue. For example, radiomics findings have been published using various open source software packages, such as CERR (Computational Environment for Radiological Research) (Apte et al 2018), IBEX (imaging biomarker explorer)(Zhang et al 2015), CGITA (Chang-Gung Image Texture Analysis) (Fang et al 2014), and Pyradiomics (Van Griethuysen et al 2017). Furthermore, other groups have opted to develop their own in-house radiomics extraction toolboxes (Parmar et al 2015a, Lafata et al 2018, Larue et al 2017, Leger et al 2017). There is currently a lack of understanding regarding the radiomics extraction consistency among these different toolboxes, and possible variation may be in part responsible for discrepancies in research findings. For example, Foy et al (Foy et al 2018a) demonstrated that features extracted from an in-house toolbox and IBEX resulted in significant differences in classifying patients with and without radiation pneumonitis.

In general, the radiomics community has certainly been aware of the importance of benchmarking feature extraction techniques. Lambin et al (Lambin et al 2017) provided reference feature values extracted from pre-processed CT data of four patients with lung cancer. Foy et al. (Foy et al 2018b) compared features extracted from 2 in-house platforms and 2 open source radiomics toolboxes based on mammography scans and head and neck cancer CT images, and they found feature value discrepancies for texture features by nonparametric Friedman tests. Both studies used patient data. While using patient data provides realistic scenarios for radiomics studies, it is not intuitive for seeking out the source of feature extraction variation (such as texture matrix computation, scalar feature definitions, normalization techniques, etc). Given this consideration, digital phantoms are more suitable for benchmarking the feature extraction techniques. For example, Zwanenburg et al. initialized the Imaging Biomarker Standardization Initiative (IBSI) (Zwanenburg, 2017) to standardize feature extraction across institutions. They built a 5 × 4 × 4 voxel digital phantom with several discrete grey levels, from which standard feature values were computed and shared to 19 institutions in 8 countries. The institutions then reported the feature value similarities and differences.

Due to the current lack of standardization regarding feature extraction, different toolboxes may produce different—yet reasonable—feature values. More specifically, different feature values could be caused by different mathematical definitions, pre-processing implementations, scalar feature computation methods, etc, but these approaches may all be reasonable. Therefore, comparison of absolute feature values may not be fully informative. Instead, we propose a new method based on a series of digital phantoms to quantify relative inconsistencies and characterize sources of variation among several radiomics toolboxes. In particular, in this study we develop a series of digital bar phantoms and a series of heteromorphic sphere phantoms, from which we characterize feature extraction inconsistency among three radiomics toolboxes: IBEX, CERR and an in-house radiomics toolbox.

The overall workflow of this study as shown in figure 1 includes: phantom design (bar phantoms and heteromorphic sphere phantoms), feature extraction using three independent calculation toolboxes, and evaluation on feature extraction consistency.

Zoom In Zoom Out Reset image size

2.1. Phantom design

We designed two types of phantoms: heteromorphic sphere phantoms and bar phantoms. There were several considerations during the phantom design: (1) the feature values of the phantom should be easy to be manual verification or intuitively perceived for the convenient final validation; (2) no pre-processing (e.g. segmentation and grey level discretization) is necessary; (3) the phantom should be comprised of a series of phantoms with regular intensity, texture, and shape variations, for the convenience of characterizing the sources of variation.

2.1.1. Heteromorphic sphere phantoms

We constructed a series of heteromorphic sphere phantoms for characterizing shape features in a stochastic and iterative mechanism. The rational for this approach was to generate local topological changes on the surface of the phantoms (e.g., as shown in figure 2(a)) and study their effect on shape feature calculation. The starting point was a perfect sphere with radius of 80 voxels, which was iteratively perturbed to simulate changes in size and structure relative to baseline. In each iteration, a sphere with ten-voxel radius was stochastically attached to the surface of the previous iteration of the phantom. The voxel size of 1 mm × 1 mm × 1 mm was used in this study. In addition to the starting point, we acquired the phantoms every 500 iterations and up to 5500 iterations.

Zoom In Zoom Out Reset image size

2.1.2. Bar phantoms

As shown in figure 2(b), a series of bar phantoms were developed for characterizing intensity and texture features. The bar phantoms consisted of 66 × 66 × 12 voxels, with the voxel size of 1 mm × 1 mm × 2 mm. The outermost one-voxel shell was of value 0. The region of interest (ROI) was specified as the non-zero regions in the phantoms. The shell was designed to mimic the real environment in common radiomics studies. The phantoms were evenly distributed into 2, 4, 8, or 64 bars, and the grey levels of the n bars were evenly distributed between 1 and 64. The construction parameters of each phantom are shown in appendix A, Table A1.

2.2. Feature extraction

2.3. Evaluation on feature extraction consistency

3.1. Shape features

The variations of shape features calculated from the heteromorphic sphere phantom with iteration number are shown in figure 3. The volume values calculated from the four calculation algorithms were comparable, while the surface area values increased at two different rates. Since all other shape features are correlated with volume and surface area, they also varied at two rates with iteration number.

Zoom In Zoom Out Reset image size

The CCC and PCC for each shape feature are provided in appendix D, Table D1 and D2. The CCC confirmed the trends shown in figure 3: the four methods produced similar volume (CCC > 0.98); surface area was highly reproducible (CCC > 0.99) between CERR_D/in-house pair or CERR_M/IBEX pair yet weakly reproducible (CCC < 0.5) between the two groups. Only 2 features, the compactness 2 between CERR_M /IBEX and volume between CERR_M/CERR_D, demonstrated perfect agreement. Meanwhile, all four methods produced PCC > 0.8 for all shape features except compactness 1. For this feature, in-house/IBEX and CERR_M /IBEX resulted in PCCs of 0.76 and 0.77, respectively. These results demonstrate that inconsistencies existed in shape feature calculations among the four methods, but the features values were usually highly linear-correlated.

3.2. Intensity and texture features

In this study, digital phantoms were developed and implemented to investigate potential inconsistencies among three different radiomics feature extraction toolboxes: CERR, IBEX, and an in-house toolbox. While all three toolboxes incorporate similar stages of computation (i.e., DICOM reading, image pre-processing, contour delineation, gray-scale re-binning, etc), the fashion in which they are implemented are often quite different. Our results demonstrate that variations exist among radiomic features derived from different toolboxes.

As mentioned above, the IBSI (Zwanenburg 2017, Zwanenburg et al 2016) greatly contributed to benchmarking radiomic extraction processes across multiple institutions. Taking a different approach to the sequential phantoms developed in our study, IBSI utilized a single cubic-like digital phantom with several discrete gray levels to achieve its goal of benchmarking. Compared to a single phantom implementation, sequential phantoms have several unique advantages: (1) they are convenient for perceptive checking, which is especially meaningful for first-order features and shape features (relative to a perfect sphere). (2) They make the task of figuring out the sources of variations relatively easy, which would be helpful for the developers of radiomics toolboxes to debug their programs. (3) They provide more information during referring other radiomics studies. Researchers could know which features could be with extraction discrepancies by the CCC evaluation, and by how much the inconsistencies are non-linear by the PCC evaluation. Beyond the merits of serial phantoms, the bar phantoms and heteromorphic sphere phantoms were designed in such shapes and textures on purpose. First, the phantoms were designed to enable easy, intuitive check on the features, considering the physical meanings and mathematical definitions of radiomics features. In particular, most shape features describe the similarity of between the ROI and a perfect sphere, hence we designed a sphere-based phantom with morphological changes. Many texture features are relevant to the homogeneity or heterogeneity of the ROI, and the computation of texture features are related to the rotational variance. Therefore, bar phantoms with a regular homogeneity change and texture variations along different axes of rotation were designed. Second, the phantoms were designed in a way that the differences will be apparent and meaningful trends will be shown. For example, the complex surface of the heteromorphic phantoms results in apparent differences caused by the different algorithms used for calculating surface area and volume. Another example is that any differences in texture feature matrix computation (e.g., GLCOM and GLRLM) will also be apparent with the bar phantoms. Third, they help in understanding mathematical interpretations of features, which is especially helpful for characterizing new feature engineering techniques. For example, our results (as shown in figure 3) demonstrated that the sphericity of the heteromorphic sphere phantoms should be maximized at a perfect sphere and should decrease as the morphology changes until saturating.

We provided two evaluation metrics- CCC and PCC. The importance of CCC and PCC depends on the specific task. For the purpose of QA, CCC is the most important evaluation criteria, since it evaluates data reproducibility. Noted that in general there is no absolute threshold indicating an acceptable degree of discrepancy, even though we scored our CCC results into different categories. PCC is informative during referring others' research findings, since researchers would want to know how much inconsistencies could exist if they used different toolboxes. If the features produced from two toolboxes do not match based on CCC but linearly correlated, they may potentially generate similar results with respect to downstream models, relative to non-linear feature inconsistences.

The use of sequential phantoms, specifically separate phantoms for different types of features, can therefore potentially provide useful information regarding feature extraction inconsistencies across toolboxes. In our study, the disagreement between feature values was driven by many different factors. First, disagreements were observed as a result of differences in mathematical definitions. For example, feature S1 (compactness1) is defined via equation (1) (Aerts et al 2014) in in-house toolbox and IBEX, yet via equation (2) (Zwanenburg et al 2016) in CERR toolbox. Feature I3 (kurtosis) by in-house toolbox and by CERR is in relationship (3). This is the reason why the I3 is in weak agreement in CCC (IH/CERR) but high Pearson correlation (PCC = 1). None of the mathematical definitions were incorrect, however, we can see the necessity of standardizing the mathematical definitions for the generalization of radiomics research findings.

$$\begin{eqnarray}&&\begin{array}{l}Compactness1\_inhouse\\ \,=\,Compactness1\_IBEX=\displaystyle \frac{V}{\sqrt{\pi }{A}^{\displaystyle \frac{2}{3}}}\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&Compactness1\_CERR=\displaystyle \frac{V}{\sqrt{\pi }{A}^{\displaystyle \frac{3}{2}}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&Kurtosis\_CERR=Kurtosis\_inhouse-3\end{eqnarray} \tag{ 3 }$$

Second, disagreements were observed due to different pre-processing steps inherent to different toolboxes. As shown in appendix F, Table F1–F6 features I1 (energy) from IBEX demonstrated weak agreement with the other two toolboxes. This is due to the specific gray-level re-binning method implemented in IBEX. Specifically, IBEX applies a re-binning algorithm of dividing the CT number 0-4096 into 25 bins for intensity features, except for entropy whose binning range and bin width are adjustable. In general, however, this type of disagreement could be avoided for clinical applications of radiomics research findings.

Lastly, disagreements were observed as a result of differences in calculation methods. For example, when calculating texture features from joint-probability matrices, the in-house toolbox applied a z-direction normalization factor when accounting for rotational invariance of anisotropic voxels, while CERR did not apply this correction. Another example is the calculation of surface area. As shown in figure 3, the surface area calculated from Delaunay algorithm is lower than that from marching cubes algorithm, and the difference increases with iterations. This difference could potentially be caused by the inherent characteristics of the algorithms, and they are sensitivity in deriving final values. For example, the Delaunay triangulation algorithm tends to incorrectly calculate the fine interfaces of the small attaching spheres, especially at high iteration. Therefore, while these algorithms resulted in similar volume measurements, other shape features that are formulated based on surface area and volume demonstrated larger differences at high iteration.

In addition to these findings, our results also demonstrated the significance of calculating texture features in a fully 3D fashion. In the radiomics toolboxes, feature values are calculated as the average of feature values along 13 unique axes of rotation. Figure 4 illustrates the value of feature 'GLRLM Run Length Percentage' along the 13 unique axes of rotation, where we can see the directional dependence of the feature. IBEX calculates GLRLM features only along 0° and 90° directions, which may cause inconsistencies due to not fully approximating a rotationally invariant system.

Zoom In Zoom Out Reset image size

We have demonstrated that there are many sources of variation in feature extraction. Although discrepancies exist, it is important to note that the different implementations are not necessarily incorrect, as long as the software is consistent when comparing results. Researchers should be aware of the differences and consider it when building radiomic models. Additionally, how these variations would affect the clinical decision support is very much a task-based problem and worth further exploration. As concluded in previous studies (Lafata et al 2018, Foy et al 2018a), variations in feature values have been shown to affect downstream modeling results. In particular, Foy et al (Foy et al 2018a) demonstrated that features from different radiomics toolboxes affected the model performance on predicting radiation pneumonitis. Therefore, for future clinical radiomic application, the standardization or QA of radiomics toolboxes will be necessary for the reproducibility and generalization of radiomics research findings. Based on our findings, we suggest a working group for the standardized of radiomic features and the development of a protocol for the QA of radiomic applications. Digital phantoms, similar to the ones used in this study, can potentially be used as a tool for this purpose. Once researchers agree on a universal feature extraction ontology to address the mentioned disagreements, feature values extracted from the digital phantoms may serve as a warranted reference calculation.

We acknowledge that the current research only studies the variability associated with a single component of a typical radiomic pipeline. Thus, no clinical outcome comparison was involved in this study. However, variability of radiomic features with respect to image acquisition parameters, imaging modalities, preprocessing steps, or modeling algorithms could also affect final clinical decision making. Eventually the community will need to understand how different sources of variations interact to understand their aggregate behavior.

We developed a series of heteromorphic sphere phantoms and a series of bar phantoms for checking the feature extraction consistency among three radiomics toolboxes: CERR, IBEX, and an in-house toolbox. We observed inconsistent feature values and discrepancies in mathematical definition, pre-processing implementation, and calculation methods among the three toolboxes, thus demonstrating the significance of radiomics software QA. The digital phantoms can be applied as useful tools for the QA.

13 searching directions