From 6461197a4e9270aefaa8e93182aac3fe0f28ed16 Mon Sep 17 00:00:00 2001 From: Joost van Griethuysen Date: Sun, 13 Nov 2016 20:54:05 -0500 Subject: [PATCH] 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 --- circle.yml | 3 - radiomics/__init__.py | 35 +++- radiomics/glcm.py | 63 +++++- radiomics/glrlm.py | 60 +++++- radiomics/glszm.py | 19 +- radiomics/shape.py | 11 +- radiomics/src/_cmatrices.c | 346 +++++++++++++++++++++++++++++++ radiomics/src/_cshape.c | 91 +++++++++ radiomics/src/cmatrices.c | 404 +++++++++++++++++++++++++++++++++++++ radiomics/src/cmatrices.h | 5 + radiomics/src/cshape.c | 239 ++++++++++++++++++++++ radiomics/src/cshape.h | 2 + setup.py | 15 +- tests/test_cmatrices.py | 77 +++++++ 14 files changed, 1342 insertions(+), 28 deletions(-) create mode 100644 radiomics/src/_cmatrices.c create mode 100644 radiomics/src/_cshape.c create mode 100644 radiomics/src/cmatrices.c create mode 100644 radiomics/src/cmatrices.h create mode 100644 radiomics/src/cshape.c create mode 100644 radiomics/src/cshape.h create mode 100644 tests/test_cmatrices.py diff --git a/circle.yml b/circle.yml index f798146d..a2046f07 100644 --- a/circle.yml +++ b/circle.yml @@ -1,4 +1,3 @@ - machine: services: - docker @@ -12,7 +11,6 @@ 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 @@ -20,4 +18,3 @@ test: - circleci-matrix: parallel: true - diff --git a/radiomics/__init__.py b/radiomics/__init__.py index b1e15617..bd67a663 100644 --- a/radiomics/__init__.py +++ b/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 diff --git a/radiomics/glcm.py b/radiomics/glcm.py index 07ca6312..89a1bd93 100644 --- a/radiomics/glcm.py +++ b/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): diff --git a/radiomics/glrlm.py b/radiomics/glrlm.py index d12e97cb..14c71a90 100644 --- a/radiomics/glrlm.py +++ b/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): diff --git a/radiomics/glszm.py b/radiomics/glszm.py index b1132c78..5b1d5e6c 100644 --- a/radiomics/glszm.py +++ b/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)) diff --git a/radiomics/shape.py b/radiomics/shape.py index 690dba95..b4800a5b 100644 --- a/radiomics/shape.py +++ b/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])) diff --git a/radiomics/src/_cmatrices.c b/radiomics/src/_cmatrices.c new file mode 100644 index 00000000..868bc454 --- /dev/null +++ b/radiomics/src/_cmatrices.c @@ -0,0 +1,346 @@ +#include +#include +#include "cmatrices.h" + +static char module_docstring[] = ("This module links to C-compiled code for efficient calculation of various matrices " + "in the pyRadiomics package. It provides fast calculation for GLCM, GLDM, NGTDM, " + "GLRLM and GLSZM. All functions are given names as follows: ""calculate_"", " + "where is the name of the matix, in lowercase. Arguments for these functions " + "are positional and start with 3 numpy arrays (image, mask and angles) and 1 integer " + "(Ng, number of gray levels). Optionally extra arguments may be required, see function " + "docstrings for detailed information."); +static char glcm_docstring[] = "Arguments: Image, Mask, Angles, Ng."; +static char glszm_docstring[] = "Arguments: Image, Mask, Angles, Ng, Ns, matrix is cropped to maximum size encountered."; +static char glrlm_docstring[] = "Arguments: Image, Mask, Angles, Ng, Nr."; + +static PyObject *cmatrices_calculate_glcm(PyObject *self, PyObject *args); +static PyObject *cmatrices_calculate_glszm(PyObject *self, PyObject *args); +static PyObject *cmatrices_calculate_glrlm(PyObject *self, PyObject *args); + +static PyMethodDef module_methods[] = { + //{"calculate_", cmatrices_, METH_VARARGS, _docstring}, + { "calculate_glcm", cmatrices_calculate_glcm, METH_VARARGS, glcm_docstring }, + { "calculate_glszm", cmatrices_calculate_glszm, METH_VARARGS, glszm_docstring }, + { "calculate_glrlm", cmatrices_calculate_glrlm, METH_VARARGS, glrlm_docstring }, + { NULL, NULL, 0, NULL } +}; + +PyMODINIT_FUNC +init_cmatrices(void) { + PyObject *m = Py_InitModule3("_cmatrices", module_methods, module_docstring); + if (m == NULL) + return; + + // Initialize numpy functionality + import_array(); +} + +static PyObject *cmatrices_calculate_glcm(PyObject *self, PyObject *args) +{ + int Ng; + PyObject *image_obj, *mask_obj, *angles_obj; + PyArrayObject *image_arr, *mask_arr, *angles_arr; + int Sx, Sy, Sz, Na; + int dims[3]; + PyArrayObject *glcm_arr; + int *image; + char *mask; + int *angles; + double *glcm; + int glcm_size, k; + + // Parse the input tuple + if (!PyArg_ParseTuple(args, "OOOi", &image_obj, &mask_obj, &angles_obj, &Ng)) + return NULL; + + // Interpret the input as numpy arrays + image_arr = (PyArrayObject *)PyArray_FROM_OTF(image_obj, NPY_INT, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + mask_arr = (PyArrayObject *)PyArray_FROM_OTF(mask_obj, NPY_BOOL, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + angles_arr = (PyArrayObject *)PyArray_FROM_OTF(angles_obj, NPY_INT, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + + if (image_arr == NULL || mask_arr == NULL || angles_arr == NULL) + { + Py_XDECREF(image_arr); + Py_XDECREF(mask_arr); + Py_XDECREF(angles_arr); + return NULL; + } + + // Check if Image and Mask have 3 dimensions + if (PyArray_NDIM(image_arr) != 3 || PyArray_NDIM(mask_arr) != 3) + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + PyErr_SetString(PyExc_RuntimeError, "Expected a 3D array for image and mask."); + return NULL; + } + + // Get sizes of the arrays + Sx = (int)PyArray_DIM(image_arr, 2); + Sy = (int)PyArray_DIM(image_arr, 1); + Sz = (int)PyArray_DIM(image_arr, 0); + + Na = (int)PyArray_DIM(angles_arr, 0); + + // Check if image and mask are the same size + if (Sx != (int)PyArray_DIM(mask_arr, 2) || Sy != (int)PyArray_DIM(mask_arr, 1) || Sz != (int)PyArray_DIM(mask_arr, 0)) + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + PyErr_SetString(PyExc_RuntimeError, "Dimensions of image and mask do not match."); + return NULL; + } + + // Initialize output array (elements not set) + dims[0] = Ng; + dims[1] = Ng; + dims[2] = Na; + glcm_arr = (PyArrayObject *)PyArray_FromDims(3, (int*)dims, NPY_DOUBLE); + + // Get arrays in Ctype + image = (int *)PyArray_DATA(image_arr); + mask = (char *)PyArray_DATA(mask_arr); + angles = (int *)PyArray_DATA(angles_arr); + glcm = (double *)PyArray_DATA(glcm_arr); + + // Set all elements to 0 + glcm_size = Ng * Ng * Na; + for (k = 0; k < glcm_size; k++) glcm[k] = 0; + + //Calculate GLCM + if (!calculate_glcm(image, mask, Sx, Sy, Sz, angles, Na, glcm, Ng)) + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + Py_DECREF(glcm_arr); + PyErr_SetString(PyExc_RuntimeError, "Calculation of GLCM Failed."); + return NULL; + } + + // Clean up + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + + return PyArray_Return(glcm_arr); +} + +static PyObject *cmatrices_calculate_glszm(PyObject *self, PyObject *args) +{ + int Ng, Ns; + PyObject *image_obj, *mask_obj, *angles_obj; + PyArrayObject *image_arr, *mask_arr, *angles_arr; + int Sx, Sy, Sz, Na; + int tempDims[2]; + PyArrayObject *temp_arr; + int *image; + signed char *mask; + int *angles; + int * tempData; + int tempData_size, k; + int maxRegion; + int dims[2]; + PyArrayObject *glszm_arr; + double *glszm; + + // Parse the input tuple + if (!PyArg_ParseTuple(args, "OOOii", &image_obj, &mask_obj, &angles_obj, &Ng, &Ns)) + return NULL; + + // Interpret the input as numpy arrays + image_arr = (PyArrayObject *)PyArray_FROM_OTF(image_obj, NPY_INT, NPY_FORCECAST | NPY_UPDATEIFCOPY |NPY_IN_ARRAY); + mask_arr = (PyArrayObject *)PyArray_FROM_OTF(mask_obj, NPY_BYTE, NPY_FORCECAST | NPY_ENSURECOPY | NPY_IN_ARRAY); + angles_arr = (PyArrayObject *)PyArray_FROM_OTF(angles_obj, NPY_INT, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + + if (image_arr == NULL || mask_arr == NULL || angles_arr == NULL) + { + Py_XDECREF(image_arr); + Py_XDECREF(mask_arr); + Py_XDECREF(angles_arr); + return NULL; + } + + // Check if Image and Mask have 3 dimensions + if (PyArray_NDIM(image_arr) != 3 || PyArray_NDIM(mask_arr) != 3) + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + PyErr_SetString(PyExc_RuntimeError, "Expected a 3D array for image and mask."); + return NULL; + } + + // Get sizes of the arrays + Sx = (int)PyArray_DIM(image_arr, 2); + Sy = (int)PyArray_DIM(image_arr, 1); + Sz = (int)PyArray_DIM(image_arr, 0); + + Na = (int)PyArray_DIM(angles_arr, 0); + + // Check if image and mask are the same size + if (Sx != (int)PyArray_DIM(mask_arr, 2) || Sy != (int)PyArray_DIM(mask_arr, 1) || Sz != (int)PyArray_DIM(mask_arr, 0)) + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + PyErr_SetString(PyExc_RuntimeError, "Dimensions of image and mask do not match."); + return NULL; + } + + // Initialize temporary output array (elements not set) + // add +1 to the first dimension so in the case every voxel represents a separate region, + // tempData still contains a -1 element at the end + tempDims[0] = Ns + 1; + tempDims[1] = 2; + temp_arr = (PyArrayObject *)PyArray_FromDims(2, (int*)tempDims, NPY_INT); + + // Get arrays in Ctype + image = (int *)PyArray_DATA(image_arr); + mask = (signed char *)PyArray_DATA(mask_arr); + angles = (int *)PyArray_DATA(angles_arr); + tempData = (int *)PyArray_DATA(temp_arr); + + // Set all elements to 0 + tempData_size = tempDims[0] * tempDims[1]; + for (k = 0; k < tempData_size; k++) tempData[k] = -1; + + //Calculate GLSZM + maxRegion = 0; + maxRegion = calculate_glszm(image, mask, Sx, Sy, Sz, angles, Na, tempData, Ng, Ns); + if (maxRegion == -1) // Error occured + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + Py_DECREF(temp_arr); + PyErr_SetString(PyExc_RuntimeError, "Calculation of GLSZM Failed."); + return NULL; + } + + // Clean up image, mask and angles arrays (not needed anymore) + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + + // Initialize temporary output array (elements not set) + if (maxRegion == 0) maxRegion = 1; + dims[0] = Ng; + dims[1] = maxRegion; + glszm_arr = (PyArrayObject *)PyArray_FromDims(2, (int*)dims, NPY_DOUBLE); + glszm = (double *)PyArray_DATA(glszm_arr); + + if (!fill_glszm(tempData, glszm, Ng, maxRegion)) + { + Py_DECREF(temp_arr); + Py_DECREF(glszm_arr); + PyErr_SetString(PyExc_RuntimeError, "Error filling GLSZM."); + return NULL; + } + + // Clean up tempData + Py_DECREF(temp_arr); + + return PyArray_Return(glszm_arr); +} + +static PyObject *cmatrices_calculate_glrlm(PyObject *self, PyObject *args) +{ + int Ng, Nr; + PyObject *image_obj, *mask_obj, *angles_obj; + PyArrayObject *image_arr, *mask_arr, *angles_arr; + int Na; + int size[3]; + int strides[3]; + int dims[3]; + PyArrayObject *glrlm_arr; + int *image; + char *mask; + int *angles; + double *glrlm; + int glrlm_size, k; + + // Parse the input tuple + if (!PyArg_ParseTuple(args, "OOOii", &image_obj, &mask_obj, &angles_obj, &Ng, &Nr)) + return NULL; + + // Interpret the input as numpy arrays + image_arr = (PyArrayObject *)PyArray_FROM_OTF(image_obj, NPY_INT, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + mask_arr = (PyArrayObject *)PyArray_FROM_OTF(mask_obj, NPY_BOOL, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + angles_arr = (PyArrayObject *)PyArray_FROM_OTF(angles_obj, NPY_INT, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + + if (image_arr == NULL || mask_arr == NULL || angles_arr == NULL) + { + Py_XDECREF(image_arr); + Py_XDECREF(mask_arr); + Py_XDECREF(angles_arr); + return NULL; + } + + // Check if Image and Mask have 3 dimensions + if (PyArray_NDIM(image_arr) != 3 || PyArray_NDIM(mask_arr) != 3) + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + PyErr_SetString(PyExc_RuntimeError, "Expected a 3D array for image and mask."); + return NULL; + } + + // Get sizes of the arrays + size[2] = (int)PyArray_DIM(image_arr, 2); + size[1] = (int)PyArray_DIM(image_arr, 1); + size[0] = (int)PyArray_DIM(image_arr, 0); + + strides[2] = (int)(image_arr->strides[2] / image_arr->descr->elsize); + strides[1] = (int)(image_arr->strides[1] / image_arr->descr->elsize); + strides[0] = (int)(image_arr->strides[0] / image_arr->descr->elsize); + + Na = (int)PyArray_DIM(angles_arr, 0); + + // Check if image and mask are the same size + if (size[2] != (int)PyArray_DIM(mask_arr, 2) || size[1] != (int)PyArray_DIM(mask_arr, 1) || size[0] != (int)PyArray_DIM(mask_arr, 0)) + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + PyErr_SetString(PyExc_RuntimeError, "Dimensions of image and mask do not match."); + return NULL; + } + + // Initialize output array (elements not set) + dims[0] = Ng; + dims[1] = Nr; + dims[2] = Na; + glrlm_arr = (PyArrayObject *)PyArray_FromDims(3, (int*)dims, NPY_DOUBLE); + + // Get arrays in Ctype + image = (int *)PyArray_DATA(image_arr); + mask = (char *)PyArray_DATA(mask_arr); + angles = (int *)PyArray_DATA(angles_arr); + glrlm = (double *)PyArray_DATA(glrlm_arr); + + // Set all elements to 0 + glrlm_size = Ng * Nr * Na; + for (k = 0; k < glrlm_size; k++) glrlm[k] = 0; + + //Calculate GLRLM + if (!calculate_glrlm(image, mask, size, strides, angles, Na, glrlm, Ng, Nr)) + { + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + Py_DECREF(glrlm_arr); + PyErr_SetString(PyExc_RuntimeError, "Calculation of GLRLM Failed."); + return NULL; + } + + // Clean up + Py_DECREF(image_arr); + Py_DECREF(mask_arr); + Py_DECREF(angles_arr); + + return PyArray_Return(glrlm_arr); +} diff --git a/radiomics/src/_cshape.c b/radiomics/src/_cshape.c new file mode 100644 index 00000000..9ec739e5 --- /dev/null +++ b/radiomics/src/_cshape.c @@ -0,0 +1,91 @@ +#include +#include +#include "cshape.h" + +static char module_docstring[] = "This module links to C-compiled code for efficient calculation of the surface area " + "in the pyRadiomics package. It provides fast calculation using a marching cubes " + "algortihm, accessed via ""calculate_surfacearea"". Arguments for this function" + "are positional and consist of two numpy arrays, mask and pixelspacing, which must " + "be supplied in this order. Pixelspacing is a 3 element vector containing the pixel" + "spacing in z, y and x dimension, respectively. All non-zero elements in mask are " + "considered to be a part of the segmentation and are included in the algorithm."; +static char surface_docstring[] = "Arguments: Mask, PixelSpacing, uses a marching cubes algorithm to calculate an " + "approximation to the total surface area. The isovalue is considered to be situated " + "midway between a voxel that is part of the segmentation and a voxel that is not."; + +static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args); + +static PyMethodDef module_methods[] = { + //{"calculate_", cmatrices_, METH_VARARGS, _docstring}, + { "calculate_surfacearea", cshape_calculate_surfacearea, METH_VARARGS, surface_docstring }, + { NULL, NULL, 0, NULL } +}; + +PyMODINIT_FUNC +init_cshape(void) { + PyObject *m = Py_InitModule3("_cshape", module_methods, module_docstring); + if (m == NULL) + return; + + // Initialize numpy functionality + import_array(); +} + +static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args) +{ + PyObject *mask_obj, *spacing_obj; + PyArrayObject *mask_arr, *spacing_arr; + int size[3]; + char *mask; + double *spacing; + double SA; + + // Parse the input tuple + if (!PyArg_ParseTuple(args, "OO", &mask_obj, &spacing_obj)) + return NULL; + + // Interpret the input as numpy arrays + mask_arr = (PyArrayObject *)PyArray_FROM_OTF(mask_obj, NPY_BYTE, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + spacing_arr = (PyArrayObject *)PyArray_FROM_OTF(spacing_obj, NPY_DOUBLE, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + + // Check if there were any errors during casting to array objects + if (mask_arr == NULL || spacing_arr == NULL) + { + Py_XDECREF(mask_arr); + Py_XDECREF(spacing_arr); + return NULL; + } + + // Check if Image and Mask have 3 dimensions + if (PyArray_NDIM(mask_arr) != 3) + { + Py_DECREF(mask_arr); + Py_DECREF(spacing_arr); + PyErr_SetString(PyExc_RuntimeError, "Expected a 3D array for mask."); + return NULL; + } + + // Get sizes of the arrays + size[2] = (int)PyArray_DIM(mask_arr, 2); + size[1] = (int)PyArray_DIM(mask_arr, 1); + size[0] = (int)PyArray_DIM(mask_arr, 0); + + // Get arrays in Ctype + mask = (char *)PyArray_DATA(mask_arr); + spacing = (double *)PyArray_DATA(spacing_arr); + + //Calculate GLCM + SA = calculate_surfacearea(mask, size, spacing); + + // Clean up + Py_DECREF(mask_arr); + Py_DECREF(spacing_arr); + + if (SA < 0) // if SA < 0, an error has occurred + { + PyErr_SetString(PyExc_RuntimeError, "Calculation of GLCM Failed."); + return NULL; + } + + return Py_BuildValue("f", SA); +} diff --git a/radiomics/src/cmatrices.c b/radiomics/src/cmatrices.c new file mode 100644 index 00000000..5ba84a1a --- /dev/null +++ b/radiomics/src/cmatrices.c @@ -0,0 +1,404 @@ +#include "cmatrices.h" + +int calculate_glcm(int *image, char *mask, int Sx, int Sy, int Sz, int *angles, int Na, double *glcm, int Ng) +{ + /* Calculate GLCM: Count the number of voxels with gray level i is neighboured by a voxel with gray level j in + * direction and distance specified by a. Returns an asymmetrical GLCM matrix for each angle/distance + * defined in angles. + */ + int glcm_idx_max = Ng * Ng * Na; + int i = 0, j = 0; + int iz, iy, ix; + int a, glcm_idx; + for (iz = 0; iz < Sz; iz ++) // Iterate over all voxels in a row - column - slice order + { + for (iy = 0; iy < Sy; iy ++) + { + for (ix = 0; ix < Sx; ix ++) + { + if (mask[i]) // Check if the current voxel is part of the segmentation + { + for (a = 0; a < Na; a++) // Iterate over angles to get the neighbours + { + // Check whether the neighbour index is not out of range (i.e. part of the image) + if (iz + angles[a * 3] >= 0 && iz + angles[a * 3] < Sz && + iy + angles[a * 3 + 1] >= 0 && iy + angles[a * 3 + 1] < Sy && + ix + angles[a * 3 + 2] >= 0 && ix + angles[a * 3 + 2] < Sx) + { + j = i + angles[a * 3 + 2] + angles[a * 3 + 1] * Sx + angles[a * 3] * Sy * Sx; + if (mask[j]) // Check whether neighbour voxel is part of the segmentation + { + glcm_idx = a + (image[j]-1) * Na + (image[i]-1) * Na * Ng; + if (glcm_idx >= glcm_idx_max) return 0; // Index out of range + glcm[glcm_idx] ++; + } + } + } + } + // Increase index to point to next item. Only one index is needed as data is stored one dimensionally in memory. + // This is in row-column-slice order. + i++; + } + } + } + return 1; +} + +int calculate_glszm(int *image, signed char *mask, int Sx, int Sy, int Sz, int *angles, int Na, int *tempData, int Ng, int Ns) +{ + int Ni = Sz * Sy * Sx; + int maxRegion = 0; + + int regionCounter = 0; + int max_region_idx = Ns * 2; + + int i = 0, j = 0, cur_idx = 0; + int iz, iy, ix; + int gl, region; + int cur_z, cur_y, cur_x, ref_idx; + int d, a; + + for (iz = 0; iz < Sz; iz++) + { + for (iy = 0; iy < Sy; iy++) + { + for (ix = 0; ix < Sx; ix++) + { + // Check if the current voxel is part of the segmentation and unprocessed + if (mask[i] > 0) + { + // Store current gray level, needed for region growing and determines index in GLSZM. + gl = image[i]; + + // Instantiate variable to hold region size at 1. Function region increases size for every found neighbour. + // The initial size of 1 signifies current voxel + region = 0; + + // Start growing the region + cur_idx = i; + + while (cur_idx > -1) + { + // Exclude current voxel to prevent reprocessing -> 0 is processed or not part of segmentation + // -1 is marked for further processing + mask[cur_idx] = 0; + // Increment region size, as number of loops corresponds to number of voxels in current region + region++; + + // Generate neighbours for current voxel + cur_z = (cur_idx / (Sx * Sy)); + cur_y = (cur_idx % (Sx * Sy)) / Sx; + cur_x = (cur_idx % (Sx * Sy)) % Sx; + for (d = 1; d >= -1; d -= 2) // Iterate over both directions for each angle + { + for (a = 0; a < Na; a++) // Iterate over angles to get the neighbours + { + // Check whether the neighbour index is not out of range (i.e. part of the image) + if (cur_z + d * angles[a * 3] >= 0 && cur_z + d * angles[a * 3] < Sz && + cur_y + d * angles[a * 3 + 1] >= 0 && cur_y + d * angles[a * 3 + 1] < Sy && + cur_x + d * angles[a * 3 + 2] >= 0 && cur_x + d * angles[a * 3 + 2] < Sx) + { + j = cur_idx + d * angles[a * 3 + 2] + + d * angles[a * 3 + 1] * Sx + + d * angles[a * 3] * Sy * Sx; + // Check whether neighbour voxel is part of the segmentation and unprocessed + if (mask[j] && (image[j] == gl)) + { + // Voxel belongs to current region, mark it for further processing + mask[j] = -1; + } + } + } // next a + } // next d + + if (mask[j] == -1) + { + cur_idx = j; + } + else + { + ref_idx = cur_idx; + cur_idx ++; // cur_idx doesn't have to be checked + // try to find the next voxel that has been marked for processing (-1) + // If none are found, current region is complete + // increment cur_idx until end of image or a voxel that has been marked is found + while (cur_idx < Ni && !(mask[cur_idx] == -1)) cur_idx++; + if (cur_idx == Ni) cur_idx = ref_idx; // no voxel found, check between ref_idx and i + { + while (cur_idx > i && !(mask[cur_idx] == -1)) cur_idx--; + if (cur_idx == i) cur_idx = -1; + } + } + } // while region + + if (regionCounter >= max_region_idx) return -1; // index out of range + + if (region > maxRegion) maxRegion = region; + + tempData[(regionCounter * 2)] = gl; + tempData[((regionCounter * 2) + 1)] = region; + + regionCounter ++; + } + // Increase index to point to next item. Only one index is needed as data is stored one dimensionally in memory. + // This is in row-column-slice order. + i++; + } + } + } +return maxRegion; +} + +int fill_glszm(int *tempData, double *glszm, int Ng, int maxRegion) +{ + int i = 0; + int glszm_idx_max = Ng * maxRegion; + int glszm_idx; + + while(tempData[i * 2] > -1) + { + glszm_idx = (tempData[i * 2] - 1) * maxRegion + tempData[i * 2 + 1] - 1; + if (glszm_idx >= glszm_idx_max) return 0; // Index out of range + + glszm[glszm_idx]++; + i++; + } + return 1; +} + +int calculate_glrlm(int *image, char *mask, int *size, int *strides, int *angles, int Na, double *glrlm, int Ng, int Nr) +{ + int a, da, nd; + int mDims[3], sDims[2]; + int d1, d2; + int jd[3]; + int glrlm_idx_max = Ng * Nr * Na; + int runVal; + char multiElement; + + for (a = 0; a < Na; a++) // Iterate over angles to get the neighbours + { + multiElement = 0; + nd = 0; + for (da = 0; da < 3; da++) + { + if (angles[a * 3 + da] == 0) + { + sDims[da - nd] = da; + } + else + { + mDims[nd] = da; + nd++; + } + } + if (nd == 1) // Moving in 1 dimension + { + // Iterate over all start voxels: All voxels, where moving dim == 0 or max (i.e. all combinations of voxels + // in static dimensions) + for (d1 = 0; d1 < size[sDims[0]]; d1++) + { + for (d2 = 0; d2 < size[sDims[1]]; d2++) + { + jd[sDims[0]] = d1; + jd[sDims[1]] = d2; + if (angles[a * 3 + mDims[0]]< 0) + { + jd[mDims[0]] = size[mDims[0]] - 1; // Moving dim is negative, set start at maximum index + } + else + { + jd[mDims[0]] = 0; // Moving dim is positive, set start at 0 + } + // run_diagonal returns -1 if error, otherwise number of elements encountered + runVal = run_diagonal(image, mask, size, strides, angles, Na, glrlm, glrlm_idx_max, Nr, jd, a); + if (runVal < 0) return runVal; // index out of range, 'raise' error. + if (!multiElement && runVal > 1) multiElement = 1; // multiple elements found + } + } + } + else if (nd == 2) // Moving in 2 dimensions + { + for (d1 = 0; d1 < size[sDims[0]]; d1++) // iterate over all voxels in the static dimension + { + for (d2 = 0; d2 < size[mDims[0]]; d2++) // Iterate over moving dimension 1, treating it as a 'static' + { + jd[sDims[0]] = d1; + jd[mDims[0]] = d2; + if (angles[a * 3 + mDims[1]] < 0) + { + jd[mDims[1]] = size[mDims[1]] - 1; // Moving dim is negative, set start at maximum index + } + else + { + jd[mDims[1]] = 0; // Moving dim is positive, set start at 0 + } + // run_diagonal returns -1 if error, otherwise number of elements encountered + runVal = run_diagonal(image, mask, size, strides, angles, Na, glrlm, glrlm_idx_max, Nr, jd, a); + if (runVal < 0) return runVal; // index out of range, 'raise' error. + if (!multiElement && runVal > 1) multiElement = 1; // multiple elements found + } + for (d2 = 1; d2 < size[mDims[1]]; d2++) // Iterate over other moving dimension, treating it as a 'static' + { + jd[sDims[0]] = d1; + jd[mDims[1]] = d2; + // Prevent calculating the true diagonal twice, which is either at index 0, or max, depending on the + // angle of the the second moving dimension: iterate 1 to max if angle is positive, 0 to max -1 if + // negative. The ignored index is handled in the previous loop. + if (angles[a * 3 + mDims[1]] < 0) jd[mDims[1]] --; + if (angles[a * 3 + mDims[0]] < 0) + { + jd[mDims[0]] = size[mDims[0]] - 1; // Moving dim is negative, set start at maximum index + } + else + { + jd[mDims[0]] = 0; // Moving dim is positive, set start at 0 + } + // run_diagonal returns -1 if error, otherwise number of elements encountered + runVal = run_diagonal(image, mask, size, strides, angles, Na, glrlm, glrlm_idx_max, Nr, jd, a); + if (runVal < 0) return runVal; // index out of range, 'raise' error. + if (!multiElement && runVal > 1) multiElement = 1; // multiple elements found + } + } + } + else if (nd == 3) + { + for (d1 = 0; d1 < size[mDims[0]]; d1++) + { + for (d2 = 0; d2 < size[mDims[1]]; d2++) + { + jd[mDims[0]] = d1; + jd[mDims[1]] = d2; + if (angles[a * 3 + mDims[2]] < 0) + { + jd[mDims[2]] = size[mDims[2]] - 1; // Moving dim is negative, set start at maximum index + } + else + { + jd[mDims[2]] = 0; // Moving dim is positive, set start at 0 + } + // run_diagonal returns -1 if error, otherwise number of elements encountered + runVal = run_diagonal(image, mask, size, strides, angles, Na, glrlm, glrlm_idx_max, Nr, jd, a); + if (runVal < 0) return runVal; // index out of range, 'raise' error. + if (!multiElement && runVal > 1) multiElement = 1; // multiple elements found + } + } + for (d1 = 0; d1 < size[mDims[1]]; d1++) + { + for (d2 = 1; d2 < size[mDims[2]]; d2++) + { + jd[mDims[1]] = d1; + jd[mDims[2]] = d2; + // Prevent calculating the true diagonal twice, which is either at index 0, or max, depending on the + // angle of the the second moving dimension: iterate 1 to max if angle is positive, 0 to max -1 if + // negative. The ignored index is handled in the previous loop. + if (angles[a * 3 + mDims[2]] < 0) jd[mDims[2]] --; + if (angles[a * 3 + mDims[0]] < 0) + { + jd[mDims[0]] = size[mDims[0]] - 1; // Moving dim is negative, set start at maximum index + } + else + { + jd[mDims[0]] = 0; // Moving dim is positive, set start at 0 + } + // run_diagonal returns -1 if error, otherwise number of elements encountered + runVal = run_diagonal(image, mask, size, strides, angles, Na, glrlm, glrlm_idx_max, Nr, jd, a); + if (runVal < 0) return runVal; // index out of range, 'raise' error. + if (!multiElement && runVal > 1) multiElement = 1; // multiple elements found + } + } + for (d1 = 1; d1 < size[mDims[2]]; d1++) + { + for (d2 = 1; d2 < size[mDims[0]]; d2++) + { + jd[mDims[2]] = d1; + jd[mDims[0]] = d2; + // Prevent calculating the true diagonal twice, which is either at index 0, or max, depending on the + // angle of the the second moving dimension: iterate 1 to max if angle is positive, 0 to max -1 if + // negative. The ignored index is handled in the previous loops. This is done for both 'statics'. + if (angles[a * 3 + mDims[2]] < 0) jd[mDims[2]] --; + if (angles[a * 3 + mDims[0]] < 0) jd[mDims[0]] --; + if (angles[a * 3 + mDims[1]] < 0) + { + jd[mDims[1]] = size[mDims[1]] - 1; // Moving dim is negative, set start at maximum index + } + else + { + jd[mDims[1]] = 0; // Moving dim is positive, set start at 0 + } + // run_diagonal returns -1 if error, otherwise number of elements encountered + runVal = run_diagonal(image, mask, size, strides, angles, Na, glrlm, glrlm_idx_max, Nr, jd, a); + if (runVal < 0) return runVal; // index out of range, 'raise' error. + if (!multiElement && runVal > 1) multiElement = 1; // multiple elements found + } + } + } // switch No of dimensions + if (!multiElement) // Segmentation is 2D for this angle, remove it. + { + // set elements at runlength index 0 to 0 for all gray levels in this angle, other runlength indices are + // already 0. + for (d1 = 0; d1 < Ng; d1++) + { + glrlm[a + d1 * Na * Nr] = 0; + } + } + } // next a + return 1; +} + +int run_diagonal(int *image, char *mask, int *size, int *strides, int *angles, int Na, double *glrlm, int glrlm_idx_max, int Nr, int *jd, int a) +{ + int j, gl, rl; + int glrlm_idx; + int elements = 0; + + gl = -1; + rl = 0; + + j = jd[0] * strides[0] + jd[1] * strides[1] + jd[2] * strides[2]; + while (jd[0] >= 0 && jd[0] < size[0] && + jd[1] >= 0 && jd[1] < size[1] && + jd[2] >= 0 && jd[2] < size[2]) + { + j = jd[0] * strides[0] + jd[1] * strides[1] + jd[2] * strides[2]; + if (mask[j]) + { + elements++; // Count the number of segmented voxels in this run + if (gl == image[j]) + { + rl++; + } + else if (gl == -1) + { + gl = image[j]; + rl = 0; + } + else + { + glrlm_idx = a + rl * Na + (gl - 1) * Na * Nr; + if (glrlm_idx >= glrlm_idx_max) return -1; + glrlm[glrlm_idx]++; + + gl = image[j]; + rl = 0; + } + } + else if (gl > -1) + { + glrlm_idx = a + rl * Na + (gl - 1) * Na * Nr; + if (glrlm_idx >= glrlm_idx_max) return -1; + glrlm[glrlm_idx]++; + gl = -1; + } + jd[0] += angles[a * 3]; + jd[1] += angles[a * 3 + 1]; + jd[2] += angles[a * 3 + 2]; + } + if (gl > -1) + { + glrlm_idx = a + rl * Na + (gl - 1) * Na * Nr; + if (glrlm_idx >= glrlm_idx_max) return -1; + glrlm[glrlm_idx]++; + } + return elements; +} diff --git a/radiomics/src/cmatrices.h b/radiomics/src/cmatrices.h new file mode 100644 index 00000000..2e9c3b04 --- /dev/null +++ b/radiomics/src/cmatrices.h @@ -0,0 +1,5 @@ +int calculate_glcm(int *image, char *mask, int Sx, int Sy, int Sz, int *angles, int Na, double *glcm, int Ng); +int calculate_glszm(int *image, signed char *mask, int Sx, int Sy, int Sz, int *angles, int Na, int *tempData, int Ng, int Ns); +int fill_glszm(int *tempData, double *glszm, int Ng, int maxRegion); +int calculate_glrlm(int *image, char *mask, int *size, int *strides, int *angles, int Na, double *glrlm, int Ng, int Nr); +int run_diagonal(int *image, char *mask, int *size, int *strides, int *angles, int Na, double *glrlm, int glrlm_idx_max, int Nr, int *jd, int a); diff --git a/radiomics/src/cshape.c b/radiomics/src/cshape.c new file mode 100644 index 00000000..553db96e --- /dev/null +++ b/radiomics/src/cshape.c @@ -0,0 +1,239 @@ +#include "cshape.h" +#include + +// Declare the look-up tables, these are filled at the bottom of this code file. +static const int gridAngles[8][3]; +static const int edgeTable[128]; +static const int triTable[128][16]; + +double calculate_surfacearea(char *mask, int *size, double *spacing) +{ + int iz, iy, ix, i, j, t, d; // iterator indices + unsigned char cube_idx; // cube identifier, 8 bits signifying which corners of the cube belong to the segmentation + int ang; // Angle index (8 'angles', one pointing to each corner of the marching cube + double vertList[12][3]; // array to hold location of isosurface vectors (points) on the marching cube edges. + double sum; + double surfaceArea = 0; // Total surface area + double a[3], b[3], c[3]; // 2 points of the triangle, relative to the third, and the cross product vector + + // Iterate over all voxels, do not include last voxels in the three dimensions, as the cube includes voxels at pos +1 + for (iz = 0; iz < (size[0] - 1); iz++) + { + for (iy = 0; iy < (size[1] - 1); iy++) + { + for (ix = 0; ix < (size[2] - 1); ix++) + { + // Get current cube_idx by analyzing each point of the current cube + cube_idx = 0; + i = ix + iy * size[2] + iz * size[2] * size[1]; // voxel at cube corner (0, 0, 0) + for (ang = 0; ang < 8; ang++) + { + j = i + gridAngles[ang][0] * size[2] * size[1] + gridAngles[ang][1] * size[2] + gridAngles[ang][2]; + if (mask[j]) cube_idx |= (1 << ang); + } + // Isosurface is symmetrical around the midpoint, flip the number if > 128 + // This enables look-up tables to be 1/2 the size. + if (cube_idx & 0x80) cube_idx ^= 0xff; + + // Exlcude cubes entirely outside or inside the segmentation. Also exclude potential invalid values (> 128). + if (cube_idx > 0 && cube_idx < 128) + { + if (edgeTable[cube_idx] & 1) interpolate(vertList[0], 1, 0, spacing); + if (edgeTable[cube_idx] & 2) interpolate(vertList[1], 2, 1, spacing); + if (edgeTable[cube_idx] & 4) interpolate(vertList[2], 2, 3, spacing); + if (edgeTable[cube_idx] & 8) interpolate(vertList[3], 3, 0, spacing); + if (edgeTable[cube_idx] & 16) interpolate(vertList[4], 5, 4, spacing); + if (edgeTable[cube_idx] & 32) interpolate(vertList[5], 6, 5, spacing); + if (edgeTable[cube_idx] & 64) interpolate(vertList[6], 6, 7, spacing); + if (edgeTable[cube_idx] & 128) interpolate(vertList[7], 7, 4, spacing); + if (edgeTable[cube_idx] & 256) interpolate(vertList[8], 4, 0, spacing); + if (edgeTable[cube_idx] & 512) interpolate(vertList[9], 5, 1, spacing); + if (edgeTable[cube_idx] & 1024) interpolate(vertList[10], 6, 2, spacing); + if (edgeTable[cube_idx] & 2048) interpolate(vertList[11], 7, 3, spacing); + + t = 0; + while (triTable[cube_idx][t*3] >= 0) // Exit loop when no more triangles are present (element at index = -1) + { + for (d = 0; d < 3; d++) + { + a[d] = vertList[triTable[cube_idx][t*3 + 1]][d] - vertList[triTable[cube_idx][t*3]][d]; + b[d] = vertList[triTable[cube_idx][t*3 + 2]][d] - vertList[triTable[cube_idx][t*3]][d]; + } + c[0] = (a[1] * b[2]) - (b[1] * a[2]); + c[1] = (a[2] * b[0]) - (b[2] * a[0]); + c[2] = (a[0] * b[1]) - (b[0] * a[1]); + + // Get the square + c[0] = c[0] * c[0]; + c[1] = c[1] * c[1]; + c[2] = c[2] * c[2]; + + // Compute the surface, which is equal to 1/2 magnitude of the cross product, where + // The magnitude is obtained by calculating the euclidean distance between (0, 0, 0) + // and the location of c + sum = c[0] + c[1] + c[2]; + sum = sqrt(sum); + sum = 0.5 * sum; + surfaceArea += sum; + t++; + } + } + } + } + } + return surfaceArea; +} + +void interpolate(double *vertEntry, int a1, int a2, double *spacing) +{ + int d; + + for (d = 0; d < 3; d++) + { + if (gridAngles[a1][d] == 1 && gridAngles[a2][d] == 1) vertEntry[d] = spacing[d]; + else if (gridAngles[a1][d] != gridAngles[a2][d] ) vertEntry[d] = 0.5 * spacing[d]; + else vertEntry[d] = 0; + } +} + +static const int gridAngles[8][3] = { { 0, 0, 0 }, { 0, 0, 1 }, { 0, 1, 1 }, {0, 1, 0}, { 1, 0, 0 }, {1, 0, 1 }, { 1, 1, 1 }, { 1, 1, 0 }}; +static const int edgeTable[128] = { + 0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, + 0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, + 0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, + 0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, + 0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, + 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, + 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, + 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0 +}; +static const int triTable[128][16] = { + { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1 }, + { 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1 }, + { 3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1 }, + { 3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1 }, + { 2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1 }, + { 8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1 }, + { 4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1 }, + { 3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1 }, + { 4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1 }, + { 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1 }, + { 5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1 }, + { 2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1 }, + { 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1 }, + { 2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1 }, + { 10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1 }, + { 4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1 }, + { 5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1 }, + { 5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1 }, + { 10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1 }, + { 8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1 }, + { 2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1 }, + { 7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1 }, + { 2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1 }, + { 11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1 }, + { 5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1 }, + { 11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1 }, + { 11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1 }, + { 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1 }, + { 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1 }, + { 2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1 }, + { 5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1 }, + { 6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1 }, + { 3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1 }, + { 6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1 }, + { 5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1 }, + { 10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1 }, + { 6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1 }, + { 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1 }, + { 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1 }, + { 3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1 }, + { 5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1 }, + { 0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1 }, + { 9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1 }, + { 8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1 }, + { 5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1 }, + { 0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1 }, + { 6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1 }, + { 10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1 }, + { 10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1 }, + { 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1 }, + { 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1 }, + { 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1 }, + { 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1 }, + { 10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1 }, + { 3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1 }, + { 6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1 }, + { 9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1 }, + { 8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1 }, + { 3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1 }, + { 6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1 }, + { 0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1 }, + { 10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1 }, + { 10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1 }, + { 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1 }, + { 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1 }, + { 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1 }, + { 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1 }, + { 2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1 }, + { 1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1 }, + { 11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1 }, + { 8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1 }, + { 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, + { 7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1 }, + { 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } +}; diff --git a/radiomics/src/cshape.h b/radiomics/src/cshape.h new file mode 100644 index 00000000..81bb5206 --- /dev/null +++ b/radiomics/src/cshape.h @@ -0,0 +1,2 @@ +double calculate_surfacearea(char *mask, int *size, double *spacing); +void interpolate(double *vertEntry, int a1, int a2, double *spacing); diff --git a/setup.py b/setup.py index d296a62f..d0014cd7 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,10 @@ #!/usr/bin/env python -from setuptools import setup +from distutils import sysconfig + +import numpy + +from setuptools import Extension, setup from setuptools.command.test import test as TestCommand import versioneer @@ -11,6 +15,8 @@ with open('requirements-dev.txt', 'r') as fp: dev_requirements = list(filter(bool, (line.strip() for line in fp))) +incDirs = [sysconfig.get_python_inc(), numpy.get_include()] + class NoseTestCommand(TestCommand): """Command to run unit tests using nose driver after in-place build""" @@ -58,6 +64,12 @@ def run_tests(self): 'pyradiomicsbatch=radiomics.scripts.commandlinebatch:main' ]}, + ext_modules=[Extension("_cmatrices", + ["radiomics/src/_cmatrices.c", "radiomics/src/cmatrices.c"], + include_dirs=incDirs, extra_compile_args=['-std=c99']), + Extension("_cshape", ["radiomics/src/_cshape.c", "radiomics/src/cshape.c"], + include_dirs=incDirs)], + description='Radiomics features library for python', license='Slicer', @@ -77,6 +89,7 @@ def run_tests(self): keywords='radiomics cancerimaging medicalresearch computationalimaging', install_requires=requirements, + setup_requires=['numpy>=1.9.2'], test_suite='nose.collector', tests_require=dev_requirements ) diff --git a/tests/test_cmatrices.py b/tests/test_cmatrices.py new file mode 100644 index 00000000..8eb29f9e --- /dev/null +++ b/tests/test_cmatrices.py @@ -0,0 +1,77 @@ +# to run this test, from directory above: +# setenv PYTHONPATH /path/to/pyradiomics/radiomics +# nosetests --nocapture -v tests/test_features.py + +import logging + +from nose_parameterized import parameterized +import numpy + +from radiomics import glcm, glrlm, glszm +from testUtils import custom_name_func, RadiomicsTestUtils + + +testUtils = RadiomicsTestUtils() +defaultTestCases = testUtils.getTestCases() +defaultFeatures = ["glcm", "glrlm", "glszm", "shape"] + +testCases = defaultTestCases +features = ["glcm", "glrlm", "glszm"] # defaultFeatures + + +class TestFeatures: + + def generate_scenarios(): + global testCases + global features + + for testCase in testCases: + for feature in features: + assert(feature is not None) + logging.debug('generate_scenarios: featureClass = %s', feature) + yield testCase, feature + + global testUtils + + @parameterized.expand(generate_scenarios(), testcase_func_name=custom_name_func) + def test_scenario(self, testCase, featureClassName): + print("") + global testUtils + + logging.debug('test_scenario: testCase = %s, featureClassName = %s', testCase, featureClassName) + + testUtils.setFeatureClassAndTestCase(featureClassName, testCase) + + testImage = testUtils.getImage() + testMask = testUtils.getMask() + + pyMat = None + cMat = None + featureClass = None + + if featureClassName == 'firstorder': + logging.debug('No C implementation in firstorder, skipping test') + elif featureClassName == 'glcm': + logging.debug('Init GLCM') + featureClass = glcm.RadiomicsGLCM(testImage, testMask, **testUtils.getKwargs()) + cMat = featureClass.P_glcm + pyMat = featureClass._calculateGLCM() + elif featureClassName == 'glrlm': + logging.debug('Init GLRLM') + featureClass = glrlm.RadiomicsGLRLM(testImage, testMask, **testUtils.getKwargs()) + cMat = featureClass.P_glrlm + pyMat = featureClass._calculateGLRLM() + elif featureClassName == 'shape': + logging.debug('No C implementation in glrlm, skipping test') + # No C implementation yet, will follow + # logging.debug('Init Shape') + # featureClass = shape.RadiomicsShape(testImage, testMask, **testUtils.getKwargs()) + elif featureClassName == 'glszm': + logging.debug('Init GLSZM') + featureClass = glszm.RadiomicsGLSZM(testImage, testMask, **testUtils.getKwargs()) + cMat = featureClass.P_glszm + pyMat = featureClass._calculateGLSZM() + + assert (featureClass is not None) + # Check if the calculated arrays match + assert numpy.max(numpy.abs(pyMat - cMat)) < 1e-3