Skip to content

Commit

Permalink
ENH: Implement extension in C for calculation of texture matrices.
Browse files Browse the repository at this point in the history
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)
JoostJM authored and jcfr committed Feb 16, 2017

Verified

This commit was signed with the committer’s verified signature. The key has expired.
jcfr Jean-Christophe Fillion-Robin
1 parent cc78841 commit fc84120
Showing 9 changed files with 903 additions and 26 deletions.
42 changes: 41 additions & 1 deletion radiomics/__init__.py
Original file line number Diff line number Diff line change
@@ -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']
46 changes: 35 additions & 11 deletions radiomics/glcm.py
Original file line number Diff line number Diff line change
@@ -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):
45 changes: 37 additions & 8 deletions radiomics/glrlm.py
Original file line number Diff line number Diff line change
@@ -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)
20 changes: 16 additions & 4 deletions radiomics/glszm.py
Original file line number Diff line number Diff line change
@@ -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))
Loading

0 comments on commit fc84120

Please sign in to comment.