Update c extensions

radiomics/featureextractor.py

@@ -428,18 +428,28 @@ def execute(self, imageFilepath, maskFilepath, label=None, voxelBased=False):
 
     if not voxelBased:
       # 4. If shape descriptors should be calculated, handle it separately here
+      Nd = mask.GetDimension()
       if 'shape' in self._enabledFeatures.keys():
-        featureVector.update(self.computeShape(image, mask, boundingBox))
-      if 'shape2D' in self._enabledFeatures.keys():
-        force2D = self.settings.get('force2D', False)
-        force2Ddimension = self.settings.get('force2Ddimension', 0)
-        if not force2D:
-          self.logger.warning('parameter force2D must be set to True to enable shape2D extraction')
-        elif not (boundingBox[1::2] - boundingBox[0::2] + 1)[force2Ddimension] > 1:
-          self.logger.warning('Size in specified 2D dimension (%i) is greater than 1, cannot calculate 2D shape',
-                              force2Ddimension)
+        if Nd == 3:
+          featureVector.update(self.computeShape(image, mask, boundingBox))
         else:
+          self.logger.warning('Shape features are only available 3D input (for 2D input, use shape2D). Found %iD input', Nd)
+      if 'shape2D' in self._enabledFeatures.keys():
+        if Nd == 3:
+          force2D = self.settings.get('force2D', False)
+          force2Ddimension = self.settings.get('force2Ddimension', 0)
+          if not force2D:
+            self.logger.warning('parameter force2D must be set to True to enable shape2D extraction')
+          elif not (boundingBox[1::2] - boundingBox[0::2] + 1)[force2Ddimension] > 1:
+            self.logger.warning('Size in specified 2D dimension (%i) is greater than 1, cannot calculate 2D shape',
+                                force2Ddimension)
+          else:
+            featureVector.update(self.computeShape(image, mask, boundingBox, 'shape2D'))
+        elif Nd == 2:
           featureVector.update(self.computeShape(image, mask, boundingBox, 'shape2D'))
+        else:
+          self.logger.warning('Shape2D features are only available for 2D and 3D (with force2D=True) input. '
+                              'Found %iD input', Nd)
 
     # (Default) Only use resegemented mask for feature classes other than shape
     # can be overridden by specifying `resegmentShape` = True

radiomics/firstorder.py

@@ -103,8 +103,7 @@ def getTotalEnergyFeatureValue(self):
       Not present in IBSI feature definitions
     """
 
-    x, y, z = self.pixelSpacing
-    cubicMMPerVoxel = x * y * z
+    cubicMMPerVoxel = numpy.multiply.reduce(self.pixelSpacing)
 
     return cubicMMPerVoxel * self.getEnergyFeatureValue()
 

radiomics/generalinfo.py

@@ -49,8 +49,9 @@ def addImageElements(self, image, prefix='original'):
 
     Adds the following:
 
-    - ImageHash: sha1 hash of the mask, which can be used to check if the same mask was used during reproducibility
+    - Hash: sha1 hash of the mask, which can be used to check if the same mask was used during reproducibility
       tests. (Only added when prefix is "original")
+    - Dimensionality: Number of dimensions (e.g. 2D, 3D) in the image. (Only added when prefix is "original")
     - Spacing: Pixel spacing (x, y, z) in mm.
     - Size: Dimensions (x, y, z) of the image in number of voxels.
     - Mean: Mean intensity value over all voxels in the image.
@@ -64,6 +65,7 @@ def addImageElements(self, image, prefix='original'):
     """
     if prefix == 'original':
       self.generalInfo[self.generalInfo_prefix + 'Image-original_Hash'] = sitk.Hash(image)
+      self.generalInfo[self.generalInfo_prefix + 'Image-original_Dimensionality'] = '%iD' % image.GetDimension()
 
     self.generalInfo[self.generalInfo_prefix + 'Image-' + prefix + '_Spacing'] = image.GetSpacing()
     self.generalInfo[self.generalInfo_prefix + 'Image-' + prefix + '_Size'] = image.GetSize()

radiomics/imageoperations.py

@@ -277,7 +277,7 @@ def _checkROI(imageNode, maskNode, **kwargs):
 
   For the second check, a tolerance of 1e-3 is allowed.
 
-  If the ROI is valid, the bounding box (lower bounds and size in 3 directions: L_X, L_Y, L_Z, S_X, S_Y, S_Z) is
+  If the ROI is valid, the bounding box (lower bounds, followd by size in all dimensions (X, Y, Z ordered)) is
   returned. Otherwise, a ValueError is raised.
   """
   global logger
@@ -293,17 +293,18 @@ def _checkROI(imageNode, maskNode, **kwargs):
   if label not in lssif.GetLabels():
     raise ValueError('Label (%d) not present in mask', label)
 
-  # LBound and size of the bounding box, as (L_X, L_Y, L_Z, S_X, S_Y, S_Z)
+  # LBound and size of the bounding box, as (L_X, L_Y, [L_Z], S_X, S_Y, [S_Z])
   bb = numpy.array(lssif.GetBoundingBox(label))
+  Nd = maskNode.GetDimension()
 
   # Determine if the ROI is within the physical space of the image
 
   logger.debug('Comparing physical space of bounding box to physical space of image')
   # Step 1: Get the origin and UBound corners of the bounding box in physical space
   # The additional 0.5 represents the difference between the voxel center and the voxel corner
-  # Upper bound index of ROI = bb[:3] + bb[3:] - 1 (LBound + Size - 1), .5 is added to get corner
-  ROIBounds = (maskNode.TransformContinuousIndexToPhysicalPoint(bb[:3] - .5),  # Origin
-               maskNode.TransformContinuousIndexToPhysicalPoint(bb[:3] + bb[3:] - 0.5))  # UBound
+  # Upper bound index of ROI = bb[:Nd] + bb[Nd:] - 1 (LBound + Size - 1), .5 is added to get corner
+  ROIBounds = (maskNode.TransformContinuousIndexToPhysicalPoint(bb[:Nd] - .5),  # Origin
+               maskNode.TransformContinuousIndexToPhysicalPoint(bb[:Nd] + bb[Nd:] - 0.5))  # UBound
   # Step 2: Translate the ROI physical bounds to the image coordinate space
   ROIBounds = (imageNode.TransformPhysicalPointToContinuousIndex(ROIBounds[0]),  # Origin
                imageNode.TransformPhysicalPointToContinuousIndex(ROIBounds[1]))
@@ -415,6 +416,11 @@ def resampleImage(imageNode, maskNode, **kwargs):
   maskSpacing = numpy.array(maskNode.GetSpacing())
   imageSpacing = numpy.array(imageNode.GetSpacing())
 
+  Nd_resampled = len(resampledPixelSpacing)
+  Nd_mask = len(maskSpacing)
+  assert Nd_resampled == Nd_mask, \
+      'Wrong dimensionality (%i-D) of resampledPixelSpacing!, %i-D required' % (Nd_resampled, Nd_mask)
+
   # If spacing for a direction is set to 0, use the original spacing (enables "only in-slice" resampling)
   logger.debug('Where resampled spacing is set to 0, set it to the original spacing (mask)')
   resampledPixelSpacing = numpy.array(resampledPixelSpacing)
@@ -426,7 +432,7 @@ def resampleImage(imageNode, maskNode, **kwargs):
 
   # Do not resample in those directions where labelmap spans only one slice.
   maskSize = numpy.array(maskNode.GetSize())
-  resampledPixelSpacing = numpy.where(bb[3:] != 1, resampledPixelSpacing, maskSpacing)
+  resampledPixelSpacing = numpy.where(bb[Nd_mask:] != 1, resampledPixelSpacing, maskSpacing)
 
   # If current spacing is equal to resampledPixelSpacing, no interpolation is needed
   # Tolerance = 1e-5 + 1e-8*abs(resampledSpacing)
@@ -440,8 +446,8 @@ def resampleImage(imageNode, maskNode, **kwargs):
   # Determine bounds of cropped volume in terms of new Index coordinate space,
   # round down for lowerbound and up for upperbound to ensure entire segmentation is captured (prevent data loss)
   # Pad with an extra .5 to prevent data loss in case of upsampling. For Ubound this is (-1 + 0.5 = -0.5)
-  bbNewLBound = numpy.floor((bb[:3] - 0.5) * spacingRatio - padDistance)
-  bbNewUBound = numpy.ceil((bb[:3] + bb[3:] - 0.5) * spacingRatio + padDistance)
+  bbNewLBound = numpy.floor((bb[:Nd_mask] - 0.5) * spacingRatio - padDistance)
+  bbNewUBound = numpy.ceil((bb[:Nd_mask] + bb[Nd_mask:] - 0.5) * spacingRatio + padDistance)
 
   # Ensure resampling is not performed outside bounds of original image
   maxUbound = numpy.ceil(maskSize * spacingRatio) - 1
@@ -741,7 +747,8 @@ def getWaveletImage(inputImage, inputMask, **kwargs):
 
   logger.debug('Generating Wavelet images')
 
-  axes = [2, 1, 0]
+  Nd = inputImage.GetDimension()
+  axes = list(range(Nd - 1, -1, -1))
   if kwargs.get('force2D', False):
     axes.remove(kwargs.get('force2Ddimension', 0))
 
@@ -773,11 +780,10 @@ def _swt3(inputImage, axes, **kwargs):  # Stationary Wavelet Transform 3D
 
   matrix = sitk.GetArrayFromImage(inputImage)  # This function gets a numpy array from the SimpleITK Image "inputImage"
   matrix = numpy.asarray(matrix) # The function np.asarray converts "matrix" (which could be also a tuple) into an array.
-  if matrix.ndim != 3:
-    raise ValueError('Expected 3D data array')
 
   original_shape = matrix.shape
   # original_shape becomes a tuple (?,?,?) containing the number of rows, columns, and slices of the image
+  # this is of course dependent on the number of dimensions, but the same principle holds
   padding = tuple([(0, 1 if dim % 2 != 0 else 0) for dim in original_shape])
   # padding is necessary because of pywt.swtn (see function Notes)
   data = matrix.copy()  # creates a modifiable copy of "matrix" and we call it "data"
@@ -786,13 +792,17 @@ def _swt3(inputImage, axes, **kwargs):  # Stationary Wavelet Transform 3D
   if not isinstance(wavelet, pywt.Wavelet):
     wavelet = pywt.Wavelet(wavelet)
 
-  for i in range(0, start_level):  # if start_level = 0 this for loop never gets executed
-    dec = pywt.swtn(data, wavelet, level=1, start_level=0, axes=axes)[0] # computes all decompositions and saves them in "dec" dict
-    data = dec['a' * len(axes)].copy()  # copies in "data" just the "aaa" decomposition (if len(axes) = 3)
+  for i in range(0, start_level):  # if start_level = 0 (default) this for loop never gets executed
+    # compute all decompositions and saves them in "dec" dict
+    dec = pywt.swtn(data, wavelet, level=1, start_level=0, axes=axes)[0]
+    # copies in "data" just the "aaa" decomposition (i.e. approximation; No of consecutive 'a's = len(axes))
+    data = dec['a' * len(axes)].copy()
 
   ret = []  # initialize empty list
   for i in range(start_level, start_level + level):
-    dec = pywt.swtn(data, wavelet, level=1, start_level=0, axes=axes)[0]  # computes the n-dimensional stationary wavelet transform
+    # compute the n-dimensional stationary wavelet transform
+    dec = pywt.swtn(data, wavelet, level=1, start_level=0, axes=axes)[0]
+    # Copy the approximation into data (approximation in output / input for next levels)
     data = dec['a' * len(axes)].copy()
 
     dec_im = {}  # initialize empty dict
@@ -984,22 +994,28 @@ def getLBP2DImage(inputImage, inputMask, **kwargs):
     logger.warning('Could not load required package "skimage", cannot implement filter LBP 2D')
     return
 
-  # Warn the user if features are extracted in 3D, as this function calculates LBP in 2D
-  if not kwargs.get('force2D', False):
-    logger.warning('Calculating Local Binary Pattern in 2D, but extracting features in 3D. Use with caution!')
-
-  lbp_axis = kwargs.get('force2Ddimension', 0)
-
   lbp_radius = kwargs.get('lbp2DRadius', 1)
   lbp_samples = kwargs.get('lbp2DSamples', 8)
   lbp_method = kwargs.get('lbp2DMethod', 'uniform')
 
   im_arr = sitk.GetArrayFromImage(inputImage)
 
-  im_arr = im_arr.swapaxes(0, lbp_axis)
-  for idx in range(im_arr.shape[0]):
-    im_arr[idx, ...] = local_binary_pattern(im_arr[idx, ...], P=lbp_samples, R=lbp_radius, method=lbp_method)
-  im_arr = im_arr.swapaxes(0, lbp_axis)
+  Nd = inputImage.GetDimension()
+  if Nd == 3:
+    # Warn the user if features are extracted in 3D, as this function calculates LBP in 2D
+    if not kwargs.get('force2D', False):
+      logger.warning('Calculating Local Binary Pattern in 2D, but extracting features in 3D. Use with caution!')
+    lbp_axis = kwargs.get('force2Ddimension', 0)
+
+    im_arr = im_arr.swapaxes(0, lbp_axis)
+    for idx in range(im_arr.shape[0]):
+      im_arr[idx, ...] = local_binary_pattern(im_arr[idx, ...], P=lbp_samples, R=lbp_radius, method=lbp_method)
+    im_arr = im_arr.swapaxes(0, lbp_axis)
+  elif Nd == 2:
+    im_arr = local_binary_pattern(im_arr, P=lbp_samples, R=lbp_radius, method=lbp_method)
+  else:
+    logger.warning('LBP 2D is only available for 2D or 3D with forced 2D extraction')
+    return
 
   im = sitk.GetImageFromArray(im_arr)
   im.CopyInformation(inputImage)
@@ -1036,6 +1052,11 @@ def getLBP3DImage(inputImage, inputMask, **kwargs):
     Science, vol 7728. Springer, Berlin, Heidelberg. doi:10.1007/978-3-642-37410-4_3
   """
   global logger
+  Nd = inputImage.GetDimension()
+  if Nd != 3:
+    logger.warning('LBP 3D only available for 3 dimensional images, found %i dimensions', Nd)
+    return
+
   try:
     from scipy.stats import kurtosis
     from scipy.ndimage.interpolation import map_coordinates

radiomics/shape.py

@@ -38,6 +38,7 @@ class RadiomicsShape(base.RadiomicsFeaturesBase):
   """
 
   def __init__(self, inputImage, inputMask, **kwargs):
+    assert inputMask.GetDimension() == 3, 'Shape features are only available in 3D. If 2D, use shape2D instead'
     super(RadiomicsShape, self).__init__(inputImage, inputMask, **kwargs)
 
   def _initVoxelBasedCalculation(self):

radiomics/shape2D.py

@@ -52,41 +52,32 @@ def _initVoxelBasedCalculation(self):
     raise NotImplementedError('Shape features are not available in pixel-based mode')
 
   def _initSegmentBasedCalculation(self):
-    if not self.settings.get('force2D', False):
-      raise ValueError('Shape2D is can only be calculated when `force2D` is True')
+    self.maskArray = (sitk.GetArrayFromImage(self.inputMask) == self.label)  # boolean array
 
-    force2DDimension = self.settings.get('force2Ddimension', 0)
-    axes = [0, 1, 2]
-    axes.remove(force2DDimension)
+    Nd = self.inputMask.GetDimension()
+    if Nd == 3:
+      if not self.settings.get('force2D', False):
+        raise ValueError('Shape2D is can only be calculated when input is 2D or 3D with `force2D=True`')
 
-    self.pixelSpacing = numpy.array(self.inputImage.GetSpacing()[::-1])[(axes,)]
+      force2DDimension = self.settings.get('force2Ddimension', 0)
+      axes = [0, 1, 2]
+      axes.remove(force2DDimension)
 
-    # Pad inputMask to prevent index-out-of-range errors
-    self.logger.debug('Padding the mask with 0s')
-
-    cpif = sitk.ConstantPadImageFilter()
-
-    padding = numpy.tile(1, 3)
-    padding[force2DDimension] = 0  # Do not pad in 2D dimension
-    padding = padding[::-1]  # Reverse order (ITK is x, y, z; Numpy is z, y, x)
-    try:
-      cpif.SetPadLowerBound(padding)
-      cpif.SetPadUpperBound(padding)
-    except TypeError:
-      # newer versions of SITK/python want a tuple or list
-      cpif.SetPadLowerBound(padding.tolist())
-      cpif.SetPadUpperBound(padding.tolist())
+      self.pixelSpacing = numpy.array(self.inputImage.GetSpacing()[::-1])[(axes,)]
 
-    self.inputMask = cpif.Execute(self.inputMask)
+      if self.maskArray.shape[force2DDimension] > 1:
+        raise ValueError('Size of the mask in dimension %i is more than 1, cannot compute 2D shape')
 
-    # Reassign self.maskArray using the now-padded self.inputMask
-    self.maskArray = numpy.array(sitk.GetArrayFromImage(self.inputMask) == self.label)
+      # Drop the 2D axis, ensuring the input is truly 2D
+      self.maskArray = numpy.squeeze(self.maskArray, axis=force2DDimension)
+    elif Nd == 2:
+      self.pixelSpacing = numpy.array(self.inputImage.GetSpacing()[::-1])
+    else:
+      raise ValueError('Shape2D is can only be calculated when input is 2D or 3D with `force2D=True`')
 
-    if self.maskArray.shape[force2DDimension] > 1:
-      raise ValueError('Size of the mask in dimension %i is more than 1, cannot compute 2D shape')
-
-    # Drop the 2D axis, ensuring the input is truly 2D
-    self.maskArray = numpy.squeeze(self.maskArray, axis=force2DDimension)
+    # Pad maskArray to prevent index-out-of-range errors
+    self.logger.debug('Padding the mask with 0s')
+    self.maskArray = numpy.pad(self.maskArray, pad_width=1, mode='constant', constant_values=0)
 
     self.labelledPixelCoordinates = numpy.where(self.maskArray != 0)