Commit b3cd2002 authored by Nicole Bussola's avatar Nicole Bussola
Browse files

adding SUVMax and SUVMean features for PET modality

parent 5ed18200
import numpy
from scipy import ndimage
from radiomics import base, deprecated, imageoperations
from oct2py import Oct2Py, Oct2PyError
class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
r"""
r"""
First-order statistics describe the distribution of voxel intensities within the image region defined by the mask
through commonly used and basic metrics.
......@@ -28,44 +30,50 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
in PyRadiomics, set ``voxelArrayShift`` to 0.
"""
def __init__(self, inputImage, inputMask, **kwargs):
super(RadiomicsFirstOrder, self).__init__(inputImage, inputMask, **kwargs)
def __init__(self, inputImage, inputMask, **kwargs):
super(RadiomicsFirstOrder, self).__init__(inputImage, inputMask, **kwargs)
self.pixelSpacing = inputImage.GetSpacing()
self.voxelArrayShift = kwargs.get('voxelArrayShift', 0)
self.pixelSpacing = inputImage.GetSpacing()
self.voxelArrayShift = kwargs.get('voxelArrayShift', 0)
def _initCalculation(self):
super(RadiomicsFirstOrder, self)._initSegmentBasedCalculation()
self.targetVoxelArray = self.imageArray[self.labelledVoxelCoordinates].astype('float')
self.discretizedTargetVoxelArray = None # Lazy instantiation
def _initCalculation(self):
super(RadiomicsFirstOrder, self)._initSegmentBasedCalculation()
self.targetVoxelArray = self.imageArray[self.labelledVoxelCoordinates].astype(
'float'
)
self.discretizedTargetVoxelArray = None # Lazy instantiation
self.logger.debug('First order feature class initialized')
self.logger.debug('First order feature class initialized')
def _moment(self, a, moment=1, axis=0):
r"""
def _moment(self, a, moment=1, axis=0):
r"""
Calculate n-order moment of an array for a given axis
"""
if moment == 1:
return numpy.float64(0.0)
else:
mn = numpy.mean(a, axis, keepdims=True)
s = numpy.power((a - mn), moment)
return numpy.mean(s, axis)
def _getDiscretizedTargetVoxelArray(self):
if self.discretizedTargetVoxelArray is None:
if self.binCount is not None:
binEdges = self.binCount
else:
binEdges = imageoperations.getBinEdges(self.binWidth, self.targetVoxelArray)
self.discretizedTargetVoxelArray = numpy.histogram(self.targetVoxelArray, binEdges)[0]
return self.discretizedTargetVoxelArray
def getEnergyFeatureValue(self):
r"""
if moment == 1:
return numpy.float64(0.0)
else:
mn = numpy.mean(a, axis, keepdims=True)
s = numpy.power((a - mn), moment)
return numpy.mean(s, axis)
def _getDiscretizedTargetVoxelArray(self):
if self.discretizedTargetVoxelArray is None:
if self.binCount is not None:
binEdges = self.binCount
else:
binEdges = imageoperations.getBinEdges(
self.binWidth, self.targetVoxelArray
)
self.discretizedTargetVoxelArray = numpy.histogram(
self.targetVoxelArray, binEdges
)[0]
return self.discretizedTargetVoxelArray
def getEnergyFeatureValue(self):
r"""
**1. Energy**
.. math::
......@@ -82,12 +90,12 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
This feature is volume-confounded, a larger value of :math:`c` increases the effect of volume-confounding.
"""
shiftedParameterArray = self.targetVoxelArray + self.voxelArrayShift
shiftedParameterArray = self.targetVoxelArray + self.voxelArrayShift
return numpy.sum(shiftedParameterArray ** 2)
return numpy.sum(shiftedParameterArray ** 2)
def getTotalEnergyFeatureValue(self):
r"""
def getTotalEnergyFeatureValue(self):
r"""
**2. Total Energy**
.. math::
......@@ -106,13 +114,13 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
Not present in IBSI feature definitions
"""
x, y, z = self.pixelSpacing
cubicMMPerVoxel = x * y * z
x, y, z = self.pixelSpacing
cubicMMPerVoxel = x * y * z
return cubicMMPerVoxel * self.getEnergyFeatureValue()
return cubicMMPerVoxel * self.getEnergyFeatureValue()
def getEntropyFeatureValue(self):
r"""
def getEntropyFeatureValue(self):
r"""
**3. Entropy**
.. math::
......@@ -127,47 +135,47 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
Defined by IBSI as Intensity Histogram Entropy.
"""
eps = numpy.spacing(1)
bins = self._getDiscretizedTargetVoxelArray()
eps = numpy.spacing(1)
bins = self._getDiscretizedTargetVoxelArray()
sumBins = bins.sum()
if sumBins == 0: # No segmented voxels
return 0
sumBins = bins.sum()
if sumBins == 0: # No segmented voxels
return 0
bins = bins + eps
bins = bins / float(sumBins)
return -1.0 * numpy.sum(bins * numpy.log2(bins))
bins = bins + eps
bins = bins / float(sumBins)
return -1.0 * numpy.sum(bins * numpy.log2(bins))
def getMinimumFeatureValue(self):
r"""
def getMinimumFeatureValue(self):
r"""
**4. Minimum**
.. math::
\textit{minimum} = \min(\textbf{X})
"""
return numpy.min(self.targetVoxelArray)
return numpy.min(self.targetVoxelArray)
def get10PercentileFeatureValue(self):
r"""
def get10PercentileFeatureValue(self):
r"""
**5. 10th percentile**
The 10\ :sup:`th` percentile of :math:`\textbf{X}`
"""
return numpy.percentile(self.targetVoxelArray, 10)
return numpy.percentile(self.targetVoxelArray, 10)
def get90PercentileFeatureValue(self):
r"""
def get90PercentileFeatureValue(self):
r"""
**6. 90th percentile**
The 90\ :sup:`th` percentile of :math:`\textbf{X}`
"""
return numpy.percentile(self.targetVoxelArray, 90)
return numpy.percentile(self.targetVoxelArray, 90)
def getMaximumFeatureValue(self):
r"""
def getMaximumFeatureValue(self):
r"""
**7. Maximum**
.. math::
......@@ -176,10 +184,10 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
The maximum gray level intensity within the ROI.
"""
return numpy.max(self.targetVoxelArray)
return numpy.max(self.targetVoxelArray)
def getMeanFeatureValue(self):
r"""
def getMeanFeatureValue(self):
r"""
**8. Mean**
.. math::
......@@ -188,19 +196,19 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
The average gray level intensity within the ROI.
"""
return numpy.mean(self.targetVoxelArray)
return numpy.mean(self.targetVoxelArray)
def getMedianFeatureValue(self):
r"""
def getMedianFeatureValue(self):
r"""
**9. Median**
The median gray level intensity within the ROI.
"""
return numpy.median(self.targetVoxelArray)
return numpy.median(self.targetVoxelArray)
def getInterquartileRangeFeatureValue(self):
r"""
def getInterquartileRangeFeatureValue(self):
r"""
**10. Interquartile Range**
.. math::
......@@ -210,10 +218,12 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
image array, respectively.
"""
return numpy.percentile(self.targetVoxelArray, 75) - numpy.percentile(self.targetVoxelArray, 25)
return numpy.percentile(self.targetVoxelArray, 75) - numpy.percentile(
self.targetVoxelArray, 25
)
def getRangeFeatureValue(self):
r"""
def getRangeFeatureValue(self):
r"""
**11. Range**
.. math::
......@@ -222,10 +232,10 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
The range of gray values in the ROI.
"""
return numpy.max(self.targetVoxelArray) - numpy.min(self.targetVoxelArray)
return numpy.max(self.targetVoxelArray) - numpy.min(self.targetVoxelArray)
def getMeanAbsoluteDeviationFeatureValue(self):
r"""
def getMeanAbsoluteDeviationFeatureValue(self):
r"""
**12. Mean Absolute Deviation (MAD)**
.. math::
......@@ -234,10 +244,12 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
Mean Absolute Deviation is the mean distance of all intensity values from the Mean Value of the image array.
"""
return numpy.mean(numpy.absolute((numpy.mean(self.targetVoxelArray) - self.targetVoxelArray)))
return numpy.mean(
numpy.absolute((numpy.mean(self.targetVoxelArray) - self.targetVoxelArray))
)
def getRobustMeanAbsoluteDeviationFeatureValue(self):
r"""
def getRobustMeanAbsoluteDeviationFeatureValue(self):
r"""
**13. Robust Mean Absolute Deviation (rMAD)**
.. math::
......@@ -249,14 +261,16 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
to the 10\ :sup:`th` and 90\ :sup:`th` percentile.
"""
prcnt10 = self.get10PercentileFeatureValue()
prcnt90 = self.get90PercentileFeatureValue()
percentileArray = self.targetVoxelArray[(self.targetVoxelArray >= prcnt10) * (self.targetVoxelArray <= prcnt90)]
prcnt10 = self.get10PercentileFeatureValue()
prcnt90 = self.get90PercentileFeatureValue()
percentileArray = self.targetVoxelArray[
(self.targetVoxelArray >= prcnt10) * (self.targetVoxelArray <= prcnt90)
]
return numpy.mean(numpy.absolute(percentileArray - numpy.mean(percentileArray)))
return numpy.mean(numpy.absolute(percentileArray - numpy.mean(percentileArray)))
def getRootMeanSquaredFeatureValue(self):
r"""
def getRootMeanSquaredFeatureValue(self):
r"""
**14. Root Mean Squared (RMS)**
.. math::
......@@ -271,16 +285,18 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
volume-confounding.
"""
# If no voxels are segmented, prevent division by 0 and return 0
if self.targetVoxelArray.size == 0:
return 0
# If no voxels are segmented, prevent division by 0 and return 0
if self.targetVoxelArray.size == 0:
return 0
shiftedParameterArray = self.targetVoxelArray + self.voxelArrayShift
return numpy.sqrt((numpy.sum(shiftedParameterArray ** 2)) / float(shiftedParameterArray.size))
shiftedParameterArray = self.targetVoxelArray + self.voxelArrayShift
return numpy.sqrt(
(numpy.sum(shiftedParameterArray ** 2)) / float(shiftedParameterArray.size)
)
@deprecated
def getStandardDeviationFeatureValue(self):
r"""
@deprecated
def getStandardDeviationFeatureValue(self):
r"""
**15. Standard Deviation**
.. math::
......@@ -297,10 +313,10 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
Not present in IBSI feature definitions (correlated with variance)
"""
return numpy.std(self.targetVoxelArray)
return numpy.std(self.targetVoxelArray)
def getSkewnessFeatureValue(self, axis=0):
r"""
def getSkewnessFeatureValue(self, axis=0):
r"""
**16. Skewness**
.. math::
......@@ -322,16 +338,16 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
value of 0 is returned.
"""
m2 = self._moment(self.targetVoxelArray, 2, axis)
m3 = self._moment(self.targetVoxelArray, 3, axis)
m2 = self._moment(self.targetVoxelArray, 2, axis)
m3 = self._moment(self.targetVoxelArray, 3, axis)
if m2 == 0: # Flat Region
return 0
if m2 == 0: # Flat Region
return 0
return m3 / m2 ** 1.5
return m3 / m2 ** 1.5
def getKurtosisFeatureValue(self, axis=0):
r"""
def getKurtosisFeatureValue(self, axis=0):
r"""
**17. Kurtosis**
.. math::
......@@ -358,16 +374,16 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
distributions. The PyRadiomics kurtosis is not corrected, yielding a value 3 higher than the IBSI kurtosis.
"""
m2 = self._moment(self.targetVoxelArray, 2, axis)
m4 = self._moment(self.targetVoxelArray, 4, axis)
m2 = self._moment(self.targetVoxelArray, 2, axis)
m4 = self._moment(self.targetVoxelArray, 4, axis)
if m2 == 0: # Flat Region
return 0
if m2 == 0: # Flat Region
return 0
return m4 / m2 ** 2.0
return m4 / m2 ** 2.0
def getVarianceFeatureValue(self):
r"""
def getVarianceFeatureValue(self):
r"""
**18. Variance**
.. math::
......@@ -377,10 +393,10 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
the spread of the distribution about the mean. By definition, :math:`\textit{variance} = \sigma^2`
"""
return numpy.std(self.targetVoxelArray) ** 2
return numpy.std(self.targetVoxelArray) ** 2
def getUniformityFeatureValue(self):
r"""
def getUniformityFeatureValue(self):
r"""
**19. Uniformity**
.. math::
......@@ -394,12 +410,203 @@ class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):
Defined by IBSI as Intensity Histogram Uniformity.
"""
eps = numpy.spacing(1)
bins = self._getDiscretizedTargetVoxelArray()
sumBins = bins.sum()
eps = numpy.spacing(1)
bins = self._getDiscretizedTargetVoxelArray()
sumBins = bins.sum()
if sumBins == 0: # No segmented voxels
return 0
bins = bins / (float(sumBins + eps))
return numpy.sum(bins ** 2)
def getSuvMaxFeatureValue(self):
r"""
**20. SUV max**
Implemented as one of the intensity features extracted by valieres et al.
Maximum SUV of the tumour region.
- input: 3D array representing the PET volume in SUV format, with
voxels outside the ROI set to NaNs.
.. note::
Extracted from PET scans and not used in the CT feature set.
"""
ROIPet = self.imageArray
mask = self.maskArray
ROIPet[~mask] = numpy.nan
return numpy.max(ROIPet[~numpy.isnan(ROIPet)])
def getSuvMeanFeatureValue(self):
r"""
**21. SUV mean**
Implemented as one of the intensity features extracted by valieres et al.
Average SUV of the tumour region.
- input: 3D array representing the PET volume in SUV format, with
voxels outside the ROI set to NaNs.
.. note::
Extracted from PET scans and not used in the CT feature set.
"""
ROIPet = self.imageArray
mask = self.maskArray
ROIPet[~mask] = numpy.nan
return numpy.mean(ROIPet[~numpy.isnan(ROIPet)])
def getTLGFeatureValue(self):
r"""
**21. TLG**
Total lesion glycolysis.
Defined as
.. math::
\textit{SUVMean}\times \textit{total volume of the tumour region}
- input: 3D array representing the PET volume in SUV format, with
voxels outside the ROI set to NaNs.
.. note::
Extracted from PET scans and not used in the CT feature set.
"""
z, y, x = self.pixelSpacing
Np = len(self.labelledVoxelCoordinates[0])
volume = Np * (z * x * y)
if sumBins == 0: # No segmented voxels
return 0
return numpy.mean(self.targetVoxelArray) * volume
bins = bins / (float(sumBins + eps))
return numpy.sum(bins ** 2)
def getInactiveVolumeFeatureValue(self):
r"""
**21. Inactive volume **
Percentage of the tumour region that is inactive. Ast defined by Valieres et al
a threshold of 0.01 × (SUVmax)^2 followed by closing and opening morphological operations were used
to differentiate active and inactive regions on FDG-PET scans
- input: 3D array representing the PET volume in SUV format, with
voxels outside the ROI set to NaNs.
.. note::
Extracted from PET scans and not used in the CT feature set.
"""
ROIPet = self.imageArray
mask = self.maskArray
ROIPet[~mask] = 0
ROIPet = numpy.pad(a, 1, mode='constant', constant_values=np.nan)
thresh = 0.01 * (numpy.max(ROIPet[~numpy.isnan(ROIPet)]) ** 2)
# mask = ROIPet[~numpy.isnan(ROIPet)]
# ROIPet[numpy.isnan(ROIPet)] = 0
mask_inactive = ROIPet > thresh
# MORPHOLOGICAL OPERATIONS
conn = numpy.zeros([5, 5, 5])
conn[:, :, 0] = numpy.array(
[
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
]
)
conn[:, :, 4] = conn[:, :, 0]
conn[:, :, 1] = numpy.array(
[
[0, 0, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 1, 1, 1, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 0, 0],
]
)
conn[:, :, 3] = conn[:, :, 1]
conn[:, :, 2] = numpy.array(
[
[0, 0, 1, 0, 0],
[0, 1, 1, 1, 0],
[1, 1, 1, 1, 1],
[0, 1, 1, 1, 0],
[0, 0, 1, 0, 0],
]
)
with Oct2Py() as oc:
try:
oc.eval('pkg load image')
perimeter = oc.bwperim(mask, 26) # oc.bwperim(mask, 26)
except Oct2PyError as e:
print(e)
oc.exit()
# mask_inactive = oc.imclose(mask_inactive, conn)
# mask_inactive = oc.imopen(mask_inactive, conn)
import ipdb
ipdb.set_trace()
new_mask = mask_inactive + perimeter
return n
def getAUCCSHFeatureValue(self):
r"""
**21. Inactive volume **
Area under the curve of the cumulative SUV-volume histogram describing the percentage
of total tumour volume above a percentage threshold of maximum SUV, as defined by van Velden et al.
- input: 3D array representing the PET volume in SUV format, with
voxels outside the ROI set to NaNs.
.. note::
Extracted from PET scans and not used in the CT feature set.
"""
nBins = 1000 # By default.
ROIPet = self.imageArray
mask = self.maskArray
ROIPet[~mask] = numpy.nan
outliers = numpy.where(
ROIPet
> (
numpy.mean(ROIPet[~numpy.isnan(ROIPet)])
+ 3 * numpy.std(ROIPet[~numpy.isnan(ROIPet)])
)
)[0]
goodVoxels = numpy.where(
ROIPet
<= (
numpy.mean(ROIPet[~numpy.isnan(ROIPet)])
+ 3 * numpy.std(ROIPet[~numpy.isnan(ROIPet)])
)
)[0]
ROIPet[outliers] = numpy.mean(ROIPet[goodVoxels])
ROIPet = ROIPet - numpy.min(ROIPet)
ROIPet = ROIPet / numpy.max(ROIPet)
volume = ROIPet[~numpy.isnan(ROIPet)]
nVoxel = len(volume)
bins = numpy.histogram(volume, nBins)
import ipdb
ipdb.set_trace()
csh = numpy.fliplr(numpy.cumsum(numpy.fliplr(bins)) / nVoxel)
return numpy.sum(csh / nBins)
# nVoxel = numel(volume);
# bins = hist(volume,nBins);
# CSH = fliplr(cumsum(fliplr(bins))./nVoxel);
# aucCSH = sum(CSH./nBins);
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment