Enhance C GLSZM

CHANGES.rst

@@ -16,6 +16,9 @@ Feature Calculation Changes
   (`#233 <https://github.com/Radiomics/pyradiomics/pull/233>`_)
 - Redefine features *Elongation* and *Flatness* as the inverse of the original definition. This prevents a returned
   value of NaN when the shape is completely flat. (`#234 <https://github.com/Radiomics/pyradiomics/pull/234>`_)
+- In certain edge cases, the calculated maximum diameters may be too small when calculating using the python
+  implementation. This is corrected by the C extension and a warning is now logged when calculating these features in
+  python. (`#255 <https://github.com/Radiomics/pyradiomics/pull/255>`_)
 
 New Features
 ############
@@ -64,6 +67,12 @@ Internal API
 ############
 
 - Add function to get or download test case (`#235 <https://github.com/Radiomics/pyradiomics/pull/235>`_)
+- Rewrite C Extension algorithm for GSLZM. Instead of searching over the image for the next voxel when
+  growing a region, store all unprocessed voxels in a stack. This yields a significant increase in performance,
+  especially in large ROIs. Requires slightly more memory (1 array, type integer, size equal to number of voxels in
+  the ROI) (`#255 <https://github.com/Radiomics/pyradiomics/pull/255>`_)
+- Implement C extension for calculation of maximum diameters.
+  (`#255 <https://github.com/Radiomics/pyradiomics/pull/255>`_)
 
 Cleanups
 ########

radiomics/shape.py

@@ -2,7 +2,7 @@
 import SimpleITK as sitk
 from six.moves import range
 
-from radiomics import base, cMatsEnabled, cShape
+from radiomics import base, cMatsEnabled, cShape, imageoperations
 
 
 class RadiomicsShape(base.RadiomicsFeaturesBase):
@@ -26,7 +26,6 @@ def __init__(self, inputImage, inputMask, **kwargs):
     self.logger.debug('Extracting simple ITK shape features')
 
     self.lssif = sitk.LabelShapeStatisticsImageFilter()
-    self.lssif.SetComputeFeretDiameter(True)
     self.lssif.Execute(inputMask)
 
     # Pad inputMask to prevent index-out-of-range errors
@@ -58,6 +57,8 @@ def __init__(self, inputImage, inputMask, **kwargs):
     else:
       self.SurfaceArea = self._calculateSurfaceArea()
 
+    self.diameters = None  # Do not precompute diameters
+
     self.logger.debug('Feature class initialized')
 
   def _calculateSurfaceArea(self):
@@ -138,38 +139,20 @@ def _calculateCSurfaceArea(self):
 
     return cShape.calculate_surfacearea(self.maskArray, self.pixelSpacing)
 
-  def _getMaximum2Ddiameter(self, dim):
-    otherDims = tuple(set([0, 1, 2]) - set([dim]))
-
-    a = numpy.array(list(zip(*self.matrixCoordinates)))
-
-    maxDiameter = 0
-    # Check maximum diameter in every slice, retain the overall maximum
-    for i in numpy.unique(a[:, dim]):
-      # Retrieve all indices of mask in current slice
-      plane = a[numpy.where(a[:, dim] == i)]
-
-      minBounds = numpy.min(plane, 0)
-      maxBounds = numpy.max(plane, 0)
-
-      # Generate 2 sets of indices: one set of indices in zSlice where at least the x or y component of the index is equal to the
-      # minimum indices in the current slice, and one set of indices where at least one element it is equal to the maximum
-      edgeVoxelsMinCoords = numpy.vstack(
-        [plane[plane[:, otherDims[0]] == minBounds[otherDims[0]]],
-         plane[plane[:, otherDims[1]] == minBounds[otherDims[1]]]]) * self.pixelSpacing
-      edgeVoxelsMaxCoords = numpy.vstack(
-        [plane[plane[:, otherDims[0]] == maxBounds[otherDims[0]]],
-         plane[plane[:, otherDims[1]] == maxBounds[otherDims[1]]]]) * self.pixelSpacing
-
-      # generate a matrix of distances for every combination of an index in edgeVoxelsMinCoords and edgeVoxelsMaxCoords
-      # By subtraction the distance between the x, y and z components are obtained. The euclidean distance is then calculated:
-      # Sum the squares of dx, dy and dz components and then take the square root; Sqrt( Sum( dx^2 + dy^2 + dz^2 ) )
-      distances = numpy.sqrt(numpy.sum((edgeVoxelsMaxCoords[:, None] - edgeVoxelsMinCoords[None, :]) ** 2, 2))
-      tempMax = numpy.max(distances)
-      if tempMax > maxDiameter:
-        maxDiameter = tempMax
+  def _calculateCDiameters(self):
+    """
+    Calculate maximum diameters in 2D and 3D using C extension. Function returns a tuple with 4 elements:
 
-    return maxDiameter
+    0. Maximum 2D diameter Slice (XY Plane, Axial)
+    1. Maximum 2D diameter Column (ZX Plane, Coronal)
+    2. Maximum 2D diameter Row (ZY Plane, Sagittal)
+    3. Maximum 3D diameter
+    """
+    self.logger.debug('Calculating Maximum 3D diameter in C')
+    Ns = self.targetVoxelArray.size
+    size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
+    angles = imageoperations.generateAngles(size)
+    return cShape.calculate_diameter(self.maskArray, self.pixelSpacing, angles, Ns)
 
   def getVolumeFeatureValue(self):
     r"""
@@ -314,38 +297,80 @@ def getMaximum3DDiameterFeatureValue(self):
     Maximum 3D diameter is defined as the largest pairwise Euclidean distance between surface voxels in the ROI.
 
     Also known as Feret Diameter.
+
+    .. warning::
+
+      This feature is only available when C Extensions are enabled
     """
-    return self.lssif.GetFeretDiameter(self.label)
+
+    if cMatsEnabled():
+      if self.diameters is None:
+        self.diameters = self._calculateCDiameters()
+      return self.diameters[3]
+    else:
+      self.logger.warning('For computational reasons, this feature is only implemented in C. Enable C extensions to '
+                          'calculate this feature.')
+      return numpy.nan
 
   def getMaximum2DDiameterSliceFeatureValue(self):
     r"""
     **9. Maximum 2D diameter (Slice)**
 
     Maximum 2D diameter (Slice) is defined as the largest pairwise Euclidean distance between tumor surface voxels in
     the row-column (generally the axial) plane.
-    """
 
-    return self._getMaximum2Ddiameter(0)
+    .. warning::
+
+      This feature is only available when C Extensions are enabled
+    """
+    if cMatsEnabled():
+      if self.diameters is None:
+        self.diameters = self._calculateCDiameters()
+      return self.diameters[0]
+    else:
+      self.logger.warning('For computational reasons, this feature is only implemented in C. Enable C extensions to '
+                          'calculate this feature.')
+      return numpy.nan
 
   def getMaximum2DDiameterColumnFeatureValue(self):
     r"""
     **10. Maximum 2D diameter (Column)**
 
     Maximum 2D diameter (Column) is defined as the largest pairwise Euclidean distance between tumor surface voxels in
     the row-slice (usually the coronal) plane.
-    """
 
-    return self._getMaximum2Ddiameter(1)
+    .. warning::
+
+      This feature is only available when C Extensions are enabled
+    """
+    if cMatsEnabled():
+      if self.diameters is None:
+        self.diameters = self._calculateCDiameters()
+      return self.diameters[1]
+    else:
+      self.logger.warning('For computational reasons, this feature is only implemented in C. Enable C extensions to '
+                          'calculate this feature.')
+      return numpy.nan
 
   def getMaximum2DDiameterRowFeatureValue(self):
     r"""
     **11. Maximum 2D diameter (Row)**
 
     Maximum 2D diameter (Row) is defined as the largest pairwise Euclidean distance between tumor surface voxels in the
     column-slice (usually the sagittal) plane.
-    """
 
-    return self._getMaximum2Ddiameter(2)
+    .. warning::
+
+      This feature is only available when C Extensions are enabled
+    """
+    if cMatsEnabled():
+      if self.diameters is None:
+        self.diameters = self._calculateCDiameters()
+      return self.diameters[2]
+    else:
+      self.logger.warning('For computational reasons, this feature is only implemented in C. Enable C extensions to '
+                          'calculate this feature.')
+      return numpy.nan
 
   def getMajorAxisFeatureValue(self):
     r"""

radiomics/src/_cmatrices.c

@@ -181,7 +181,7 @@ static PyObject *cmatrices_calculate_glszm(PyObject *self, PyObject *args)
   PyArrayObject *image_arr, *mask_arr, *angles_arr;
   int Sx, Sy, Sz, Na;
   int *image;
-  signed char *mask;
+  char *mask;
   int *angles;
   int *tempData;
   int maxRegion;
@@ -248,7 +248,7 @@ static PyObject *cmatrices_calculate_glszm(PyObject *self, PyObject *args)
 
   // Get arrays in Ctype
   image = (int *)PyArray_DATA(image_arr);
-  mask = (signed char *)PyArray_DATA(mask_arr);
+  mask = (char *)PyArray_DATA(mask_arr);
   angles = (int *)PyArray_DATA(angles_arr);
 
   //Calculate GLSZM

radiomics/src/_cshape.c

@@ -1,3 +1,4 @@
+#include <stdlib.h>
 #include <Python.h>
 #include <numpy/arrayobject.h>
 #include "cshape.h"
@@ -12,12 +13,15 @@ static char module_docstring[] = "This module links to C-compiled code for effic
 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 char diameter_docstring[] = "Arguments: Mask, PixelSpacing, Angles, ROI size.";
 
 static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args);
+static PyObject *cshape_calculate_diameter(PyObject *self, PyObject *args);
 
 static PyMethodDef module_methods[] = {
   //{"calculate_", cmatrices_, METH_VARARGS, _docstring},
   { "calculate_surfacearea", cshape_calculate_surfacearea, METH_VARARGS, surface_docstring },
+  { "calculate_diameter", cshape_calculate_diameter,METH_VARARGS, diameter_docstring},
   { NULL, NULL, 0, NULL }
 };
 
@@ -83,7 +87,6 @@ static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args)
   char *mask;
   double *spacing;
   double SA;
-
   // Parse the input tuple
   if (!PyArg_ParseTuple(args, "OO", &mask_obj, &spacing_obj))
     return NULL;
@@ -118,7 +121,7 @@ static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args)
   mask = (char *)PyArray_DATA(mask_arr);
   spacing = (double *)PyArray_DATA(spacing_arr);
 
-  //Calculate GLCM
+  //Calculate Surface Area
   SA = calculate_surfacearea(mask, size, spacing);
 
   // Clean up
@@ -127,9 +130,84 @@ static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args)
 
   if (SA < 0) // if SA < 0, an error has occurred
   {
-    PyErr_SetString(PyExc_RuntimeError, "Calculation of GLCM Failed.");
+    PyErr_SetString(PyExc_RuntimeError, "Calculation of Surface Area Failed.");
     return NULL;
   }
 
   return Py_BuildValue("f", SA);
 }
+
+static PyObject *cshape_calculate_diameter(PyObject *self, PyObject *args)
+{
+  PyObject *mask_obj, *spacing_obj, *angles_obj;
+  PyArrayObject *mask_arr, *spacing_arr, *angles_arr;
+  int Ns, Na;
+  int size[3];
+  char *mask;
+  double *spacing;
+  int *angles;
+  double *diameters;
+  PyObject *rslt;
+
+  // Parse the input tuple
+  if (!PyArg_ParseTuple(args, "OOOi", &mask_obj, &spacing_obj, &angles_obj, &Ns))
+    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);
+  angles_arr = (PyArrayObject *)PyArray_FROM_OTF(angles_obj, NPY_INT, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY);
+
+  if (mask_arr == NULL || spacing_arr == NULL || angles_arr == NULL)
+  {
+    Py_XDECREF(mask_arr);
+    Py_XDECREF(spacing_arr);
+    Py_XDECREF(angles_arr);
+    return NULL;
+  }
+
+  if (PyArray_NDIM(mask_arr) != 3)
+  {
+    Py_XDECREF(mask_arr);
+    Py_XDECREF(spacing_arr);
+    Py_XDECREF(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(mask_arr, 2);
+  size[1] = (int)PyArray_DIM(mask_arr, 1);
+  size[0] = (int)PyArray_DIM(mask_arr, 0);
+
+  Na = (int)PyArray_DIM(angles_arr, 0);
+
+  // Get arrays in Ctype
+  mask = (char *)PyArray_DATA(mask_arr);
+  spacing = (double *)PyArray_DATA(spacing_arr);
+  angles = (int *)PyArray_DATA(angles_arr);
+
+  // Initialize output array (elements not set)
+  diameters = (double *)calloc(4, sizeof(double));
+
+  // Calculating Max 3D Diameter
+  if (!calculate_diameter(mask, size, spacing, angles, Na, Ns, diameters))
+  {
+    Py_XDECREF(mask_arr);
+    Py_XDECREF(spacing_arr);
+    Py_XDECREF(angles_arr);
+    free(diameters);
+    PyErr_SetString(PyExc_RuntimeError, "Calculation of maximum 3D diameter failed.");
+    return NULL;
+  }
+
+  rslt = Py_BuildValue("ffff", diameters[0], diameters[1], diameters[2], diameters[3]);
+
+  // Clean up
+  Py_XDECREF(mask_arr);
+  Py_XDECREF(spacing_arr);
+  Py_XDECREF(angles_arr);
+  free(diameters);
+
+  return rslt;
+}