ENH: Improve search strategy in C

Improve the search that finds a new index to process in the region growing in C calculation of glszm.

circle.yml

@@ -8,5 +8,5 @@ dependencies:
     - pip install PyWavelets==0.4.0
 test:
   override:
-    - nosetests --with-xunit --logging-level=DEBUG --verbosity=3 tests/test_features.py
+    - nosetests --with-xunit --logging-level=DEBUG --verbosity=3 tests/test_features.py tests/test_cmatrices.py
     - cp nosetests.xml $CIRCLE_TEST_REPORTS

radiomics/glcm.py

@@ -107,16 +107,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.targetVoxelArray, self.matrix, self.matrixCoordinates)
     self.coefficients['Ng'] = self.histogram[1].shape[0] - 1
 
     if radiomics.debugging:
-      self._calculateGLCM()
+      self.P_glcm = self._calculateGLCM()
     else:
-      self._calculateCGLCM()
+      self.P_glcm = self._calculateCGLCM()
     self._calculateCoefficients()
 
   def _calculateGLCM(self):
@@ -134,7 +134,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')
 
@@ -159,12 +159,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:
@@ -181,28 +181,28 @@ 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(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
 
   def _calculateCGLCM(self):
     size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
     angles = imageoperations.generateAngles(size)
     Ng = self.coefficients['Ng']
 
-    self.P_glcm = _cmatrices.calculate_glcm(self.matrix, self.maskArray, angles, Ng)
+    P_glcm = _cmatrices.calculate_glcm(self.matrix, self.maskArray, angles, Ng)
 
     # 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:
@@ -219,17 +219,17 @@ def _calculateCGLCM(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(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):

radiomics/gldm.py

@@ -2,6 +2,8 @@
 import radiomics
 from radiomics import base, imageoperations
 import _cmatrices
+if radiomics.debugging:  # Only import tqdm if using python matrix calculation
+  from tqdm import trange
 
 class RadiomicsGLDM(base.RadiomicsFeaturesBase):
   r"""
@@ -57,16 +59,16 @@ def __init__(self, inputImage, inputMask, **kwargs):
     self.gldm_a = kwargs.get('gldm_a', 0)
 
     self.coefficients = {}
-    self.P_gldm = {}
+    self.P_gldm = None
 
     # binning
     self.matrix, self.histogram = imageoperations.binImage(self.binWidth, self.targetVoxelArray, self.matrix, self.matrixCoordinates)
     self.coefficients['Ng'] = self.histogram[1].shape[0] - 1
 
     if radiomics.debugging:
-      self._calculateGLDM()
+      self.P_gldm = self._calculateGLDM()
     else:
-      self._calculateCGLDM()
+      self.P_gldm = self._calculateCGLDM()
     self._calculateCoefficients()
 
   def _calculateGLDM(self):
@@ -81,15 +83,15 @@ def _calculateGLDM(self):
     size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
     angles = imageoperations.generateAngles(size)
     angles = numpy.concatenate((angles, angles * -1))
-    Nd = len(angles) + 1
 
-    diffMat = numpy.empty((self.matrix.shape + (Nd,)))
+    depMat = numpy.zeros(self.matrix.shape, dtype='int')
 
+    if self.verbose: bar = trange(len(angles), desc='Calculate shifted matrices (GLDM)')
     for a_idx, a in enumerate(angles):
+      if self.verbose: bar.update()
       # create shifted array (by angle), so that for an index idx, angMat[idx] is the neigbour of self.matrix[idx]
       # for the current angle.
-      angMat = diffMat[:, :, :, a_idx]
-      angMat[:, :, :] = numpy.roll(numpy.roll(numpy.roll(self.matrix, -a[0], 0), -a[1], 1), -a[2], 2)
+      angMat = numpy.abs(numpy.roll(numpy.roll(numpy.roll(self.matrix, -a[0], 0), -a[1], 1), -a[2], 2) - self.matrix)
       if a[0] > 0:
         angMat[-a[0]:, :, :] = padVal
       elif a[0] < 0:
@@ -105,35 +107,31 @@ def _calculateGLDM(self):
       elif a[2] < 0:
         angMat[:, :, :-a[2]] = padVal
 
-    # Create boolean array with shape equal to self.matrix, where elements are true if |depMat-self.matrix| <=
-    # numpy.where needed to prevent warnings caused by comparison with NaN
-    # sum over angles axis to get total number of dependent voxels for all matrix voxels
-    diffMat = numpy.abs(diffMat - self.matrix[:, :, :, None])
-    nanMask = numpy.isnan(diffMat)
-
-    depMat = numpy.zeros((self.matrix.shape + (Nd,)), dtype='int')
-    depMat[~nanMask] = (diffMat[~nanMask] <= self.gldm_a)
-    depMat = numpy.sum(depMat, 3)
-
-    P_gldm = numpy.zeros((Ng, Nd))
-    for i in xrange(1, Ng + 1):
+      nanMask = numpy.isnan(angMat)
+      depMat[~nanMask] += (angMat[~nanMask] <= self.gldm_a)
+
+    if self.verbose: bar.close()
+
+    Nd = numpy.max(depMat)
+    P_gldm = numpy.zeros((Ng, Nd + 1))
+    grayLevels = numpy.unique(self.matrix[self.matrixCoordinates])
+
+    if self.verbose: bar = trange(len(grayLevels), desc= 'calculate GLDM')
+    for i in grayLevels:
+        if self.verbose: bar.update()
         i_mat = (self.matrix == i)
-        for d in numpy.unique(depMat):
+        for d in numpy.unique(depMat[i_mat]):
             # By multiplying i_mat and depMat == d, a boolean area is obtained,
             # where the number of elements that are true (1) is equal to the number of voxels
             # with gray level i and dependence d.
             P_gldm[i-1, d] = numpy.sum(i_mat * (depMat == d))
+    if self.verbose: bar.close()
 
-    # Crop gray-level axis of GLDM matrix to between minimum and maximum observed gray-levels
-    # Crop dependence axis of GLDM matrix up to maximum observed dependence
-    P_gldm_bounds = numpy.argwhere(P_gldm)
-    (xstart, ystart), (xstop, ystop) = P_gldm_bounds.min(0), P_gldm_bounds.max(0) + 1
-    self.P_gldm = P_gldm[xstart:xstop,:ystop]
+    return P_gldm
 
   def _calculateCGLDM(self):
     size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
     angles = imageoperations.generateAngles(size)
-    Nd = len(angles)
     Ng = self.coefficients['Ng']
 
     P_gldm = _cmatrices.calculate_gldm(self.matrix, self.maskArray, angles, Ng, self.gldm_a)
@@ -142,7 +140,7 @@ def _calculateCGLDM(self):
     # Crop dependence axis of GLDM matrix up to maximum observed dependence
     P_gldm_bounds = numpy.argwhere(P_gldm)
     (xstart, ystart), (xstop, ystop) = P_gldm_bounds.min(0), P_gldm_bounds.max(0) + 1
-    self.P_gldm = P_gldm[xstart:xstop, :ystop]
+    return P_gldm[xstart:xstop, :ystop]
 
   def _calculateCoefficients(self):
     sumP_gldm = numpy.sum(self.P_gldm, (0, 1))

radiomics/glszm.py

@@ -51,17 +51,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.targetVoxelArray, self.matrix, self.matrixCoordinates)
     self.coefficients['Ng'] = self.histogram[1].shape[0] - 1
     self.coefficients['Np'] = self.targetVoxelArray.size
 
     if radiomics.debugging:
-      self._calculateGLSZM()
+      self.P_glszm = self._calculateGLSZM()
     else:
-      self._calculateCGLSZM()
+      self.P_glszm = self._calculateCGLSZM()
     self._calculateCoefficients()
 
   def _calculateGLSZM(self):
@@ -127,21 +127,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
-    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']
 
-    P_glszm = _cmatrices.calculate_glszm(self.matrix, self.maskArray, angles, Ng, Ns)
-
-    # Crop gray-level axis of GLSZM matrix to between minimum and maximum observed gray-levels
-    # 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
-    self.P_glszm = P_glszm[xstart:xstop, :ystop]
+    return _cmatrices.calculate_glszm(self.matrix, self.maskArray, angles, Ng, Ns)
 
   def _calculateCoefficients(self):
     sumP_glszm = numpy.sum(self.P_glszm, (0, 1) )

radiomics/ngtdm.py

@@ -2,6 +2,8 @@
 import radiomics
 from radiomics import base, imageoperations
 import _cmatrices
+if radiomics.debugging:  # Only import tqdm if using python matrix calculation
+  from tqdm import trange
 
 
 class RadiomicsNGTDM(base.RadiomicsFeaturesBase):
@@ -82,6 +84,7 @@ def __init__(self, inputImage, inputMask, **kwargs):
     super(RadiomicsNGTDM, self).__init__(inputImage, inputMask, **kwargs)
 
     self.coefficients = {}
+    self.P_ngtdm = None
 
     # binning
     self.matrix, self.histogram = imageoperations.binImage(self.binWidth, self.targetVoxelArray, self.matrix,
@@ -90,9 +93,11 @@ def __init__(self, inputImage, inputMask, **kwargs):
     self.coefficients['Np'] = self.targetVoxelArray.size
 
     if radiomics.debugging:
-      self._calculateNGTDM()
+      self.P_ngtdm = self._calculateNGTDM()
     else:
-      self._calculateCNGTDM()
+      self.P_ngtdm = self._calculateCNGTDM()
+
+    self._calculateCoefficients()
 
   def _calculateNGTDM(self):
     Ng = self.coefficients['Ng']
@@ -108,13 +113,15 @@ def _calculateNGTDM(self):
     angles = numpy.concatenate((angles, angles * -1))
     Nd = len(angles)
 
-    angMats = numpy.empty((self.matrix.shape + (Nd,)), dtype='float')
+    dataTemp = numpy.zeros(self.matrix.shape, dtype='float')
+    countMat = numpy.zeros(self.matrix.shape, dtype='int')
 
+    if self.verbose: bar = trange(Nd-1, desc='Calculate shifted matrices (NGTDM)')
     for a_idx, a in enumerate(angles):
+      if self.verbose: bar.update()
       # create shifted array (by angle), so that for an index idx, angMat[idx] is the neigbour of self.matrix[idx]
       # for the current angle.
-      angMat = angMats[:, :, :, a_idx]
-      angMat[:, :, :] = numpy.roll(numpy.roll(numpy.roll(self.matrix, -a[0], 0), -a[1], 1), -a[2], 2)
+      angMat = numpy.roll(numpy.roll(numpy.roll(self.matrix, -a[0], 0), -a[1], 1), -a[2], 2)
       if a[0] > 0:
           angMat[-a[0]:, :, :] = padVal
       elif a[0] < 0:
@@ -129,21 +136,34 @@ def _calculateNGTDM(self):
           angMat[:, :, -a[2]:] = padVal
       elif a[2] < 0:
           angMat[:, :, :-a[2]] = padVal
+      nanmask = numpy.isnan(angMat)
+      dataTemp[~nanmask] += angMat[~nanmask]
+      countMat[~nanmask] += 1
+
+    if self.verbose: bar.close()
+
+    # Average neighbourhood is the dataTemp (which is the sum of gray levels of neighbours that are non-NaN) divided by
+    # the countMat (which is the number of neighbours that are non-NaN)
+    nZeroMask = countMat > 0  # Prevent division by 0
+    dataTemp[nZeroMask] = dataTemp[nZeroMask] / countMat[nZeroMask]
 
-    # Create a difference matrix by taking the absolut of self.matrix - the neighbourhood mean, excluding NaNs
-    # in the calculation of the mean
-    diffMat = numpy.abs(self.matrix - numpy.nanmean(angMats, axis=3))
+    dataTemp = numpy.abs(dataTemp - self.matrix)  # Calculate the absolute difference
+
+    P_ngtdm = numpy.zeros((Ng, 3), dtype='float')
 
     # For each gray level present in self.matrix:
     # element 0 = probability of gray level (p_i),
     # element 1 = sum of the absolute differences (s_i),
     # element 2 = gray level (i)
-    P_ngtdm = numpy.zeros((Ng, 3), dtype='float')
-    for i in numpy.unique(self.matrix):
+    grayLevels = numpy.unique(self.matrix[self.matrixCoordinates])
+    if self.verbose: bar = trange(len(grayLevels), desc='Calculate NGTDM')
+    for i in grayLevels:
+        if self.verbose: bar.update()
         if not numpy.isnan(i):
           i_ind = numpy.where(self.matrix == i)
           P_ngtdm[int(i-1), 0] = len(i_ind[0])
-          P_ngtdm[int(i-1), 1] = numpy.sum(diffMat[i_ind])
+          P_ngtdm[int(i-1), 1] = numpy.sum(dataTemp[i_ind])
+    if self.verbose: bar.close()
 
     # Fill in gray levels (needed as empty gray level slices will be deleted)
     P_ngtdm[:, 2] = numpy.arange(1, Ng+1)
@@ -152,13 +172,9 @@ def _calculateNGTDM(self):
     P_ngtdm = numpy.delete(P_ngtdm, numpy.where(P_ngtdm[:, 0] == 0), 0)
 
     # Normalize P_ngtdm[:, 0] (= p_i)
-    self.coefficients['p_i'] = P_ngtdm[:, 0] / self.coefficients['Np']
-
-    self.coefficients['s_i'] = P_ngtdm[:, 1]
-    self.coefficients['ivector'] = P_ngtdm[:, 2]
+    P_ngtdm[:, 0] = P_ngtdm[:, 0] / self.coefficients['Np']
 
-    # Ngp = number of graylevels, for which p_i > 0
-    self.coefficients['Ngp'] = P_ngtdm.shape[0]
+    return P_ngtdm
 
   def _calculateCNGTDM(self):
     size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1
@@ -171,13 +187,18 @@ def _calculateCNGTDM(self):
     P_ngtdm = numpy.delete(P_ngtdm, numpy.where(P_ngtdm[:, 0] == 0), 0)
 
     # Normalize P_ngtdm[:, 0] (= p_i)
-    self.coefficients['p_i'] = P_ngtdm[:, 0] / self.coefficients['Np']
+    P_ngtdm[:, 0] = P_ngtdm[:, 0] / self.coefficients['Np']
+
+    return P_ngtdm
+
+  def _calculateCoefficients(self):
+    self.coefficients['p_i'] = self.P_ngtdm[:, 0]
 
-    self.coefficients['s_i'] = P_ngtdm[:, 1]
-    self.coefficients['ivector'] = P_ngtdm[:, 2]
+    self.coefficients['s_i'] = self.P_ngtdm[:, 1]
+    self.coefficients['ivector'] = self.P_ngtdm[:, 2]
 
     # Ngp = number of graylevels, for which p_i > 0
-    self.coefficients['Ngp'] = P_ngtdm.shape[0]
+    self.coefficients['Ngp'] = self.P_ngtdm.shape[0]
 
   def getCoarsenessFeatureValue(self):
     r"""