ENH: Define return values for 0-division errors

For some features the formula would return NaN in specific cases (e.g. a flat region). Document and return predefined values in these cases.

radiomics/firstorder.py

@@ -62,7 +62,7 @@ def getEnergyFeatureValue(self):
 
     shiftedParameterArray = self.targetVoxelArray + self.voxelArrayShift
 
-    return (numpy.sum(shiftedParameterArray ** 2))
+    return numpy.sum(shiftedParameterArray ** 2)
 
   def getTotalEnergyFeatureValue(self):
     r"""
@@ -81,7 +81,7 @@ def getTotalEnergyFeatureValue(self):
     x, y, z = self.pixelSpacing
     cubicMMPerVoxel = x * y * z
 
-    return (cubicMMPerVoxel * self.getEnergyFeatureValue())
+    return cubicMMPerVoxel * self.getEnergyFeatureValue()
 
   def getEntropyFeatureValue(self):
     r"""
@@ -98,12 +98,14 @@ def getEntropyFeatureValue(self):
     eps = numpy.spacing(1)
     binEdges = imageoperations.getBinEdges(self.binWidth, self.targetVoxelArray)
     bins = numpy.histogram(self.targetVoxelArray, binEdges)[0]
-    try:
-      bins = bins + eps
-      bins = bins / float(bins.sum())
-      return (-1.0 * numpy.sum(bins * numpy.log2(bins)))
-    except ZeroDivisionError:
-      return numpy.core.nan
+
+    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))
 
   def getMinimumFeatureValue(self):
     r"""
@@ -112,21 +114,21 @@ def getMinimumFeatureValue(self):
     :math:`minimum = \min(\textbf{X})`
     """
 
-    return (numpy.min(self.targetVoxelArray))
+    return numpy.min(self.targetVoxelArray)
 
   def get10PercentileFeatureValue(self):
     r"""
     Calculate the 10\ :sup:`th` percentile in the image array.
     """
 
-    return (numpy.percentile(self.targetVoxelArray, 10))
+    return numpy.percentile(self.targetVoxelArray, 10)
 
   def get90PercentileFeatureValue(self):
     r"""
     Calculate the 90\ :sup:`th` percentile in the image array.
     """
 
-    return (numpy.percentile(self.targetVoxelArray, 90))
+    return numpy.percentile(self.targetVoxelArray, 90)
 
   def getMaximumFeatureValue(self):
     r"""
@@ -135,7 +137,7 @@ def getMaximumFeatureValue(self):
     :math:`maximum = \max(\textbf{X})`
     """
 
-    return (numpy.max(self.targetVoxelArray))
+    return numpy.max(self.targetVoxelArray)
 
   def getMeanFeatureValue(self):
     r"""
@@ -144,14 +146,14 @@ def getMeanFeatureValue(self):
     :math:`mean = \frac{1}{N}\displaystyle\sum^{N}_{i=1}{\textbf{X}(i)}`
     """
 
-    return (numpy.mean(self.targetVoxelArray))
+    return numpy.mean(self.targetVoxelArray)
 
   def getMedianFeatureValue(self):
     r"""
     Calculate the Median Value for the image array.
     """
 
-    return (numpy.median(self.targetVoxelArray))
+    return numpy.median(self.targetVoxelArray)
 
   def getInterquartileRangeFeatureValue(self):
     r"""
@@ -170,7 +172,7 @@ def getRangeFeatureValue(self):
     :math:`range = \max(\textbf{X}) - \min(\textbf{X})`
     """
 
-    return (numpy.max(self.targetVoxelArray) - numpy.min(self.targetVoxelArray))
+    return numpy.max(self.targetVoxelArray) - numpy.min(self.targetVoxelArray)
 
   def getMeanAbsoluteDeviationFeatureValue(self):
     r"""
@@ -181,7 +183,7 @@ def getMeanAbsoluteDeviationFeatureValue(self):
     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"""
@@ -215,8 +217,12 @@ def getRootMeanSquaredFeatureValue(self):
     volume-confounding.
     """
 
+    # 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)))
+    return numpy.sqrt((numpy.sum(shiftedParameterArray ** 2)) / float(shiftedParameterArray.size))
 
   def getStandardDeviationFeatureValue(self):
     r"""
@@ -227,7 +233,7 @@ def getStandardDeviationFeatureValue(self):
     Standard Deviation measures the amount of variation or dispersion from the Mean Value.
     """
 
-    return (numpy.std(self.targetVoxelArray))
+    return numpy.std(self.targetVoxelArray)
 
   def getSkewnessFeatureValue(self, axis=0):
     r"""
@@ -245,12 +251,18 @@ def getSkewnessFeatureValue(self, axis=0):
     Related links:
 
     https://en.wikipedia.org/wiki/Skewness
+
+    .. note::
+
+      In case of a flat region, the standard deviation and 4\ :sup:`rd` central moment will be both 0. In this case, a
+      value of 0 is returned.
     """
 
     m2 = self._moment(self.targetVoxelArray, 2, axis)
     m3 = self._moment(self.targetVoxelArray, 3, axis)
 
-    if (m2 == 0): return numpy.core.nan
+    if m2 == 0:  # Flat Region
+      return 0
 
     return m3 / m2 ** 1.5
 
@@ -271,12 +283,18 @@ def getKurtosisFeatureValue(self, axis=0):
     Related links:
 
     https://en.wikipedia.org/wiki/Kurtosis
+
+    .. note::
+
+      In case of a flat region, the standard deviation and 4\ :sup:`rd` central moment will be both 0. In this case, a
+      value of 0 is returned.
     """
 
     m2 = self._moment(self.targetVoxelArray, 2, axis)
     m4 = self._moment(self.targetVoxelArray, 4, axis)
 
-    if (m2 == 0): return numpy.core.nan
+    if m2 == 0:  # Flat Region
+      return 0
 
     return m4 / m2 ** 2.0
 
@@ -290,7 +308,7 @@ def getVarianceFeatureValue(self):
     the spread of the distribution about the mean.
     """
 
-    return (numpy.std(self.targetVoxelArray) ** 2)
+    return numpy.std(self.targetVoxelArray) ** 2
 
   def getUniformityFeatureValue(self):
     r"""
@@ -306,8 +324,10 @@ def getUniformityFeatureValue(self):
     eps = numpy.spacing(1)
     binEdges = imageoperations.getBinEdges(self.binWidth, self.targetVoxelArray)
     bins = numpy.histogram(self.targetVoxelArray, binEdges)[0]
-    try:
-      bins = bins / (float(bins.sum() + eps))
-      return (numpy.sum(bins ** 2))
-    except ZeroDivisionError:
-      return numpy.core.nan
+    sumBins = bins.sum()
+
+    if sumBins == 0:  # No segmented voxels
+      return 0
+
+    bins = bins / (float(sumBins + eps))
+    return numpy.sum(bins ** 2)