myshape.py 4.01 KB
 Alessia Marcolini committed Mar 16, 2020 1 2 3 4 5 6 7 8 9 10 11 import numpy from scipy.spatial import ConvexHull from radiomics.shape import RadiomicsShape from .mybase import MyRadiomicsFeaturesBase class MyRadiomicsShape(MyRadiomicsFeaturesBase): def __init__(self, inputImage, inputMask, **kwargs): super().__init__(inputImage, inputMask, **kwargs) def _initVoxelBasedCalculation(self):  Alessia Marcolini committed Mar 16, 2020 12  return RadiomicsShape._initVoxelBasedCalculation(self)  Alessia Marcolini committed Mar 16, 2020 13 14  def _initSegmentBasedCalculation(self):  Alessia Marcolini committed Mar 16, 2020 15  return RadiomicsShape._initSegmentBasedCalculation(self)  Alessia Marcolini committed Mar 16, 2020 16 17  def _calculateSurfaceArea(self):  Alessia Marcolini committed Mar 16, 2020 18  return RadiomicsShape._calculateSurfaceArea(self)  Alessia Marcolini committed Mar 16, 2020 19 20  def calculateDiameters(self):  Alessia Marcolini committed Mar 16, 2020 21  return RadiomicsShape.calculateDiameters(self)  Alessia Marcolini committed Mar 16, 2020 22 23  def _interpolate(self, grid, p1, p2):  Alessia Marcolini committed Mar 16, 2020 24  return RadiomicsShape._interpolate(self, grid, p1, p2)  Alessia Marcolini committed Mar 16, 2020 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45  def getCompactness2FeatureValue(self): r""" **6. Compactness 2** .. math:: \textit{compactness 2} = 36 \pi \frac{V^2}{A^3} Similar to Sphericity and Compactness 1, Compactness 2 is a measure of how compact the shape of the tumor is relative to a sphere (most compact). It is a dimensionless measure, independent of scale and orientation. The value range is :math:0 < compactness\ 2 \leq 1, where a value of 1 indicates a perfect sphere. By definition, :math:compactness\ 2 = (sphericity)^3 .. note:: 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. """ if self.SurfaceArea is None:  Alessia Marcolini committed Mar 16, 2020 46  self.SurfaceArea = self._calculateSurfaceArea()  Alessia Marcolini committed Mar 16, 2020 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105  return (36.0 * numpy.pi) * (self.Volume ** 2.0) / (self.SurfaceArea ** 3.0) def getEccentricityFeatureValue(self): r""" **17. Eccentricity** This is a custom shape feature, implemented for comparison with Valieres et al. Eccentricity is defined as: .. math:: \textit{Eccentricity} = \sqrt{1-\frac{a\times b}{c^2}} Here, where c is the longest semi-principal axes of the ellipsoid, and a and b are the second and third longest semi-principal axes of the ellipsoid. """ major_axis = numpy.sqrt(self.eigenValues[2]) * 4 minor_axis = numpy.sqrt(self.eigenValues[1]) * 4 least_axis = numpy.sqrt(self.eigenValues[0]) * 4 return numpy.sqrt(1.0 - (minor_axis * least_axis / (major_axis ** 2))) def getSolidityFeatureValue(self): r""" **18. Solidity** This is a custom shape feature, implemented for comparison with Valieres et al. Eccentricity is defined as: Solidity is defined as: Ratio of the number of voxels in the tumour region to the number of voxels in the 3D convex hull of the tumour region (smallest polyhedron containing the tumour region) """ Np = len(self.labelledVoxelCoordinates[0]) coordinates = numpy.array(self.labelledVoxelCoordinates, dtype='int').transpose( (1, 0) ) # Transpose equals zip(*a) physicalCoordinates = coordinates * self.pixelSpacing[None, :] physicalCoordinates -= numpy.mean(physicalCoordinates, axis=0) # Centered at 0 physicalCoordinates /= numpy.sqrt(Np) ch = ConvexHull(physicalCoordinates) simplices = numpy.column_stack( (numpy.repeat(ch.vertices[0], ch.nsimplex), ch.simplices) ) tets = ch.points[simplices] a, b, c, d = tets[:, 0], tets[:, 1], tets[:, 2], tets[:, 3] tetrahedron_vol = ( numpy.abs(numpy.einsum('ij,ij->i', a - d, numpy.cross(b - d, c - d))) / 6 ) ROIVolume = self.Volume return ROIVolume / numpy.sum(tetrahedron_vol)