Merge pull request #302 from JoostJM/add-resegementation

Add optional resegmentation of mask

docs/customization.rst

@@ -154,6 +154,17 @@ Feature Extractor Level
 .. note::
     Resampling is disabled when either `resampledPixelSpacing` or `interpolator` is set to `None`
 
+*Resegmentation*
+
+- resegmentRange [None]: List of 2 floats, specifies the lower and upper threshold, respectively. Segmented voxels
+  outside this range are removed from the mask prior to feature calculation. When the value is None (default), no
+  resegmentation is performed. Resegemented size is checked (using parameter ``minimumROISize``, default 1) and upon
+  fail, an error is logged and the mask is reset to the original mask.
+
+.. note::
+    This only affects first order and texture classes. No resegmentation is performed prior to calculating shape
+    features.
+
 *Mask validation*
 
 - minimumROIDimensions [1]: Integer, range 1-3, specifies the minimum dimensions (1D, 2D or 3D, respectively).
@@ -238,6 +249,7 @@ Feature Class Level
     For more information on how weighting is applied, see the documentation on :ref:`GLCM <radiomics-glcm-label>` and
     :ref:`GLRLM <radiomics-glszm-label>`.
 
+
 Feature Class Specific Settings
 +++++++++++++++++++++++++++++++
 

radiomics/featureextractor.py

@@ -89,6 +89,7 @@ class specific are defined in the respective feature classes and and not include
             'distances': [1],
             'force2D': False,
             'force2Ddimension': 0,
+            'resegmentRange': None,  # No resegmentation by default
             'label': 1,
             'enableCExtensions': True,
             'additionalInfo': True}
@@ -291,10 +292,12 @@ def execute(self, imageFilepath, maskFilepath, label=None):
        bounding box.
     3. If enabled, provenance information is calculated and stored as part of the result.
     4. Shape features are calculated on a cropped (no padding) version of the original image.
-    5. Other enabled feature classes are calculated using all specified image types in ``_enabledImageTypes``. Images
+    5. If enabled, resegment the mask based upon the range specified in ``resegmentRange`` (default None: resegmentation
+       disabled).
+    6. Other enabled feature classes are calculated using all specified image types in ``_enabledImageTypes``. Images
        are cropped to tumor mask (no padding) after application of any filter and before being passed to the feature
        class.
-    6. The calculated features is returned as ``collections.OrderedDict``.
+    7. The calculated features is returned as ``collections.OrderedDict``.
 
     :param imageFilepath: SimpleITK Image, or string pointing to image file location
     :param maskFilepath: SimpleITK Image, or string pointing to labelmap file location
@@ -359,7 +362,25 @@ def execute(self, imageFilepath, maskFilepath, label=None):
         newFeatureName = 'original_shape_%s' % featureName
         featureVector[newFeatureName] = featureValue
 
-    # 5. Calculate other enabled feature classes using enabled image types
+    # 5. Resegment the mask if enabled (parameter regsegmentMask is not None)
+    resegmentRange = self.settings.get('resegmentRange', None)
+    if resegmentRange is not None:
+      resegmentedMask = imageoperations.resegmentMask(image, mask, resegmentRange, self.settings['label'])
+
+      # Recheck to see if the mask is still valid
+      boundingBox, correctedMask = imageoperations.checkMask(image, resegmentedMask, **self.settings)
+      # Update the mask if it had to be resampled
+      if correctedMask is not None:
+        resegmentedMask = correctedMask
+
+      if boundingBox is None:
+        # Mask checks failed, do not extract features and return the empty featureVector
+        return featureVector
+
+      # Resegmentation successful
+      mask = resegmentedMask
+
+    # 6. Calculate other enabled feature classes using enabled image types
     # Make generators for all enabled image types
     self.logger.debug('Creating image type iterator')
     imageGenerators = []

radiomics/imageoperations.py

@@ -251,15 +251,15 @@ def checkMask(imageNode, maskNode, **kwargs):
     correctedMask = _correctMask(imageNode, maskNode, label)
     if correctedMask is None:  # Resampling failed (ROI outside image physical space
       logger.error('Image/Mask correction failed, ROI invalid (not found or outside of physical image bounds)')
-      return (boundingBox, correctedMask)
+      return boundingBox, correctedMask
 
     # Resampling succesful, try to calculate boundingbox
     try:
       lsif.Execute(imageNode, correctedMask)
     except RuntimeError:
       logger.error('Calculation of bounding box failed, for more information run with DEBUG logging and check log')
       logger.debug('Bounding box calculation with resampled mask failed', exc_info=True)
-      return (boundingBox, correctedMask)
+      return boundingBox, correctedMask
 
   # LBound and UBound of the bounding box, as (L_X, U_X, L_Y, U_Y, L_Z, U_Z)
   boundingBox = numpy.array(lsif.GetBoundingBox(label))
@@ -268,16 +268,16 @@ def checkMask(imageNode, maskNode, **kwargs):
   ndims = numpy.sum((boundingBox[1::2] - boundingBox[0::2] + 1) > 1)  # UBound - LBound + 1 = Size
   if ndims <= minDims:
     logger.error('mask has too few dimensions (number of dimensions %d, minimum required %d)', ndims, minDims)
-    return (boundingBox, correctedMask)
+    return None, correctedMask
 
   if minSize is not None:
     logger.debug('Checking minimum size requirements (minimum size: %d)', minSize)
     roiSize = lsif.GetCount(label)
     if roiSize <= minSize:
       logger.error('Size of the ROI is too small (minimum size: %g, ROI size: %g', minSize, roiSize)
-      return (boundingBox, correctedMask)
+      return None, correctedMask
 
-  return (boundingBox, correctedMask)
+  return boundingBox, correctedMask
 
 
 def _correctMask(imageNode, maskNode, label):
@@ -575,19 +575,41 @@ def normalizeImage(image, scale=1, outliers=None):
   return image
 
 
-def applyThreshold(inputImage, lowerThreshold, upperThreshold, insideValue=None, outsideValue=0):
-  # this mode is useful to generate the mask of thresholded voxels
-  if insideValue:
-    tif = sitk.BinaryThresholdImageFilter()
-    tif.SetInsideValue(insideValue)
-    tif.SetLowerThreshold(lowerThreshold)
-    tif.SetUpperThreshold(upperThreshold)
-  else:
-    tif = sitk.ThresholdImageFilter()
-    tif.SetLower(lowerThreshold)
-    tif.SetUpper(upperThreshold)
-  tif.SetOutsideValue(outsideValue)
-  return tif.Execute(inputImage)
+def resegmentMask(imageNode, maskNode, resegmentRange, label=1):
+  """
+  Resegment the Mask based on the range specified in ``resegmentRange``. All voxels with a gray level outside the
+  range specified are removed from the mask. The resegmented mask is therefore always equal or smaller in size than
+  the original mask. The resegemented mask is then checked for size (as specified by parameter ``minimumROISize``,
+  defaults to minimum size 1). When this check fails, an error is logged and maskArray is reset to the original mask.
+  """
+  global logger
+
+  if resegmentRange is None:
+    return maskNode
+
+  logger.debug('Resegmenting mask (range %s)', resegmentRange)
+
+  im_arr = sitk.GetArrayFromImage(imageNode)
+  ma_arr = (sitk.GetArrayFromImage(maskNode) == label)  # boolean array
+
+  oldSize = numpy.sum(ma_arr)
+  # Get a boolean array specifying per voxel if its gray value is inside the range specified by resegmentRange.
+  intensitySegmentation = numpy.logical_and(im_arr >= resegmentRange[0],
+                                            im_arr <= resegmentRange[1])
+  # Then, take the intersection between voxels segmented (ma_arr) and
+  # voxels inside specified range (intensitySegmentation)
+  ma_arr = numpy.logical_and(ma_arr, intensitySegmentation)
+  roiSize = numpy.sum(ma_arr)
+
+  # Transform the boolean array back to an image with the correct voxels set to the label value
+  newMask_arr = numpy.zeros(ma_arr.shape, dtype='int')
+  newMask_arr[ma_arr] = label
+
+  newMask = sitk.GetImageFromArray(newMask_arr)
+  newMask.CopyInformation(maskNode)
+  logger.debug('Resegmentation complete, new size: %d voxels (excluded %d voxels)', roiSize, oldSize - roiSize)
+
+  return newMask
 
 
 def getOriginalImage(inputImage, **kwargs):

radiomics/schemas/paramSchema.yaml

@@ -65,6 +65,9 @@ mapping:
         range:
           min: 0
           max: 2
+      resegmentRange:
+        seq:
+          - type: float
       sigma:
         seq:
           - type: float