C matrices

.gitignore

@@ -8,6 +8,7 @@ __pycache__/
 
 # Distribution / packaging
 .Python
+_skbuild/
 env/
 build/
 develop-eggs/
@@ -59,3 +60,7 @@ docs/_build/
 
 # PyBuilder
 target/
+
+# scikit-build
+_skbuild
+MANIFEST

CMakeLists.txt

@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 3.7)
+
+project(radiomics)
+
+find_package(PythonInterp REQUIRED)
+find_package(PythonLibs REQUIRED)
+find_package(PythonExtensions REQUIRED)
+
+add_subdirectory(radiomics/src)
+

circle.yml

@@ -1,11 +1,9 @@
 dependencies:
-  pre:
-    - pip install numpy==1.11.0
-    - pip install --trusted-host www.simpleitk.org -f http://www.simpleitk.org/SimpleITK/resources/software.html SimpleITK==0.9.1
-    - pip install nose-parameterized==0.5.0
-    - pip install tqdm==4.7.1
-    - pip install PyWavelets==0.4.0
+  override:
+    - pip install -r requirements.txt
+    - pip install -U cmake scikit-build
 test:
   override:
-    - nosetests --with-xunit --logging-level=DEBUG --verbosity=3 tests/test_features.py
+    - pip install .
+    - python setup.py test  --args="--with-xunit --logging-level=DEBUG"
     - cp nosetests.xml $CIRCLE_TEST_REPORTS

data/schemaFuncs.py

@@ -1,7 +1,7 @@
 import pywt
-from radiomics.featureextractor import RadiomicsFeaturesExtractor
+from radiomics import getFeatureClasses
 
-featureClasses = RadiomicsFeaturesExtractor.getFeatureClasses()
+featureClasses = getFeatureClasses()
 
 
 def checkWavelet(value, rule_obj, path):

radiomics/__init__.py

@@ -1,11 +1,18 @@
 import sys
+import traceback
+import pkgutil
+import inspect
+import os
+import importlib
 
 if sys.version_info < (2, 6, 0):
   raise ImportError("pyradiomics > 0.9.7 requires python 2.6 or later")
 in_py3 = sys.version_info[0] > 2
 
 import logging
 
+from . import base
+
 
 def debug(debug_on=True):
   """
@@ -32,6 +39,59 @@ def debug(debug_on=True):
     debugging = False
 
 
+def getFeatureClasses():
+  """
+  Iterates over all modules of the radiomics package using pkgutil and subsequently imports those modules.
+
+  Return a dictionary of all modules containing featureClasses, with modulename as key, abstract
+  class object of the featureClass as value. Assumes only one featureClass per module
+
+  This is achieved by inspect.getmembers. Modules are added if it contains a member that is a class,
+  with name starting with 'Radiomics' and is inherited from :py:class:`radiomics.base.RadiomicsFeaturesBase`.
+
+  This iteration only runs once, subsequent calls return the dictionary created by the first call.
+  """
+  global _featureClasses
+  if _featureClasses is not None:
+    return _featureClasses
+
+  featureClasses = {}
+  for _, mod, _ in pkgutil.iter_modules([os.path.dirname(base.__file__)]):
+    if str(mod).startswith('_'):  # Do not load _version, _cmatrices and _cshape here
+      continue
+    __import__('radiomics.' + mod)
+    module = sys.modules['radiomics.' + mod]
+    attributes = inspect.getmembers(module, inspect.isclass)
+    for a in attributes:
+      if a[0].startswith('Radiomics'):
+        if base.RadiomicsFeaturesBase in inspect.getmro(a[1])[1:]:
+          featureClasses[mod] = a[1]
+  _featureClasses = featureClasses
+
+  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
+
+
 debugging = True
 logger = logging.getLogger(__name__)
 handler = logging.StreamHandler()
@@ -41,10 +101,26 @@ def debug(debug_on=True):
 logger.addHandler(handler)
 debug(False)  # force level=WARNING, in case logging default is set differently (issue 102)
 
+_featureClasses = None
+
+try:
+  from . import _cmatrices as cMatrices
+  from . 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
+
 # For convenience, import the most used packages into the "pyradiomics" namespace
 import collections, numpy
 
 from ._version import get_versions
 
 __version__ = get_versions()['version']
 del get_versions
+
+getFeatureClasses()

radiomics/base.py

@@ -1,9 +1,9 @@
+import inspect
 import logging
 import traceback
-import SimpleITK as sitk
+
 import numpy
-import inspect
-from radiomics import imageoperations
+import SimpleITK as sitk
 
 
 class RadiomicsFeaturesBase(object):

radiomics/featureextractor.py

@@ -1,15 +1,15 @@
 # -*- coding: utf-8 -*-
-import os
-import logging
 import collections
-from itertools import chain
+import logging
+import os
+
 import numpy
-import SimpleITK as sitk
-import pkgutil
-import inspect
 import pykwalify.core
-import radiomics
-from radiomics import base, imageoperations, generalinfo
+import SimpleITK as sitk
+
+from itertools import chain
+
+from . import imageoperations, generalinfo, getFeatureClasses
 
 
 class RadiomicsFeaturesExtractor:
@@ -64,7 +64,7 @@ class RadiomicsFeaturesExtractor:
   def __init__(self, *args, **kwargs):
     self.logger = logging.getLogger(__name__)
 
-    self.featureClasses = self.getFeatureClasses()
+    self.featureClasses = getFeatureClasses()
 
     self.kwargs = {}
     self.provenance_on = True
@@ -136,7 +136,7 @@ def loadParams(self, paramsFile):
 
     If supplied params file does not match the requirements, a pykwalify error is raised.
     """
-    dataDir = os.path.abspath(os.path.join(radiomics.__path__[0], '..', 'data'))
+    dataDir = os.path.abspath(os.path.join(self.__path__[0], '..', 'data'))
     schemaFile = os.path.join(dataDir, 'paramSchema.yaml')
     schemaFuncs = os.path.join(dataDir, 'schemaFuncs.py')
     c = pykwalify.core.Core(source_file=paramsFile, schema_files=[schemaFile], extensions=[schemaFuncs])
@@ -561,28 +561,6 @@ def getInputImageTypes(self):
     """
     return [member[9:] for member in dir(self) if member.startswith('generate_')]
 
-  @classmethod
-  def getFeatureClasses(cls):
-    """
-    Iterates over all modules of the radiomics package using pkgutil and subsequently imports those modules.
-
-    Return a dictionary of all modules containing featureClasses, with modulename as key, abstract
-    class object of the featureClass as value. Assumes only one featureClass per module
-
-    This is achieved by inspect.getmembers. Modules are added if it contains a memeber that is a class,
-    with name starting with 'Radiomics' and is inherited from :py:class:`radiomics.base.RadiomicsFeaturesBase`.
-    """
-    featureClasses = {}
-    for _, mod, _ in pkgutil.iter_modules(radiomics.__path__):
-      __import__('radiomics.' + mod)
-      attributes = inspect.getmembers(eval('radiomics.' + mod), inspect.isclass)
-      for a in attributes:
-        if a[0].startswith('Radiomics'):
-          if radiomics.base.RadiomicsFeaturesBase in inspect.getmro(a[1])[1:]:
-            featureClasses[mod] = a[1]
-
-    return featureClasses
-
   def getFeatureClassNames(self):
     return self.featureClasses.keys()
 

radiomics/firstorder.py

@@ -1,7 +1,6 @@
 import numpy
-import collections
-from radiomics import base, imageoperations
-import SimpleITK as sitk
+
+from . import base, imageoperations
 
 
 class RadiomicsFirstOrder(base.RadiomicsFeaturesBase):

radiomics/generalinfo.py

@@ -1,8 +1,9 @@
 import collections
 import logging
-import numpy
+
 import SimpleITK as sitk
-import radiomics
+
+from . import __version__
 
 
 class GeneralInfo():
@@ -115,7 +116,7 @@ def getVersionValue(self):
     """
     Return the current version of this package.
     """
-    return radiomics.__version__
+    return __version__
 
   def getVolumeNumValue(self):
     """

radiomics/glcm.py

@@ -1,9 +1,8 @@
 import numpy
-import collections
-from radiomics import base, imageoperations
-import SimpleITK as sitk
+
 from tqdm import trange
 
+from . import base, cMatrices, cMatsEnabled, imageoperations
 
 class RadiomicsGLCM(base.RadiomicsFeaturesBase):
   r"""
@@ -116,13 +115,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 cMatsEnabled:
+      self.P_glcm = self._calculateCGLCM()
+    else:
+      self.P_glcm = self._calculateGLCM()
     self._calculateCoefficients()
 
   def _calculateGLCM(self):
@@ -140,7 +142,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')
 
@@ -165,12 +167,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 not self.weightingNorm is None:
@@ -189,17 +191,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)
+
+    sumGlcm = 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(sumGlcm == 0), 2)
+      sumGlcm = numpy.delete(sumGlcm, numpy.where(sumGlcm == 0), 0)
+
+    # Normalize each glcm
+    return P_glcm/sumGlcm
+
+  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 = 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 not self.weightingNorm is 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)
 
-    sumGlcm = numpy.sum(self.P_glcm, (0, 1), keepdims=True)  # , keepdims=True)
+    sumGlcm = 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(sumGlcm == 0), 2)
+    if P_glcm.shape[2] > 1:
+      P_glcm = numpy.delete(P_glcm, numpy.where(sumGlcm == 0), 2)
       sumGlcm = numpy.delete(sumGlcm, numpy.where(sumGlcm == 0), 0)
 
     # Normalize each glcm
-    self.P_glcm = self.P_glcm / sumGlcm
+    return P_glcm / sumGlcm
 
   # check if ivector and jvector can be replaced
   def _calculateCoefficients(self):