Merge pull request #209 from JoostJM/add-normalize-function

Implement normalization function in imageoperations

Author: JoostJM

data/paramSchema.yaml

@@ -18,6 +18,16 @@ mapping:
         type: float
         range:
           min-ex: 0
+      normalize:
+        type: bool
+      normalizeScale:
+        type: float
+        range:
+          min-ex: 0
+      removeOutliers:
+        type: float
+        range:
+          min-ex: 0
       resampledPixelSpacing:
         seq:
           - type: int

radiomics/featureextractor.py

@@ -15,7 +15,7 @@
 
 
 class RadiomicsFeaturesExtractor:
-  """
+  r"""
   Wrapper class for calculation of a radiomics signature.
   At and after initialisation various settings can be used to customize the resultant signature.
   This includes which classes and features to use, as well as what should be done in terms of preprocessing the image
@@ -28,14 +28,24 @@ class RadiomicsFeaturesExtractor:
   It initialisation, a parameters file can be provided containing all necessary settings. This is done by passing the
   location of the file as the single argument in the initialization call, without specifying it as a keyword argument.
   If such a file location is provided, any additional kwargs are ignored.
-  Alternatively, at initialisation, the following general settings can be specified in kwargs:
+  Alternatively, at initialisation, the following general settings can be specified in the parameter file or ``kwargs``,
+  with default values in brackets:
 
   - verbose [False]: Boolean, set to False to disable status update printing.
   - enableCExtensions [True]: Boolean, set to False to force calculation to full-python mode. See also
     :py:func:`~radiomics.enableCExtensions()`.
   - additionalInfo [True]: boolean, set to False to disable inclusion of additional information on the extraction in the
     output. See also :py:func:`~addProvenance()`.
   - binWidth [25]: Float, size of the bins when making a histogram and for discretization of the image gray level.
+  - normalize [False]: Boolean, set to True to enable normalizing of the image before any resampling. See also
+    :py:func:`~imageoperations.normalizeImage`.
+  - normalizeScale [1]: Float, determines the scale after normalizing the image. If normalizing is disabled, this has
+    no effect.
+  - removeOutliers [None]: Float, defines the outliers to remove from the image. An outlier is defined as values that
+    differ more than :math:`n\sigma_x` from the mean, where :math:`n>0` and equal to the value of this setting. If this
+    parameter is omitted (providing it without a value (i.e. None) in the parameter file will throw an error), no
+    outliers are removed. If normalizing is disabled, this has no effect. See also
+    :py:func:`~imageoperations.normalizeImage`.
   - resampledPixelSpacing [None]: List of 3 floats, sets the size of the voxel in (x, y, z) plane when resampling.
   - interpolator [sitkBSpline]: Simple ITK constant or string name thereof, sets interpolator to use for resampling.
     Enumerated value, possible values:
@@ -82,7 +92,10 @@ def __init__(self, *args, **kwargs):
       self.loadParams(args[0])
     else:
       # Set default settings and update with and changed settings contained in kwargs
-      self.kwargs = {'resampledPixelSpacing': None,  # No resampling by default
+      self.kwargs = {'normalize': False,
+                     'normalizeScale': 1,
+                     'removeOutliers': None,
+                     'resampledPixelSpacing': None,  # No resampling by default
                      'interpolator': sitk.sitkBSpline,
                      'padDistance': 5,
                      'label': 1,
@@ -171,7 +184,10 @@ def loadParams(self, paramsFile):
       self.enabledFeatures = enabledFeatures
 
     # Set default settings and update with and changed settings contained in kwargs
-    self.kwargs = {'resampledPixelSpacing': None,  # No resampling by default
+    self.kwargs = {'normalize': False,
+                   'normalizeScale': 1,
+                   'removeOutliers': None,
+                   'resampledPixelSpacing': None,  # No resampling by default
                    'interpolator': sitk.sitkBSpline,
                    'padDistance': 5,
                    'label': 1,
@@ -287,11 +303,11 @@ def enableFeaturesByName(self, **enabledFeatures):
   def execute(self, imageFilepath, maskFilepath, label=None):
     """
     Compute radiomics signature for provide image and mask combination.
-    First, image and mask are loaded and resampled if necessary. Second, if enabled, provenance information is
-    calculated and stored as part of the result. Next, shape features are calculated on a cropped (no padding) version
-    of the original image. Then other featureclasses are calculated using all specified input images in ``inputImages``.
-    Images are cropped to tumor mask (no padding) after application of any filter and before being passed to the feature
-    class. Finally, the dictionary containing all calculated features is returned.
+    First, image and mask are loaded and normalized/resampled if necessary. Second, if enabled, provenance information
+    is calculated and stored as part of the result. Next, shape features are calculated on a cropped (no padding)
+    version of the original image. Then other featureclasses are calculated using all specified input images in
+    ``inputImages``. Images are cropped to tumor mask (no padding) after application of any filter and before being
+    passed to the feature class. Finally, the dictionary containing all calculated features is returned.
 
     :param imageFilepath: SimpleITK Image, or string pointing to image file location
     :param maskFilepath: SimpleITK Image, or string pointing to labelmap file location
@@ -360,7 +376,9 @@ def loadImage(self, ImageFilePath, MaskFilePath):
     All other cases are ignored (nothing calculated).
     Equal approach is used for assignment of mask using MaskFilePath.
 
-    If resampling is enabled, both image and mask are resampled and cropped to the tumormask (with additional
+    If normalizing is enabled image is first normalized before any resampling is applied.
+
+    If resampling is enabled, both image and mask are resampled and cropped to the tumor mask (with additional
     padding as specified in padDistance) after assignment of image and mask.
     """
     if isinstance(ImageFilePath, six.string_types) and os.path.exists(ImageFilePath):
@@ -381,6 +399,9 @@ def loadImage(self, ImageFilePath, MaskFilePath):
       if self.kwargs['verbose']: print("Error reading mask Filepath or SimpleITK object")
       mask = None
 
+    if self.kwargs['normalize']:
+      image = imageoperations.normalizeImage(image, self.kwargs['normalizeScale'], self.kwargs['removeOutliers'])
+
     if self.kwargs['interpolator'] is not None and self.kwargs['resampledPixelSpacing'] is not None:
       image, mask = imageoperations.resampleImage(image, mask,
                                                   self.kwargs['resampledPixelSpacing'],
@@ -436,6 +457,9 @@ def computeFeatures(self, image, mask, inputImageName, **kwargs):
     return featureVector
 
   def getFeatureClassNames(self):
+    """
+    Returns a list of all possible feature classes.
+    """
     return self.featureClasses.keys()
 
   def getFeatureNames(self, featureClassName):

radiomics/imageoperations.py

@@ -138,7 +138,8 @@ def cropToTumorMask(imageNode, maskNode, label=1, boundingBox=None):
 
 
 def resampleImage(imageNode, maskNode, resampledPixelSpacing, interpolator=sitk.sitkBSpline, label=1, padDistance=5):
-  """Resamples image or label to the specified pixel spacing (The default interpolator is Bspline)
+  """
+  Resamples image or label to the specified pixel spacing (The default interpolator is Bspline)
 
   'imageNode' is a SimpleITK Object, and 'resampledPixelSpacing' is the output pixel spacing (list of 3 elements).
 
@@ -231,6 +232,39 @@ def resampleImage(imageNode, maskNode, resampledPixelSpacing, interpolator=sitk.
 
   return resampledImageNode, resampledMaskNode
 
+def normalizeImage(image, scale=1, outliers=None):
+  r"""
+  Normalizes the image by centering it at the mean with standard deviation. Normalization is based on all gray values in
+  the image, not just those inside the segementation.
+
+  :math:`f(x) = \frac{s(x - \mu_x)}{\sigma_x}`
+
+  Where:
+
+  - :math:`x` and :math:`f(x)` are the original and normalized intensity, respectively.
+  - :math:`\mu_x` and :math:`\sigma_x` are the mean and standard deviation of the image instensity values.
+  - :math:`s` is an optional scaling defined by ``scale``. By default, it is set to 1.
+
+  Optionally, outliers can be removed, in which case values for which :math:`x > \mu_x + n\sigma_x` or
+  :math:`x < \mu_x - n\sigma_x` are set to :math:`\mu_x + n\sigma_x` and :math:`\mu_x - n\sigma_x`, respectively.
+  Here, :math:`n>0` and defined by ``outliers``. This, in turn, is controlled by the ``removeOutliers`` parameter.
+  Removal of outliers is done after the values of the image are normalized, but before ``scale`` is applied.
+  """
+  image = sitk.Normalize(image)
+
+  if outliers is not None:
+    imageArr = sitk.GetArrayFromImage(image)
+
+    imageArr[imageArr > outliers] = outliers
+    imageArr[imageArr < -outliers] = -outliers
+
+    newImage = sitk.GetImageFromArray(imageArr)
+    newImage.CopyInformation(image)
+
+  image *= scale
+
+  return image
+
 
 def applyThreshold(inputImage, lowerThreshold, upperThreshold, insideValue=None, outsideValue=0):
   # this mode is useful to generate the mask of thresholded voxels
@@ -309,32 +343,32 @@ def getLoGImage(inputImage, **kwargs):
 
 def getWaveletImage(inputImage, **kwargs):
   """
-      Apply wavelet filter to image and compute signature for each filtered image.
+  Apply wavelet filter to image and compute signature for each filtered image.
 
-      Following settings are possible:
+  Following settings are possible:
 
-      - start_level [0]: integer, 0 based level of wavelet which should be used as first set of decompositions
-        from which a signature is calculated
-      - level [1]: integer, number of levels of wavelet decompositions from which a signature is calculated.
-      - wavelet ["coif1"]: string, type of wavelet decomposition. Enumerated value, validated against possible values
-        present in the ``pyWavelet.wavelist()``. Current possible values (pywavelet version 0.4.0) (where an
-        aditional number is needed, range of values is indicated in []):
+  - start_level [0]: integer, 0 based level of wavelet which should be used as first set of decompositions
+    from which a signature is calculated
+  - level [1]: integer, number of levels of wavelet decompositions from which a signature is calculated.
+  - wavelet ["coif1"]: string, type of wavelet decomposition. Enumerated value, validated against possible values
+    present in the ``pyWavelet.wavelist()``. Current possible values (pywavelet version 0.4.0) (where an
+    aditional number is needed, range of values is indicated in []):
 
-        - haar
-        - dmey
-        - sym[2-20]
-        - db[1-20]
-        - coif[1-5]
-        - bior[1.1, 1.3, 1.5, 2.2, 2.4, 2.6, 2.8, 3.1, 3.3, 3.5, 3.7, 3.9, 4.4, 5.5, 6.8]
-        - rbio[1.1, 1.3, 1.5, 2.2, 2.4, 2.6, 2.8, 3.1, 3.3, 3.5, 3.7, 3.9, 4.4, 5.5, 6.8]
+    - haar
+    - dmey
+    - sym[2-20]
+    - db[1-20]
+    - coif[1-5]
+    - bior[1.1, 1.3, 1.5, 2.2, 2.4, 2.6, 2.8, 3.1, 3.3, 3.5, 3.7, 3.9, 4.4, 5.5, 6.8]
+    - rbio[1.1, 1.3, 1.5, 2.2, 2.4, 2.6, 2.8, 3.1, 3.3, 3.5, 3.7, 3.9, 4.4, 5.5, 6.8]
 
-      Returned filter name reflects wavelet type:
-      wavelet[level]-<decompositionName>
+  Returned filter name reflects wavelet type:
+  wavelet[level]-<decompositionName>
 
-      N.B. only levels greater than the first level are entered into the name.
+  N.B. only levels greater than the first level are entered into the name.
 
-      :return: Yields each wavelet decomposition and final approximation, corresponding filter name and ``kwargs``
-      """
+  :return: Yields each wavelet decomposition and final approximation, corresponding filter name and ``kwargs``
+  """
   global logger
 
   approx, ret = _swt3(inputImage, kwargs.get('wavelet', 'coif1'), kwargs.get('level', 1), kwargs.get('start_level', 0))