ENH: Support for faster C-based matrix and shape surface computation

Stats (Profiling 5 testcases):
  GLCM 9653 ms -> 7 ms
  GLDM 348 ms -> 7 ms
  NGTDM 133 ms -> 4ms

This commit adds two C Python extensions and associated tests:

* _cmatrices: C implementation for matrix computation associayed
  with GLCM, GLDM, GLDZM, GLRLM, GLSZM. and NGTDM features

* _cshape: C implementation for Shape surface computation

* tests/test_cmatrices: testing for matrix equality: Tests
  whether the python generated matrix is equal to the matrix
  generated by the C implementations (allows for machine precision errors)

Details
-------

* Add docstring to C modules

* Use of C implementation optional. At initialization, the package
  tries to use C, but if loading has failed, or calculation is forced
  to python, python is used. Note that the import of _cmatrices is done
  after initialization of logger. This ensures error will be logged
  if import fails.

* GLDM, NGTDM: C implementation accepts only one
  direction set of angles (i.e. 13 angles in a 3D image,
  4 angles in a 2D).

* GLSZM: Use "char" datatype for mask. (It is signed char
  in GLSZM).

* C code is consistent with C89 convention. All variables (pointers
  for python objects) are initialized at top of each block.

Optimizations
-------------

* GLSZM:

  - Uses "while" loops. This allows to reduce the memory usage.
    We observed that with recursive functions it was 'unexpectedly'
    failing.

  - Optimized search that finds a new index to process in
    the region growing.

Co-authored-by: Jean-Christophe Fillion-Robin <[email protected]>

circle.yml

@@ -1,4 +1,3 @@
-
 machine:
   services:
     - docker
@@ -12,12 +11,10 @@ dependencies:
     - pip install scikit-ci-addons==0.11.0
     - ci_addons docker load-pull-save dockcross/manylinux-x64
     - ci_addons docker load-pull-save dockcross/manylinux-x86
-
 test:
   override:
      # Workaround https://github.com/pypa/setuptools/pull/971
      - (unlink venv > /dev/null 2>&1 && echo "symlink removed") || echo "no symlink"
 
      - circleci-matrix:
         parallel: true
-

radiomics/__init__.py

@@ -8,6 +8,7 @@
 import os
 import pkgutil
 import sys
+import traceback
 
 import numpy  # noqa: F401
 
@@ -16,7 +17,6 @@
 if sys.version_info < (2, 6, 0):
   raise ImportError("pyradiomics > 0.9.7 requires python 2.6 or later")
 
-
 def debug(debug_on=True):
   """
   Set up logging system for the whole package.
@@ -72,6 +72,27 @@ class object of the featureClass as value. Assumes only one featureClass per mod
   return _featureClasses
 
 
+def pythonMatrixCalculation(usePython=False):
+  """
+  By default, calculation of matrices 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.
+  Calculation in python can be forced by calling this function, specifying ``pyMatEnabled = True``
+
+  Re-enabling use of C implementation is also done by this function, but if extension is not loaded correctly,
+  a warning is logged and matrix calculation uses python.
+  """
+  global cMatsEnabled
+  if usePython:
+    cMatsEnabled = False
+  elif _cMatLoaded:
+    cMatsEnabled = True
+  else:
+    logger.warning("C Matrices not loaded correctly, cannot calculate matrices in C")
+    cMatsEnabled = False
+
+
 def getInputImageTypes():
   """
   Returns a list of possible input image types. This function finds the image types dynamically by matching the
@@ -99,6 +120,18 @@ def getInputImageTypes():
 logger.addHandler(handler)
 debug(False)  # force level=WARNING, in case logging default is set differently (issue 102)
 
+try:
+  import _cmatrices as cMatrices
+  import _cshape as cShape
+  cMatsEnabled = True
+  _cMatLoaded = True
+except Exception:
+  logger.warning("Error loading C Matrices, switching to python calculation\n%s", traceback.format_exc())
+  cMatrices = None
+  cShape = None
+  cMatsEnabled = False
+  _cMatLoaded = False
+
 _featureClasses = None
 _inputImages = None
 

radiomics/glcm.py

@@ -2,9 +2,9 @@
 from six.moves import range
 from tqdm import trange
 
+import radiomics
 from radiomics import base, imageoperations
 
-
 class RadiomicsGLCM(base.RadiomicsFeaturesBase):
   r"""
   A Gray Level Co-occurrence Matrix (GLCM) of size :math:`N_g \times N_g` describes the second-order joint probability function of an image region
@@ -118,13 +118,16 @@ def __init__(self, inputImage, inputMask, **kwargs):
     self.weightingNorm = kwargs.get('weightingNorm', None)  # manhattan, euclidean, infinity
 
     self.coefficients = {}
-    self.P_glcm = {}
+    self.P_glcm = None
 
     # binning
     self.matrix, self.histogram = imageoperations.binImage(self.binWidth, self.matrix, self.matrixCoordinates)
     self.coefficients['Ng'] = self.histogram[1].shape[0] - 1
 
-    self._calculateGLCM()
+    if radiomics.cMatsEnabled:
+      self.P_glcm = self._calculateCGLCM()
+    else:
+      self.P_glcm = self._calculateGLCM()
     self._calculateCoefficients()
 
   def _calculateGLCM(self):
@@ -142,7 +145,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 +170,12 @@ 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()
 
     # 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 +194,55 @@ 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(P_glcm, (0, 1), keepdims=True)  # , keepdims=True)
+
+    # Delete empty angles if no weighting is applied
+    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
+    return P_glcm / sumP_glcm
+
+  def _calculateCGLCM(self):
+    size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
+    angles = imageoperations.generateAngles(size)
+    Ng = self.coefficients['Ng']
+
+    P_glcm = radiomics.cMatrices.calculate_glcm(self.matrix, self.maskArray, angles, Ng)
+
+    # Optionally make GLCMs symmetrical for each angle
+    if self.symmetricalGLCM:
+      P_glcm += numpy.transpose(P_glcm, (1, 0, 2))
+
+    # Optionally apply a weighting factor
+    if self.weightingNorm is not None:
+      pixelSpacing = self.inputImage.GetSpacing()[::-1]
+      weights = numpy.empty(len(angles))
+      for a_idx, a in enumerate(angles):
+        if self.weightingNorm == 'infinity':
+          weights[a_idx] = numpy.exp(-max(numpy.abs(a) * pixelSpacing) ** 2)
+        elif self.weightingNorm == 'euclidean':
+          weights[a_idx] = numpy.exp(-numpy.sum((numpy.abs(a) * pixelSpacing) ** 2))  # sqrt ^ 2 = 1
+        elif self.weightingNorm == 'manhattan':
+          weights[a_idx] = numpy.exp(-numpy.sum(numpy.abs(a) * pixelSpacing) ** 2)
+        else:
+          self.logger.warning('weigthing norm "%s" is unknown, W is set to 1', self.weightingNorm)
+          weights[a_idx] = 1
+
+      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)  # , 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,6 +3,7 @@
 import numpy
 from six.moves import range
 
+import radiomics
 from radiomics import base, imageoperations
 
 
@@ -93,7 +94,10 @@ def __init__(self, inputImage, inputMask, **kwargs):
     self.coefficients['Nr'] = numpy.max(self.matrix.shape)
     self.coefficients['Np'] = self.targetVoxelArray.size
 
-    self._calculateGLRLM()
+    if radiomics.cMatsEnabled:
+      self.P_glrlm = self._calculateCGLRLM()
+    else:
+      self.P_glrlm = self._calculateGLRLM()
     self._calculateCoefficients()
 
   def _calculateGLRLM(self):
@@ -156,7 +160,7 @@ def _calculateGLRLM(self):
     # 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,16 +179,60 @@ 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 _calculateCGLRLM(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 = radiomics.cMatrices.calculate_glrlm(self.matrix, self.maskArray, angles, Ng, Nr)
+
+    # 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
+    P_glrlm = P_glrlm[xstart:xstop, :ystop, :]
+
+    # Optionally apply a weighting factor
+    if self.weightingNorm is not None:
+      pixelSpacing = self.inputImage.GetSpacing()[::-1]
+      weights = numpy.empty(len(angles))
+      for a_idx, a in enumerate(angles):
+        if self.weightingNorm == 'infinity':
+          weights[a_idx] = max(numpy.abs(a) * pixelSpacing)
+        elif self.weightingNorm == 'euclidean':
+          weights[a_idx] = numpy.sqrt(numpy.sum((numpy.abs(a) * pixelSpacing) ** 2))
+        elif self.weightingNorm == 'manhattan':
+          weights[a_idx] = numpy.sum(numpy.abs(a) * pixelSpacing)
+        else:
+          self.logger.warning('weigthing norm "%s" is unknown, weighting factor is set to 1', self.weightingNorm)
+          weights[a_idx] = 1
+
+      P_glrlm = numpy.sum(P_glrlm * weights[None, None, :], 2, keepdims=True)
+
+    sumP_glrlm = numpy.sum(P_glrlm, (0, 1))
+
+    # Delete empty angles if no weighting is applied
+    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):
 

radiomics/glszm.py

@@ -2,9 +2,9 @@
 from six.moves import range
 from tqdm import trange
 
+import radiomics
 from radiomics import base, imageoperations
 
-
 class RadiomicsGLSZM(base.RadiomicsFeaturesBase):
   r"""
   A Gray Level Size Zone (GLSZM) quantifies gray level zones in an image.
@@ -59,14 +59,17 @@ def __init__(self, inputImage, inputMask, **kwargs):
     super(RadiomicsGLSZM, self).__init__(inputImage, inputMask, **kwargs)
 
     self.coefficients = {}
-    self.P_glszm = {}
+    self.P_glszm = None
 
     # binning
     self.matrix, self.histogram = imageoperations.binImage(self.binWidth, self.matrix, self.matrixCoordinates)
     self.coefficients['Ng'] = self.histogram[1].shape[0] - 1
     self.coefficients['Np'] = self.targetVoxelArray.size
 
-    self._calculateGLSZM()
+    if radiomics.cMatsEnabled:
+      self.P_glszm = self._calculateCGLSZM()
+    else:
+      self.P_glszm = self._calculateGLSZM()
     self._calculateCoefficients()
 
   def _calculateGLSZM(self):
@@ -131,7 +134,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 _calculateCGLSZM(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 radiomics.cMatrices.calculate_glszm(self.matrix, self.maskArray, angles, Ng, Ns)
 
   def _calculateCoefficients(self):
     sumP_glszm = numpy.sum(self.P_glszm, (0, 1))

radiomics/shape.py

@@ -2,6 +2,7 @@
 import SimpleITK as sitk
 from six.moves import range
 
+import radiomics
 from radiomics import base
 
 
@@ -14,7 +15,7 @@ class RadiomicsShape(base.RadiomicsFeaturesBase):
   def __init__(self, inputImage, inputMask, **kwargs):
     super(RadiomicsShape, self).__init__(inputImage, inputMask, **kwargs)
 
-    self.pixelSpacing = inputImage.GetSpacing()[::-1]
+    self.pixelSpacing = numpy.array(inputImage.GetSpacing()[::-1])
     z, x, y = self.pixelSpacing
     self.cubicMMPerVoxel = z * x * y
 
@@ -43,7 +44,10 @@ def __init__(self, inputImage, inputMask, **kwargs):
 
     # Volume and Surface Area are pre-calculated
     self.Volume = self.lssif.GetPhysicalSize(self.label)
-    self.SurfaceArea = self._calculateSurfaceArea()
+    if radiomics.cMatsEnabled:
+      self.SurfaceArea = self._calculateCSurfaceArea()
+    else:
+      self.SurfaceArea = self._calculateSurfaceArea()
 
   def _calculateSurfaceArea(self):
     # define relative locations of the 8 voxels of a sampling cube
@@ -116,6 +120,9 @@ def _calculateSurfaceArea(self):
 
     return S_A
 
+  def _calculateCSurfaceArea(self):
+    return radiomics.cShape.calculate_surfacearea(self.maskArray, self.pixelSpacing)
+
   def _getMaximum2Ddiameter(self, dim):
     otherDims = tuple(set([0, 1, 2]) - set([dim]))