Merge pull request #202 from Radiomics/c-matrices_revised-jcfr-20170216

C matrices revised jcfr 20170216

MANIFEST.in

@@ -1,2 +1,17 @@
-include versioneer.py
+include CONTRIBUTING.md
+include LICENSE.txt
 include radiomics/_version.py
+include README.md
+include requirements.txt
+include requirements-dev.txt
+include radiomics/src/*.[ch]
+include versioneer.py
+
+recursive-include data *
+exclude data/PyradiomicsFeatures.csv data/Baseline2PyradiomicsFeaturesDiff.csv
+
+recursive-include tests *
+recursive-include docs *.rst conf.py Makefile make.bat
+
+recursive-exclude * __pycache__
+recursive-exclude * *.py[co]

data/paramSchema.yaml

@@ -8,6 +8,10 @@ mapping:
     mapping:
       verbose:
         type: bool
+      enableCExtensions:
+        type: bool
+      additionalInfo:
+        type: bool
       label:
         type: int
       binWidth:

radiomics/__init__.py

@@ -1,4 +1,3 @@
-
 # For convenience, import collection and numpy
 # into the "pyradiomics" namespace
 
@@ -8,6 +7,7 @@
 import os
 import pkgutil
 import sys
+import traceback
 
 import numpy  # noqa: F401
 
@@ -42,6 +42,36 @@ def debug(debug_on=True):
     debugging = False
 
 
+def enableCExtensions(enabled=True):
+  """
+  By default, calculation of GLCM, GLRLM and GLSZM is done in C, using extension ``_cmatrices.py``
+
+  If an error occurs during loading of this extension, a warning is logged and the extension is disabled,
+  matrices are then calculated in python.
+  The C extension can be disabled by calling this function as ``enableCExtensions(False)``, which forces the calculation
+  of the matrices to full-python mode.
+
+  Re-enabling use of C implementation is also done by this function, but if the extension is not loaded correctly,
+  a warning is logged and matrix calculation is forced to full-python mode.
+  """
+  global _cMatsState, logger
+  if enabled:
+    # If extensions are not yet enabled (_cMatsState == 2), check whether they are loaded (_cMatsState == 1) and if so,
+    # enable them. Otherwise, log a warning.
+    if _cMatsState == 1:  # Extension loaded but not enabled
+      logger.info("Enabling C extensions")
+      _cMatsState = 2  # Enables matrix calculation in C
+    elif _cMatsState == 0:  # _Extension not loaded correctly, do not enable matrix calculation in C and log warning
+      logger.warning("C Matrices not loaded correctly, cannot calculate matrices in C")
+  elif _cMatsState == 2:  # enabled = False, _cMatsState = 2: extensions currently enabled, disable them
+    logger.info("Disabling C extensions")
+    _cMatsState = 1
+
+
+def cMatsEnabled():
+  return _cMatsState == 2
+
+
 def getFeatureClasses():
   """
   Iterates over all modules of the radiomics package using pkgutil and subsequently imports those modules.
@@ -102,6 +132,20 @@ def getInputImageTypes():
 _featureClasses = None
 _inputImages = None
 
+# Indicates status of C extensions: 0 = not loaded, 1 = loaded but not enabled, 2 = enabled
+_cMatsState = 0
+
+try:
+  logger.debug("Loading C extensions")
+  from radiomics import _cmatrices as cMatrices
+  from radiomics import _cshape as cShape
+  _cMatsState = 1
+  enableCExtensions()
+except Exception:
+  logger.warning("Error loading C extensions, switching to python calculation:\n%s", traceback.format_exc())
+  cMatrices = None  # set cMatrices to None to prevent an import error in the feature classes.
+  cShape = None
+
 from ._version import get_versions
 
 __version__ = get_versions()['version']

radiomics/featureextractor.py

@@ -31,6 +31,8 @@ class RadiomicsFeaturesExtractor:
   Alternatively, at initialisation, the following general settings can be specified in kwargs:
 
   - 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.
@@ -85,6 +87,7 @@ def __init__(self, *args, **kwargs):
                      'padDistance': 5,
                      'label': 1,
                      'verbose': False,
+                     'enableCExtensions': True,
                      'additionalInfo': True}
       self.kwargs.update(kwargs)
 
@@ -173,6 +176,7 @@ def loadParams(self, paramsFile):
                    'padDistance': 5,
                    'label': 1,
                    'verbose': False,
+                   'enableCExtensions': True,
                    'additionalInfo': True}
     self.kwargs.update(kwargs)
 
@@ -295,6 +299,10 @@ def execute(self, imageFilepath, maskFilepath, label=None):
         is used. Default label is 1.
     :returns: dictionary containing calculated signature ("<filter>_<featureClass>_<featureName>":value).
     """
+    # Enable or disable C extensions for high performance matrix calculation. Only logs a message (INFO) when setting is
+    # successfully changed. If an error occurs, full-python mode is forced and a warning is logged.
+    radiomics.enableCExtensions(self.kwargs['enableCExtensions'])
+
     if label is not None:
       self.kwargs.update({'label': label})
 
@@ -303,7 +311,7 @@ def execute(self, imageFilepath, maskFilepath, label=None):
     featureVector = collections.OrderedDict()
     image, mask = self.loadImage(imageFilepath, maskFilepath)
 
-    if self.kwargs.get('additionalInfo'):
+    if self.kwargs['additionalInfo']:
       featureVector.update(self.getProvenance(imageFilepath, maskFilepath, mask))
 
     # Bounding box only needs to be calculated once after resampling, store the value, so it doesn't get calculated

radiomics/glcm.py

@@ -2,7 +2,7 @@
 from six.moves import range
 from tqdm import trange
 
-from radiomics import base, imageoperations
+from radiomics import base, cMatrices, cMatsEnabled, imageoperations
 
 
 class RadiomicsGLCM(base.RadiomicsFeaturesBase):
@@ -124,10 +124,14 @@ def __init__(self, inputImage, inputMask, **kwargs):
     self.matrix, self.histogram = imageoperations.binImage(self.binWidth, self.matrix, self.matrixCoordinates)
     self.coefficients['Ng'] = self.histogram[1].shape[0] - 1
 
-    self._calculateGLCM()
+    if cMatsEnabled():
+      self.P_glcm = self._calculateCMatrix()
+    else:
+      self.P_glcm = self._calculateMatrix()
+
     self._calculateCoefficients()
 
-  def _calculateGLCM(self):
+  def _calculateMatrix(self):
     r"""
     Compute GLCMs for the input image for every direction in 3D.
     Calculated GLCMs are placed in array P_glcm with shape (i/j, a)
@@ -142,7 +146,7 @@ def _calculateGLCM(self):
     size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
     angles = imageoperations.generateAngles(size)
 
-    self.P_glcm = numpy.zeros((Ng, Ng, int(angles.shape[0])), dtype='float64')
+    P_glcm = numpy.zeros((Ng, Ng, int(angles.shape[0])), dtype='float64')
 
     if self.verbose: bar = trange(Ng, desc='calculate GLCM')
 
@@ -167,12 +171,32 @@ def _calculateGLCM(self):
           # that are also a neighbour of a voxel with gray level i for angle a.
           # The number of indices is then equal to the total number of pairs with gray level i and j for angle a
           count = len(neighbour_indices.intersection(j_indices))
-          self.P_glcm[i - 1, j - 1, a_idx] = count
+          P_glcm[i - 1, j - 1, a_idx] = count
     if self.verbose: bar.close()
 
+    P_glcm = self._applyMatrixOptions(P_glcm, angles)
+
+    return P_glcm
+
+  def _calculateCMatrix(self):
+    size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
+    angles = imageoperations.generateAngles(size)
+    Ng = self.coefficients['Ng']
+
+    P_glcm = cMatrices.calculate_glcm(self.matrix, self.maskArray, angles, Ng)
+    P_glcm = self._applyMatrixOptions(P_glcm, angles)
+
+    return P_glcm
+
+  def _applyMatrixOptions(self, P_glcm, angles):
+    """
+    Further process calculated matrix by optionally making it symmetrical and/or applying a weighting factor.
+    Finally, delete empty angles and normalize the GLCM by dividing it by the sum of its elements.
+    """
+
     # Optionally make GLCMs symmetrical for each angle
     if self.symmetricalGLCM:
-      self.P_glcm += numpy.transpose(self.P_glcm, (1, 0, 2))
+      P_glcm += numpy.transpose(P_glcm, (1, 0, 2))
 
     # Optionally apply a weighting factor
     if self.weightingNorm is not None:
@@ -191,17 +215,17 @@ def _calculateGLCM(self):
           self.logger.warning('weigthing norm "%s" is unknown, W is set to 1', self.weightingNorm)
           weights[a_idx] = 1
 
-      self.P_glcm = numpy.sum(self.P_glcm * weights[None, None, :], 2, keepdims=True)
+      P_glcm = numpy.sum(P_glcm * weights[None, None, :], 2, keepdims=True)
 
-    sumP_glcm = numpy.sum(self.P_glcm, (0, 1), keepdims=True)  # , keepdims=True)
+    sumP_glcm = numpy.sum(P_glcm, (0, 1), keepdims=True)
 
     # Delete empty angles if no weighting is applied
-    if self.P_glcm.shape[2] > 1:
-      self.P_glcm = numpy.delete(self.P_glcm, numpy.where(sumP_glcm == 0), 2)
+    if P_glcm.shape[2] > 1:
+      P_glcm = numpy.delete(P_glcm, numpy.where(sumP_glcm == 0), 2)
       sumP_glcm = numpy.delete(sumP_glcm, numpy.where(sumP_glcm == 0), 2)
 
     # Normalize each glcm
-    self.P_glcm = self.P_glcm / sumP_glcm
+    return P_glcm / sumP_glcm
 
   # check if ivector and jvector can be replaced
   def _calculateCoefficients(self):

radiomics/glrlm.py

@@ -3,7 +3,7 @@
 import numpy
 from six.moves import range
 
-from radiomics import base, imageoperations
+from radiomics import base, cMatrices, cMatsEnabled, imageoperations
 
 
 class RadiomicsGLRLM(base.RadiomicsFeaturesBase):
@@ -93,10 +93,14 @@ def __init__(self, inputImage, inputMask, **kwargs):
     self.coefficients['Nr'] = numpy.max(self.matrix.shape)
     self.coefficients['Np'] = self.targetVoxelArray.size
 
-    self._calculateGLRLM()
+    if cMatsEnabled():
+      self.P_glrlm = self._calculateCMatrix()
+    else:
+      self.P_glrlm = self._calculateMatrix()
+
     self._calculateCoefficients()
 
-  def _calculateGLRLM(self):
+  def _calculateMatrix(self):
     Ng = self.coefficients['Ng']
     Nr = self.coefficients['Nr']
 
@@ -152,11 +156,34 @@ def _calculateGLRLM(self):
       if not isMultiElement:
         P[:] = 0
 
+    P_glrlm = self._applyMatrixOptions(P_glrlm, angles)
+
+    return P_glrlm
+
+  def _calculateCMatrix(self):
+    Ng = self.coefficients['Ng']
+    Nr = self.coefficients['Nr']
+
+    size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
+    angles = imageoperations.generateAngles(size)
+
+    P_glrlm = cMatrices.calculate_glrlm(self.matrix, self.maskArray, angles, Ng, Nr)
+    P_glrlm = self._applyMatrixOptions(P_glrlm, angles)
+
+    return P_glrlm
+
+  def _applyMatrixOptions(self, P_glrlm, angles):
+    """
+    Further process the calculated matrix by cropping the matrix to between minimum and maximum observed gray-levels and
+    up to maximum observed run-length. Optionally apply a weighting factor. Finally delete empty angles and store the
+    sum of the matrix in ``self.coefficients``.
+    """
+
     # Crop gray-level axis of GLRLMs to between minimum and maximum observed gray-levels
     # Crop run-length axis of GLRLMs up to maximum observed run-length
     P_glrlm_bounds = numpy.argwhere(P_glrlm)
     (xstart, ystart, zstart), (xstop, ystop, zstop) = P_glrlm_bounds.min(0), P_glrlm_bounds.max(0) + 1  # noqa: F841
-    self.P_glrlm = P_glrlm[xstart:xstop, :ystop, :]
+    P_glrlm = P_glrlm[xstart:xstop, :ystop, :]
 
     # Optionally apply a weighting factor
     if self.weightingNorm is not None:
@@ -175,17 +202,19 @@ def _calculateGLRLM(self):
           self.logger.warning('weigthing norm "%s" is unknown, weighting factor is set to 1', self.weightingNorm)
           weights[a_idx] = 1
 
-      self.P_glrlm = numpy.sum(self.P_glrlm * weights[None, None, :], 2, keepdims=True)
+      P_glrlm = numpy.sum(P_glrlm * weights[None, None, :], 2, keepdims=True)
 
-    sumP_glrlm = numpy.sum(self.P_glrlm, (0, 1))
+    sumP_glrlm = numpy.sum(P_glrlm, (0, 1))
 
     # Delete empty angles if no weighting is applied
-    if self.P_glrlm.shape[2] > 1:
-      self.P_glrlm = numpy.delete(self.P_glrlm, numpy.where(sumP_glrlm == 0), 2)
+    if P_glrlm.shape[2] > 1:
+      P_glrlm = numpy.delete(P_glrlm, numpy.where(sumP_glrlm == 0), 2)
       sumP_glrlm = numpy.delete(sumP_glrlm, numpy.where(sumP_glrlm == 0), 0)
 
     self.coefficients['sumP_glrlm'] = sumP_glrlm
 
+    return P_glrlm
+
   def _calculateCoefficients(self):
 
     pr = numpy.sum(self.P_glrlm, 0)

radiomics/glszm.py

@@ -2,7 +2,7 @@
 from six.moves import range
 from tqdm import trange
 
-from radiomics import base, imageoperations
+from radiomics import base, cMatrices, cMatsEnabled, imageoperations
 
 
 class RadiomicsGLSZM(base.RadiomicsFeaturesBase):
@@ -66,10 +66,14 @@ def __init__(self, inputImage, inputMask, **kwargs):
     self.coefficients['Ng'] = self.histogram[1].shape[0] - 1
     self.coefficients['Np'] = self.targetVoxelArray.size
 
-    self._calculateGLSZM()
+    if cMatsEnabled():
+      self.P_glszm = self._calculateCMatrix()
+    else:
+      self.P_glszm = self._calculateMatrix()
+
     self._calculateCoefficients()
 
-  def _calculateGLSZM(self):
+  def _calculateMatrix(self):
     """
     Number of times a region with a
     gray level and voxel count occurs in an image. P_glszm[level, voxel_count] = # occurrences
@@ -131,7 +135,15 @@ def _calculateGLSZM(self):
     # Crop size-zone area axis of GLSZM matrix up to maximum observed size-zone area
     P_glszm_bounds = numpy.argwhere(P_glszm)
     (xstart, ystart), (xstop, ystop) = P_glszm_bounds.min(0), P_glszm_bounds.max(0) + 1  # noqa: F841
-    self.P_glszm = P_glszm[xstart:xstop, :ystop]
+    return P_glszm[xstart:xstop, :ystop]
+
+  def _calculateCMatrix(self):
+    size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
+    angles = imageoperations.generateAngles(size)
+    Ng = self.coefficients['Ng']
+    Ns = self.coefficients['Np']
+
+    return cMatrices.calculate_glszm(self.matrix, self.maskArray, angles, Ng, Ns)
 
   def _calculateCoefficients(self):
     sumP_glszm = numpy.sum(self.P_glszm, (0, 1))