ENH: Implement extension in C for calculation of texture matrices.

Affects GLCM, GLRLM and GLSZM calculation. On initialization,
extension is loaded and enabled by default. If an error occurs,
a warning is logged and matrices will be calculated in
full-python mode.

Alternatively, full-python mode can be forced by a call
to `radiomics.enableCExtensions(False)`.

Added `requirements-setup.txt` to specify dependencies needed
during setup (requires numpy and cython)

Author: JoostJM

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,34 @@ 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 _cMatsState == 1:  # Extension loaded but not enabled
+      logger.info("Enabling C extensions")
+      _cMatsState = 2  # Enables matrix calculation in C
+    else:  # _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:
+    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 +130,18 @@ 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
+  _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.
+
 from ._version import get_versions
 
 __version__ = get_versions()['version']

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))