The complexity of tumor shape, spiculatedness, correlates with tumor radiomic shape features

Radiomics extracts high-throughput quantitative data from medical images to contribute to precision medicine. Radiomic shape features have been shown to correlate with patient outcomes. However, how radiomic shape features vary in function of tumor complexity and tumor volume, as well as with method used for meshing and voxel resampling, remains unknown. The aims of this study are to create tumor models with varying degrees of complexity, or spiculatedness, and evaluate their relationship with quantitatively extracted shape features. Twenty-eight tumor models were mathematically created using spherical harmonics with the spiculatedness degree d being increased by increments of 3 (d = 11 to d = 92). Models were 3D printed with identical bases of 5 cm, imaged with a CT scanner with two different slice thicknesses, and semi-automatically delineated. Resampling of the resulting masks on a 1 × 1 × 1 mm3 grid was performed, and the voxel size of each model was then calculated to eliminate volume differences. Four MATLAB-based algorithms (isosurface (M1), isosurface filter (M2), isosurface remeshing (M3), and boundary (M4)) were used to extract nine 3D features (Volume, Surface area, Surface-to-volume, Compactness1, Compactness2, Compactness3, Spherical Disproportion, Sphericity and Fractional Concavity). To quantify the impact of 3D printing, acquisition, segmentation and meshing, features were computed directly from the stereolithography (STL) file format that was used for 3D printing, and compared to those computed. Changes in feature values between 0.6 and 2 mm slice acquisitions were also compared. Spearman’s rank-order correlation coefficients were computed to determine the relationship of each shape feature with spiculatedness for each of the four meshing algorithms. Percent changes were calculated between shape features extracted from the original and resampled contoured images to evaluate the influence of spatial resampling. Finally, the percent change in shape features when the volume was changed from 25% to 150% of their original volume was quantified for three distinct tumor models and compared to the percent change observed when modifying the spiculatedness of the model from d = 11 to d = 92. Values extracted using isosurface remeshing method are the closest to the STL reference ones, with mean differences less than 10.8% (Compactness2) for all features. Seven of the eight features had strong significant correlations with tumor model complexity irrespective of the meshing algorithm (r > 0.98, p < 10-4), with fractional concavity having the lowest correlation coefficient (r = 0.83, p < 10-4, M2). Comparisons of features extracted from the 0.6 and 2 mm slice thicknesses showed that mean differences were from 2.1% (Compactness3) to 12.7% (Compactness2) for the isosurface remeshing method. Resampling on a 1 × 1 × 1 mm3 grid resulted in between 1.3% (Compactness3) to 9.5% (Fractional Concavity) mean changes in feature values. Compactness2, Compactness3, Spherical Disproportion, Sphericity and Fractional Concavity were the features least affected by volume changes. Compactness1 had a 90.4% change with volume, which was greater than the change between the least and most spiculated models. This is the first methodological study that directly demonstrates the relationship of tumor spiculatedness with radiomic shape features, that also produced 3D tumor models, which may serve as reference phantoms for future radiomic studies. Surface Area, Surface-to-volume, and Spherical Disproportion had direct relationships with spiculatedness while the three formulas for Compactness, Sphericity and Fractional Concavity had inverse relationships. The features Compactness2, Compactness3, Spherical Disproportion, and Sphericity should be prioritized as these have minimal variations with volume changes, slice thickness and resampling.

Impact of the meshing algorithm on feature values. Pixel dimensions varied from 0.643 mm to 0.835 mm (X and Y directions) and from 0.565 to 0.733 mm (Z direction) after volume equalization for the masks extracted from the 0.6 mm slice thickness images. These values ranged between 0.639 mm and 0.840 mm (X and Y directions) and between 1.870 to 2.459 mm (Z direction) for the masks extracted from the 2 mm slice thickness images. Figure 2 shows the meshes obtained from the 0.6 mm slice thickness acquisitions for a representative model d = 47 for the four algorithms. Figure 3 illustrates the impact of the meshing algorithm and slice  www.nature.com/scientificreports www.nature.com/scientificreports/ thickness on feature values for all the 28 tumor models after volume equalization. Outliers are observed using the Boundary M4 method. Differences ranged from 4.1% (Compactness3, d = 77) to 43.2% (Compactness2, d = 11) between M1 and M2 methods. Surface area, Surface-to-volume, and Spherical Disproportion had direct relationships with spiculatedness (increasing value with increasing tumor spiculatedness). The three formulas for Compactness, Sphericity and Fractional Concavity had inverse relationships. All features exhibit large variations between d = 11 to d = 41. Compactness2 is able to highlight shape differences even for the least spiculated models. For this feature, no slope breaking is observed until d = 65 irrespective of the slice thickness and meshing method. Subsequent analyses were performed on the 2 mm slice thickness acquisitions, as this is more coherent with what is used in the clinics.
Correlations between tumor complexity and features. Spearman's rho correlation coefficients between each shape feature and tumor complexity were computed for the M1, M2 and M3 meshing algorithms ( Table 3). Seven of the eight features had strong significant correlations with tumor model complexity irrespective of the meshing algorithm (r > 0.98, p < 10 -4 ). Fractional Concavity showed the lowest correlation coefficient (r = 0.83, p < 10 -4 , M2).
Correlations between shape features. Correlations among the eight shape features, with the exception of volume that was fixed at a constant value, were calculated using Spearman's rho coefficients for the M1, M2 and M3 meshing algorithms. Almost all of the features were highly correlated with each other with r = 1, as seen in the correlation matrix plots (Supplementary Figure S1). Only Fractional Concavity was slightly less correlated with the others, with r values from 0.85 to 0.99, with M2 having the lowest correlation.
Effect of grid resampling. All Figure 4). The ranking of the models was not influenced by the resampling.  www.nature.com/scientificreports www.nature.com/scientificreports/ Impact of change in volume on shape features. The absolute percent changes of each feature with regards to the change in volume from 25% to 150% were compared to the percent change observed when modifying the spiculatedness of the model from d = 11 to d = 92 for M3 (Table 4). Compactness2, Compactness3, Spherical Disproportion, Sphericity and Fractional Concavity were less affected by volume changes. Surface-to-volume and Compactness1 were more affected by change in volume than tumor complexity, with the Surface-to-volume feature having a 69.9% change from the least to most spiculated models versus 54.2% change with volume; and Compactness1 having a 47.5% change with spiculatedness compared to 90.4% change with volume. Figure 5 illustrates the change in feature values with changes in volume for the three representative phantoms d = 11, 47, 92. Figure 6 summarizes the effect of the technical parameters slice thickness, resampling, and change in volume on the radiomic shape features, reiterating that Compactness1 has important variations with changes both in slice thickness and in volume, and Surface-to-volume with change in volume.
Comparison with actual tumors. To confirm the clinical relevance and applicability of our method, three representative patients from the RIDER database were contoured, and radiomic shape features were extracted thereafter. The range of values for each of the features among the three patients are within the range of the values from the shape phantoms, even if the volume after equalization was greater in the phantoms (range:

Discussion
This is the first methodological study that directly demonstrates the relationship of tumor spiculatedness with radiomic shape features. Models with increasing degrees of spiculatedness were created to examine the behavior of quantitative shape features with known incremental degrees of tumor border complexity. It was seen that specific features increase monotonically with increasing tumor spiculatedness, in particular Surface Area, Surface-to-volume, and Spherical Disproportion. Conversely, certain features exhibit a monotone decreasing correlation with increasing spiculatedness (Compactness, Sphericity, Fractional Concavity). Quantitative extracted shape features have already been demonstrated to give insights on tumor behavior, underlying their importance in radiomic analysis. Based on CT scans, several publications have shown that shape features differentiate between benign and malignant nodules [6][7][8] as well as correlate with patient outcomes 2,9,28 . In addition, a radiomic study of pre-treatment contrast-enhanced T1 MRI images in glioblastoma showed that tumor surface regularity was a powerful predictor of survival in the discovery (p = 0.005, hazard ratio [HR] = 1.61) and validation groups (p = 0.05, HR = 1.84) 29 .
It can be seen that many of the features exhibit strong correlations with each other, either positive or negative. If the behaviors of certain features are known to depend on specific parameters, calculating all may not necessarily give complementary information but instead redundant ones. In particular, Surface-to-volume, Compactness, Spherical Disproportion and Sphericity are all calculated from tumor volume and surface area 17 , which explains the strong relationships among them. In this regard, features might eventually be grouped into clusters instead of being analyzed individually 12 . For instance, different formulas for compactness have been previously published and used in radiomic studies 8,17 . In clinical studies, it needs to be determined whether correlations seen with radiomic shape features are inherent, or if tumor volume is a confounding factor. In our study, it is seen that Surface-to-volume and Compactness1 are affected with volume changes, and should thus be used with caution when comparing tumors with differing volumes. Indeed, even if it has been widely used in previous publications, the result obtained from the formula for Compactness1 is not dimensionless, and thus is not ideal in feature analysis. Compactness2, Compactness3, Spherical Disproportion, and Sphericity's percent changes between d = 11 and d = 92 were noticeably higher than the percent change with volume variations, which may make these features more useful in analyses of patient tumors as they are not volume dependent. In addition, it was also seen that in general, the features highlight differences in complexity better in more spiculated tumors. For instance, in this study, the slope of the relationship between feature and spiculatedness was steep until d = 41, and thereafter relatively flattened out for the less spiculated models.   www.nature.com/scientificreports www.nature.com/scientificreports/ In this study, four meshing algorithms have been used and the differences between the features computed directly from the STL files and those computed after meshing have been compared (Fig. 1). Values extracted using M3 meshing method are the closest to the STL reference ones, with mean differences lower than 10.8% (Compactness2) for all features. The decrease of the standard deviation when comparing the values from STL files versus from M1 and M3 with 0.6 or 2 mm slice thickness (Fig. 1) validates the fact that a 2 mm slice thickness should be preferred for shape-based radiomic analysis. With a 2 mm slice thickness, shape features are impacted in a more homogeneous way by the entire radiomics process. The associations of the deduced shape features with change in tumor complexity have also been analyzed for each method (Fig. 3). Using the MATLAB Boundary function, a non-monotone behavior of the features with phantom spiculatedness was observed, with the presence of outliers. The use of this function is thus not recommended in in-house MATLAB-based softwares. Comparison of the M1, M2 and M3 methods shows that different meshing implementations can lead to different quantitative values. As a consequence, thresholds determined in the literature should be used with caution. Numerical phantoms such as the ones developed in this study can be also of major interest for the evaluation of the pertinence of meshing algorithms as well as for the development of new shape features. Notable is that in this study, we chose to use meshing for volume extraction, which is not performed in most of the radiomic software that typically www.nature.com/scientificreports www.nature.com/scientificreports/ multiplies the voxel size by the number of voxels in the volume of interest 30 . This choice is of importance for maintaining consistency between surface and volume quantities.
Comparisons of feature values between scans acquired with 0.6 mm versus 2 mm slice thickness reveal that this parameter affects all radiomic shape features, with changes of up to 12.7% (Compactness1, M3). Resampling the CT images on a 1 × 1 × 1 mm 3 grid likewise resulted in small differences of between 1.3 to 9.6% changes in extracted features for the isosurface remesher M3 method. In a phantom study that computed the differences between original features and those resampled on a 1 × 1 × 2 mm 3 grid with original pixel sizes ranging from 0.39 to 0.98 mm, shape features belonged to the group that were generally not significantly affected by resampling 31 . In this study, the Credence Cartridge Radiomics (CCR) phantoms used were rectangular in shape and created primarily for texture analysis, whereas ours had fine spiculations specifically created for shape analysis. Another study using the same phantoms showed that radiomic features were affected by slice thickness, but that this effect could be reduced by resampling the images before feature extraction. However, this study focused on 114 first order and textural features and did not include shape 32 . In yet another phantom study using spherical, elliptical, lobulated and spiculated forms, it was shown that shape features were significantly different between 1.25 and 5 mm slice thickness scans 24 , from which we can infer that voxel size affects results of feature extraction. At present, we therefore recommend not to constitute a cohort with images having too different slice thicknesses particularly if the Compactness1, Surface-To-Volume and Fractional Concavity indices are computed, given their dependence on slice thickness. However, the ideal is prospective studies with homogenous acquisition parameters, as resampling alone does not completely eliminate bias resulting from differences in acquisition such as slice thickness.
There are disadvantages to this study. First, only 3D features were calculated as the tumor models were contoured on axial CT slices and had discontinuous islets on some slices (usually at the top and edges of the tumor) because of the spiculations. In addition, the tumor phantoms were printed with a flat base, instead of a spherical-based shape with no flat edges due to technical considerations for 3D printing. However, all the phantoms were created in the same manner (with a flat base) such that all shape feature variations are expected to be comparable. Another limitation is that although the shape phantoms have increasing degrees of complexity, the variations of these are all based on the formula of spherical harmonics and thus have a consistent mathematical progression. Actual tumors are rarely symmetric and regularly shaped. However, theoretical knowledge of how radiomic shape features vary remains of value in deducing the complexity of actual tumors. Also, in studying variations with volume, the volumes were modified mathematically by recomputing the pixel sizes, which are inherently correct; but another way would have been to do a 3D reprint of each model with each corresponding volume change. Another limitation in the conduct of radiomics studies in general is that there is no generally accepted and universally utilized meshing method, and as illustrated in this study, different methods do not result to identical values.
In summary, majority of radiomic shape features have strong monotone direct or inverse correlations with tumor spiculatedness. However, we have shown that quantitative values of these features can vary with slice thickness, volume, and resampling; and depend on the meshing algorithm used for surface and volume extraction. The radiomic shape features Compactness2, Compactness3, Spherical Disproportion, and Sphericity have been shown to have minimal variations with the aforementioned parameters, and should thus be prioritized in future studies. It is clear that quantitative radiomic shape features provide important information on tumor characteristics, underlining the importance of their integration into future radiomic models and notably their combination with clinical, textural and genomic features. Refinements in the methodology of conducting radiomic studies as well as transparency in the exact nomenclature and formula used for each feature are indispensable to enable its eventual translatability to clinical utility. Effects of the different parameters were compared to the ability of each feature to distinguish change in spiculatedness. Green cases correspond to a ratio of less than 5% between the effect of the technical parameter to the percent change observed when modifying the spiculatedness of the model from d = 11 to d = 92. Orange cases correspond to ratios ranging from 10 to 20% and red cases to ratios superior to 20%.

Material and Methods
shape phantoms. Spherical harmonics were used to create mathematical tumor models with increasing degrees of complexity 28 . In the presented work, ℓ was arbitrarily set to 10, m to 5 and A to 10. The degree of the spherical harmonic, d, was increased in increments of 3 from 11 to 92 to create a total of 28 models. Model "11" corresponded to the most spiculated model and model "92" to the least one. The 3D models were cut in the middle of the horizontal plane which permitted these to have a flat base for printing. Then, models were set with identical bases of 5 cm by adjusting the height ratio of the original models to the new base. Each of the tumor models was created using a 3D printer (Discoeasy200, dagoma.fr), using a polylactic acid filament (ρ = 1.25 g.cm -3 ) with standard printing speed (Supplementary Figure S3).
The models were then scanned using a Siemens Sensation Open CT scan (Siemens Healthineers, Erlangen, Germany) with 0.6 mm and 2 mm slice thicknesses, 100 kVp tube voltage, 300 mAs and a 350 mm reconstruction field of view. The phantoms were scanned on top of a cardboard box, with only the bases being in contact with a surface. Original pixel size was 0.68 mm × 0.68 mm in the transverse planes. Scans were contoured using the thresholding function of 3DSlicer (http://www.slicer.org), with lower limit at -700 Hounsfield Units (HU) and upper limit at 3000 HU, resulting in binary masks.  Table S4). The first method (M1) used the Isosurface function of MATLAB. This method connects points having the same value to generate the mask. The isovalue was set to 0.9. The second method (M2), Isosurface filter, consisted of smoothing the triangulated mesh generated with the first method by using the normalized curvature operators as weights. The mesh was mainly smoothed in the normal direction to preserve the original ratio in length between edges. One smoothing iteration was used and the smoothing quantity was set to 5. The third strategy (M3), Isosurface remesher, is an iterative triangle optimization for meshing. In this method, all the closed meshes obtained with the first method are cleaned according

Volume
Compute the enclosed volume of the object of interest. The enclosed volume is evaluated by triangulation (i.e. dividing the surface into connected triangles) Green-Ostrogradski formula: where F corresponds to the vector field deduced from the triangulation Surface area 17 Area of the surface encompassing the volume of interest, calculated by triangulation Surface-to-volume ratio 17 Ratio of surface to volume = Surface to volume ratio A V Compactness1 17 Describes how much the shape of a tumor resembles that of a sphere/can be encompassed by a sphere Compactness of a sphere = 1 = π .
Spherical disproportion 17 The ratio of the surface area of the tumor to the surface area of a sphere with the same volume as the tumor Fractional concavity 34 The ratio between the surface of the convex hull encompassing the tumor, and the actual surface of the tumor. = -

Fcc D 3
Surface of the convex hull A www.nature.com/scientificreports www.nature.com/scientificreports/ to a targeted edge length. The edge length was set to 2 and only one iteration was used. The last method (M4) used the Boundary function of MATLAB that returns a triangulation corresponding to a single conforming 3D boundary around the points. A shrinking factor of 1 was used to obtain the concave hull of the shape of interest.
Nine three-dimensional (3D) shape features were deduced from the surface and volume values extracted using the four meshing methods previously described. These included Volume, Surface Area, Surface-to-volume, three formulas for Compactness 8 , Spherical Disproportion, Sphericity 17 and Fractional Concavity 34 . Table 5 shows the description and formulas of the computed features. The surface of the convex hull included in the fractional concavity formula was obtained using the Boundary function and a shrinking value of 0.
To quantify the impact of 3D printing, acquisition, segmentation and meshing on the radiomic shape features, features were computed directly from the stereolithography (STL) file format that was used for 3D printing, and compared to those computed using the aforementioned Matlab functions.
To study the effect of resampling, masks extracted from the original images were resampled on a 1 × 1 × 1 mm 3 grid using a 3D linear interpolation. To remove the inherent variation on volume between objects, a homothetic transformation was then applied to bring back all the volumes to the value corresponding to the average of all volume values calculated for the 28 shapes for the M1 method. Finally, three representative phantoms (d = 11, 47, 92) were resampled to have 25, 50, 75, 125 and 150% of their original volume.
To validate the clinical relevance of the phantoms, CT-scans of lung tumors of three patients from the publicly available RIDER database 35 were contoured. Radiomic features were extracted and compared to the range of the values extracted from the 28 printed phantoms. statistical analysis. All statistical analyses were performed with R version 3.3.2 (https://www.r-project.org/). Differences between features computed directly from the STL files and those computed after meshing were compared. Percent changes between scans acquired with 0.6 and 2 mm slice thicknesses were quantified. To determine the relationship of each shape feature with tumor complexity, Spearman's rank-order correlation coefficients were computed for each of the four meshing methods. Complexity was considered as an ordinal variable with numeric values ranging from 11 (most spiculated) to 92 (least spiculated). Pairwise correlations among variables were also computed. To evaluate the effect of resampling on feature variation, percent changes were computed between features extracted from original and resampled (1 × 1 × 1 mm 3 grid) 3D masks for the M1 and M3 meshing methods. To evaluate the effect of changes in volume, the percent change in shape features when the volume varied from 25% to 150% was computed for d = 11, 47 and 92 and compared to the percent change observed when modifying the spiculatedness of the model from d = 11 to d = 92 for M3. Figure 7 summarizes the general schema of the methodology.

Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.