From 263e8b36b7a0b8bbeee9af8bdf4c98f333c2103d Mon Sep 17 00:00:00 2001 From: Joost van Griethuysen Date: Thu, 27 Sep 2018 14:59:35 +0200 Subject: [PATCH] STYL: Update shape documentation and feature names To better reflect computation method, explicitly name the two volume features `MeshVolume` and `VoxelVolume`. Add documentation specifying the unless otherwise mentioned, features are based on the mesh approximation of the shape. Only exceptions to this are the `VoxelVolume` and PCA-derived features (`MajorAxisLength`, `MinorAxisLength`,`LeastAxisLength`, `Elongation` and `Flatness`. For these features, it is explicitly stated in the documentation that these do not make use of the mesh approximation, but the segmented voxel(center)s themselves. --- data/baseline/baseline_shape.csv | 10 +- ...hape_params.yaml => example_allShape.yaml} | 33 ++-- radiomics/shape.py | 172 ++++++++++++------ 3 files changed, 134 insertions(+), 81 deletions(-) rename examples/exampleSettings/{Full_shape_params.yaml => example_allShape.yaml} (72%) diff --git a/data/baseline/baseline_shape.csv b/data/baseline/baseline_shape.csv index bb520b75..92904340 100644 --- a/data/baseline/baseline_shape.csv +++ b/data/baseline/baseline_shape.csv @@ -16,14 +16,14 @@ general_info_VoxelNum,4137,453,143,837,24644 original_shape_Maximum3DDiameter,69.6009903059,24.4445322363,9.08101860833,21.3449467202,68.9824710303 original_shape_Maximum2DDiameterSlice,53.5939777692,16.3129789202,8.39979371027,16.6273237626,55.8023479973 original_shape_Sphericity,0.480242186182,0.75130449808,0.695903713518,0.748037658194,0.672409930059 -original_shape_MinorAxis,34.849701666854706,13.383589687312323,5.4115477584202925,12.488622781409092,41.45975188669023 +original_shape_MinorAxisLength,34.849701666854706,13.383589687312323,5.4115477584202925,12.488622781409092,41.45975188669023 original_shape_Elongation,0.5621171627174109,0.7407691177548188,0.6999838102754056,0.7187910312752432,0.7433464635239126 original_shape_SurfaceVolumeRatio,0.398386088279,0.535570612898,1.39322476873,0.587720147764,0.197495275935 -original_shape_Volume,16149.4954427,1736.01786296,124.091304832,1330.97607931,48292.5216046 -original_shape_MajorAxis,61.99722046980879,18.067153943831187,7.730961315078334,17.374483317150464,55.77446577218634 +original_shape_MeshVolume,16149.4954427,1736.01786296,124.091304832,1330.97607931,48292.5216046 +original_shape_MajorAxisLength,61.99722046980879,18.067153943831187,7.730961315078334,17.374483317150464,55.77446577218634 original_shape_SurfaceArea,6433.7343171,929.760150865,172.887079476,782.241457999,9537.54487991 original_shape_Flatness,0.46105975346582634,0.6188162226844603,0.6842246761542068,0.5143357681770733,0.569289235382734 -original_shape_LeastAxis,28.584423185376522,11.180247958180265,5.289714502170174,8.93631822360633,31.75180297332843 +original_shape_LeastAxisLength,28.584423185376522,11.180247958180265,5.289714502170174,8.93631822360633,31.75180297332843 original_shape_Maximum2DDiameterColumn,49.4908549791,23.593253306,7.12140846294,19.5556659888,65.3030786496 original_shape_Maximum2DDiameterRow,65.8890595172,20.7289284953,8.24091097816,17.5574948384,62.1026289395 -original_shape_ApproximateVolume,16412.658691406243,1797.180175781249,132.42579545515324,1361.1978149414062,48434.108762463955 +original_shape_VoxelVolume,16412.658691406243,1797.180175781249,132.42579545515324,1361.1978149414062,48434.108762463955 diff --git a/examples/exampleSettings/Full_shape_params.yaml b/examples/exampleSettings/example_allShape.yaml similarity index 72% rename from examples/exampleSettings/Full_shape_params.yaml rename to examples/exampleSettings/example_allShape.yaml index ff03080d..e821abd9 100644 --- a/examples/exampleSettings/Full_shape_params.yaml +++ b/examples/exampleSettings/example_allShape.yaml @@ -14,19 +14,20 @@ setting: # for that class. Otherwise, the specified features are calculated, or, if none are specified, all are calculated (excluding redundant/deprecated features). featureClass: shape: # disable redundant Compactness 1 and Compactness 2 features by specifying all other shape features - - 'Volume' - - 'SurfaceArea' - - 'SurfaceVolumeRatio' - - 'Compactness1' - - 'Compactness2' - - 'Sphericity' - - 'SphericalDisproportion' - - 'Maximum3DDiameter' - - 'Maximum2DDiameterSlice' - - 'Maximum2DDiameterColumn' - - 'Maximum2DDiameterRow' - - 'MajorAxis' - - 'MinorAxis' - - 'LeastAxis' - - 'Elongation' - - 'Flatness' + - VoxelVolume + - MeshVolume + - SurfaceArea + - SurfaceVolumeRatio + - Compactness1 + - Compactness2 + - Sphericity + - SphericalDisproportion + - Maximum3DDiameter + - Maximum2DDiameterSlice + - Maximum2DDiameterColumn + - Maximum2DDiameterRow + - MajorAxisLength + - MinorAxisLength + - LeastAxisLength + - Elongation + - Flatness diff --git a/radiomics/shape.py b/radiomics/shape.py index fec52326..0e02fd81 100644 --- a/radiomics/shape.py +++ b/radiomics/shape.py @@ -10,10 +10,31 @@ class RadiomicsShape(base.RadiomicsFeaturesBase): are independent from the gray level intensity distribution in the ROI and are therefore only calculated on the non-derived image and mask. + Unless otherwise specified, features are derived from the approximated shape defined by the triangle mesh. To build + this mesh, vertices (points) are first defined as points halfway on an edge between a voxel included in the ROI and + one outside the ROI. By connecting these vertices a mesh of connected triangles is obtained, with each triangle + defined by 3 adjacent vertices, which shares each side with exactly one other triangle. + + This mesh is generated using a marching cubes algorithm. In this algorithm, a 2x2 cube is moved through the mask + space. For each position, the corners of the cube are then marked 'segmented' (1) or 'not segmented' (0). Treating the + corners as specific bits in a binary number, a unique cube-index is obtained (0-255). This index is then used to + determine which triangles are present in the cube, which are defined in a lookup table. + + These triangles are defined in such a way, that the normal (obtained from the cross product of vectors describing 2 + out of 3 edges) are always oriented in the same direction. For PyRadiomics, the calculated normals are always pointing + outward. This is necessary to obtain the correct signed volume used in calculation of ``MeshVolume``. + Let: - - :math:`V` the volume of the ROI in mm\ :sup:`3` - - :math:`A` the surface area of the ROI in mm\ :sup:`2` + - :math:`N_v` represent the number of voxels included in the ROI + - :math:`N_f` represent the number of faces (triangles) defining the Mesh. + - :math:`V` the volume of the mesh in mm\ :sup:`3`, calculated by :py:func:`getMeshVolumeFeatureValue` + - :math:`A` the surface area of the mesh in mm\ :sup:`2`, calculated by :py:func:`getMeshSurfaceAreaFeatureValue` + + References: + + - Lorensen WE, Cline HE. Marching cubes: A high resolution 3D surface construction algorithm. ACM SIGGRAPH Comput + Graph `Internet `_. 1987;21:163-9. """ def __init__(self, inputImage, inputMask, **kwargs): @@ -72,36 +93,42 @@ def _initSegmentBasedCalculation(self): self.logger.debug('Shape feature class initialized') - def getVolumeFeatureValue(self): + def getMeshVolumeFeatureValue(self): r""" - **1. Volume** + **1. Mesh Volume** .. math:: - V_{f} = \displaystyle\frac{a \dot (b \times c)}{6} (1) + V_i = \displaystyle\frac{Oa_i \cdot (Ob_i \times Oc_i)}{6} \text{ (1)} + + V = \displaystyle\sum^{N_f}_{i=1}{V_i} \text{ (2)} - V = \displaystyle\sum^{N_f}_{f=1}{V_f} (2) + The volume of the ROI :math:`V` is calculated from the triangle mesh of the ROI. + For each face :math:`i` in the mesh, defined by points :math:`a_i, b_i` and :math:`c_i`, the (signed) volume + :math:`V_f` of the tetrahedron defined by that face and the origin of the image (:math:`O`) is calculated. (1) + The sign of the volume is determined by the sign of the normal, which must be consistently defined as either facing + outward or inward of the ROI. - The volume of the ROI :math:`V` is calculated from the triangle mesh of the ROI. First the volume is calculated of - the tetrahedron that is defined by face f and the origin of the image (1) for all :math:`N_f` faces that make up the - triangle mesh. This mesh is also used for the calculation of surface area. By ensuring that the normals (obtained - from the cross product) have a consistent orientation (outward or inward facing), the volumes obtained are signed. - All sub-volumes are then summed (2), yielding the volume of the ROI. + Then taking the sum of all :math:`V_i`, the total volume of the ROI is obtained (2) .. note:: - For more extensive documentation on how the volume is obtained using the surface mesh, see the IBSI document on - the calculation of Volume. + For more extensive documentation on how the volume is obtained using the surface mesh, see the IBSI document, + where this feature is defined as ``Volume``. """ return self.Volume - def getApproximateVolumeFeatureValue(self): + def getVoxelVolumeFeatureValue(self): r""" - **1. Volume** + **2. Voxel Volume** .. math:: - V_{approx} = \displaystyle\sum^{N}_{i=1}{V_i} + V_{voxel} = \displaystyle\sum^{N_v}_{k=1}{V_k} - The volume of the ROI :math:`V` is approximated by multiplying the number of voxels in the ROI by the volume of a - single voxel :math:`V_i`. This is a less precise approximation of the volume and is not used in subsequent features. + The volume of the ROI :math:`V_{voxel}` is approximated by multiplying the number of voxels in the ROI by the volume + of a single voxel :math:`V_k`. This is a less precise approximation of the volume and is not used in subsequent + features. This feature does not make use of the mesh and is not used in calculation of other shape features. + + .. note:: + Defined in IBSI as ``Approximate Volume``. """ z, y, x = self.pixelSpacing Np = len(self.labelledVoxelCoordinates[0]) @@ -109,31 +136,29 @@ def getApproximateVolumeFeatureValue(self): def getSurfaceAreaFeatureValue(self): r""" - **2. Surface Area** + **3. Surface Area** .. math:: - A = \displaystyle\sum^{N}_{i=1}{\frac{1}{2}|\text{a}_i\text{b}_i \times \text{a}_i\text{c}_i|} + A_i = \frac{1}{2}|\text{a}_i\text{b}_i \times \text{a}_i\text{c}_i| \text{ (1)} - Where: + A = \displaystyle\sum^{N_f}_{i=1}{A_i} \text{ (2)} - :math:`N` is the number of triangles forming the surface mesh of the volume (ROI) + where: - :math:`\text{a}_i\text{b}_i` and :math:`\text{a}_i\text{c}_i` are the edges of the :math:`i^{\text{th}}` triangle - formed by points :math:`\text{a}_i`, :math:`\text{b}_i` and :math:`\text{c}_i` + :math:`\text{a}_i\text{b}_i` and :math:`\text{a}_i\text{c}_i` are edges of the :math:`i^{\text{th}}` triangle in the + mesh, formed by vertices :math:`\text{a}_i`, :math:`\text{b}_i` and :math:`\text{c}_i`. - Surface Area is an approximation of the surface of the ROI in mm2, calculated using a the triangle mesh generated - using the marching cubes algorithm. + To calculate the surface area, first the surface area :math:`A_i` of each triangle in the mesh is calculated (1). + The total surface area is then obtained by taking the sum of all calculated sub-areas (2). - References: - - - Lorensen WE, Cline HE. Marching cubes: A high resolution 3D surface construction algorithm. ACM SIGGRAPH Comput - Graph `Internet `_. 1987;21:163-9. + .. note:: + Defined in IBSI as ``Surface Area``. """ return self.SurfaceArea def getSurfaceVolumeRatioFeatureValue(self): r""" - **3. Surface Area to Volume ratio** + **4. Surface Area to Volume ratio** .. math:: \textit{surface to volume ratio} = \frac{A}{V} @@ -145,7 +170,7 @@ def getSurfaceVolumeRatioFeatureValue(self): def getSphericityFeatureValue(self): r""" - **4. Sphericity** + **5. Sphericity** .. math:: \textit{sphericity} = \frac{\sqrt[3]{36 \pi V^2}}{A} @@ -165,7 +190,7 @@ def getSphericityFeatureValue(self): @deprecated def getCompactness1FeatureValue(self): r""" - **5. Compactness 1** + **6. Compactness 1** .. math:: \textit{compactness 1} = \frac{V}{\sqrt{\pi A^3}} @@ -182,14 +207,15 @@ def getCompactness1FeatureValue(self): This feature is correlated to Compactness 2, Sphericity and Spherical Disproportion. Therefore, this feature is marked, so it is not enabled by default (i.e. this feature will not be enabled if no individual features are specified (enabling 'all' features), but will be enabled when individual features are - specified, including this feature). To include this feature in the extraction, specify it by name in the enabled features. + specified, including this feature). To include this feature in the extraction, specify it by name in the enabled + features. """ return self.Volume / (self.SurfaceArea ** (3.0 / 2.0) * numpy.sqrt(numpy.pi)) @deprecated def getCompactness2FeatureValue(self): r""" - **6. Compactness 2** + **7. Compactness 2** .. math:: \textit{compactness 2} = 36 \pi \frac{V^2}{A^3} @@ -204,14 +230,15 @@ def getCompactness2FeatureValue(self): This feature is correlated to Compactness 1, Sphericity and Spherical Disproportion. Therefore, this feature is marked, so it is not enabled by default (i.e. this feature will not be enabled if no individual features are specified (enabling 'all' features), but will be enabled when individual features are - specified, including this feature). To include this feature in the extraction, specify it by name in the enabled features. + specified, including this feature). To include this feature in the extraction, specify it by name in the enabled + features. """ return (36.0 * numpy.pi) * (self.Volume ** 2.0) / (self.SurfaceArea ** 3.0) @deprecated def getSphericalDisproportionFeatureValue(self): r""" - **7. Spherical Disproportion** + **8. Spherical Disproportion** .. math:: \textit{spherical disproportion} = \frac{A}{4\pi R^2} = \frac{A}{\sqrt[3]{36 \pi V^2}} @@ -227,13 +254,14 @@ def getSphericalDisproportionFeatureValue(self): This feature is correlated to Compactness 2, Compactness2 and Sphericity. Therefore, this feature is marked, so it is not enabled by default (i.e. this feature will not be enabled if no individual features are specified (enabling 'all' features), but will be enabled when individual features are - specified, including this feature). To include this feature in the extraction, specify it by name in the enabled features. + specified, including this feature). To include this feature in the extraction, specify it by name in the enabled + features. """ return self.SurfaceArea / (36 * numpy.pi * self.Volume ** 2) ** (1.0 / 3.0) def getMaximum3DDiameterFeatureValue(self): r""" - **8. Maximum 3D diameter** + **9. Maximum 3D diameter** Maximum 3D diameter is defined as the largest pairwise Euclidean distance between tumor surface mesh vertices. @@ -244,7 +272,7 @@ def getMaximum3DDiameterFeatureValue(self): def getMaximum2DDiameterSliceFeatureValue(self): r""" - **9. Maximum 2D diameter (Slice)** + **10. Maximum 2D diameter (Slice)** Maximum 2D diameter (Slice) is defined as the largest pairwise Euclidean distance between tumor surface mesh vertices in the row-column (generally the axial) plane. @@ -253,7 +281,7 @@ def getMaximum2DDiameterSliceFeatureValue(self): def getMaximum2DDiameterColumnFeatureValue(self): r""" - **10. Maximum 2D diameter (Column)** + **11. Maximum 2D diameter (Column)** Maximum 2D diameter (Column) is defined as the largest pairwise Euclidean distance between tumor surface mesh vertices in the row-slice (usually the coronal) plane. @@ -262,46 +290,61 @@ def getMaximum2DDiameterColumnFeatureValue(self): def getMaximum2DDiameterRowFeatureValue(self): r""" - **11. Maximum 2D diameter (Row)** + **12. Maximum 2D diameter (Row)** Maximum 2D diameter (Row) is defined as the largest pairwise Euclidean distance between tumor surface mesh vertices in the column-slice (usually the sagittal) plane. """ return self.diameters[2] - def getMajorAxisFeatureValue(self): + def getMajorAxisLengthFeatureValue(self): r""" - **12. Major Axis** + **13. Major Axis Length** .. math:: - \textit{major axis} = 4 \sqrt{\lambda_{\text{major}}} + \textit{major axis} = 4 \sqrt{\lambda_{major}} + This feature yield the largest axis length of the ROI-enclosing ellipsoid and is calculated using the largest + principal component :math:`\lambda_{major}`. + + The principal component analysis is performed using the physical coordinates of the voxel centers defining the ROI. + It therefore takes spacing into account, but does not make use of the shape mesh. """ if self.eigenValues[2] < 0: self.logger.warning('Major axis eigenvalue negative! (%g)', self.eigenValues[2]) return numpy.nan return numpy.sqrt(self.eigenValues[2]) * 4 - def getMinorAxisFeatureValue(self): + def getMinorAxisLengthFeatureValue(self): r""" - **13. Minor Axis** + **14. Minor Axis Length** .. math:: - \textit{minor axis} = 4 \sqrt{\lambda_{\text{minor}}} + \textit{minor axis} = 4 \sqrt{\lambda_{minor}} + + This feature yield the second-largest axis length of the ROI-enclosing ellipsoid and is calculated using the largest + principal component :math:`\lambda_{minor}`. + The principal component analysis is performed using the physical coordinates of the voxel centers defining the ROI. + It therefore takes spacing into account, but does not make use of the shape mesh. """ if self.eigenValues[1] < 0: self.logger.warning('Minor axis eigenvalue negative! (%g)', self.eigenValues[1]) return numpy.nan return numpy.sqrt(self.eigenValues[1]) * 4 - def getLeastAxisFeatureValue(self): + def getLeastAxisLengthFeatureValue(self): r""" - **14. Least Axis** + **15. Least Axis Length** .. math:: - \textit{least axis} = 4 \sqrt{\lambda_{\text{least}}} + \textit{least axis} = 4 \sqrt{\lambda_{least}} + This feature yield the smallest axis length of the ROI-enclosing ellipsoid and is calculated using the largest + principal component :math:`\lambda_{least}`. In case of a 2D segmentation, this value will be 0. + + The principal component analysis is performed using the physical coordinates of the voxel centers defining the ROI. + It therefore takes spacing into account, but does not make use of the shape mesh. """ if self.eigenValues[0] < 0: self.logger.warning('Least axis eigenvalue negative! (%g)', self.eigenValues[0]) @@ -310,17 +353,21 @@ def getLeastAxisFeatureValue(self): def getElongationFeatureValue(self): r""" - **15. Elongation** + **16. Elongation** - Elongation is calculated using its implementation in SimpleITK, and is defined as: + Elongation shows the relationship between the two largest principal components in the ROI shape. + For computational reasons, this feature is defined as the inverse of true elongation. .. math:: - \textit{elongation} = \sqrt{\frac{\lambda_{\text{minor}}}{\lambda_{\text{major}}}} + \textit{elongation} = \sqrt{\frac{\lambda_{minor}}{\lambda_{major}}} Here, :math:`\lambda_{\text{major}}` and :math:`\lambda_{\text{minor}}` are the lengths of the largest and second largest principal component axes. The values range between 1 (where the cross section through the first and second - largest principal moments is circle-like (non-elongated)) and 0 (where the object is a single point or 1 dimensional - line). + largest principal moments is circle-like (non-elongated)) and 0 (where the object is a maximally elongated: i.e. a 1 + dimensional line). + + The principal component analysis is performed using the physical coordinates of the voxel centers defining the ROI. + It therefore takes spacing into account, but does not make use of the shape mesh. """ if self.eigenValues[1] < 0 or self.eigenValues[2] < 0: self.logger.warning('Elongation eigenvalue negative! (%g, %g)', self.eigenValues[1], self.eigenValues[2]) @@ -329,15 +376,20 @@ def getElongationFeatureValue(self): def getFlatnessFeatureValue(self): r""" - **16. Flatness** + **17. Flatness** - Flatness is calculated using its implementation in SimpleITK, and is defined as: + Flatness shows the relationship between the largest and smallest principal components in the ROI shape. + For computational reasons, this feature is defined as the inverse of true flatness. .. math:: - \textit{flatness} = \sqrt{\frac{\lambda_{\text{least}}}{\lambda_{\text{major}}}} + \textit{flatness} = \sqrt{\frac{\lambda_{least}}{\lambda_{major}}} Here, :math:`\lambda_{\text{major}}` and :math:`\lambda_{\text{least}}` are the lengths of the largest and smallest - principal component axes. The values range between 1 (non-flat, sphere-like) and 0 (a flat object). + principal component axes. The values range between 1 (non-flat, sphere-like) and 0 (a flat object, or single-slice + segmentation). + + The principal component analysis is performed using the physical coordinates of the voxel centers defining the ROI. + It therefore takes spacing into account, but does not make use of the shape mesh. """ if self.eigenValues[0] < 0 or self.eigenValues[2] < 0: self.logger.warning('Elongation eigenvalue negative! (%g, %g)', self.eigenValues[0], self.eigenValues[2])