ENH: Add Neighbouring Gray Tone Difference Matrix (NGTDM)

data/baseline/baseline_ngtdm.csv

@@ -0,0 +1,6 @@
+Patient ID,general_info_BoundingBox,general_info_GeneralSettings,general_info_ImageHash,general_info_ImageSpacing,general_info_InputImages,general_info_MaskHash,general_info_Version,general_info_VolumeNum,general_info_VoxelNum,Coarseness,Complexity,Strength,Busyness,Contrast
+brain1,"(162, 84, 11, 47, 70, 7)","{'distances': [1], 'voxelArrayShift': 2000, 'additionalInfo': True, 'enableCExtensions': True, 'weightingNorm': None, 'force2D': False, 'interpolator': None, 'resampledPixelSpacing': None, 'gldm_a': 0, 'label': 1, 'normalizeScale': 1, 'normalize': False, 'force2Ddimension': 0, 'removeOutliers': None, 'minimumROISize': None, 'binWidth': 25, 'minimumROIDimensions': 1, 'symmetricalGLCM': False, 'padDistance': 5}",5c9ce3ca174f0f8324aa4d277e0fef82dc5ac566,"(0.7812499999999999, 0.7812499999999999, 6.499999999999998)",{},9dc2c3137b31fd872997d92c9a92d5178126d9d3,1.2.0.post17.dev0+g46c7fb8,2,4137,0.0017120352713566309,0.34796054586526076,0.63073148273486113,1.4467306835644076,5.9914259435373462e-05
+brain2,"(205, 155, 8, 20, 15, 3)","{'distances': [1], 'voxelArrayShift': 2000, 'additionalInfo': True, 'enableCExtensions': True, 'weightingNorm': None, 'force2D': False, 'interpolator': None, 'resampledPixelSpacing': None, 'gldm_a': 0, 'label': 1, 'normalizeScale': 1, 'normalize': False, 'force2Ddimension': 0, 'removeOutliers': None, 'minimumROISize': None, 'binWidth': 25, 'minimumROIDimensions': 1, 'symmetricalGLCM': False, 'padDistance': 5}",f2b8fbc4d5d1da08a1a70e79a301f8a830139438,"(0.7812499999999999, 0.7812499999999999, 6.499999999999998)",{},b41049c71633e194bee4891750392b72eabd8800,1.2.0.post17.dev0+g46c7fb8,1,453,0.01286160865655888,1.6825211314863389,3.9663129914847497,0.15164513751048295,0.0002689553152317465
+breast1,"(21, 64, 8, 9, 12, 3)","{'distances': [1], 'voxelArrayShift': 2000, 'additionalInfo': True, 'enableCExtensions': True, 'weightingNorm': None, 'force2D': False, 'interpolator': None, 'resampledPixelSpacing': None, 'gldm_a': 0, 'label': 1, 'normalizeScale': 1, 'normalize': False, 'force2Ddimension': 0, 'removeOutliers': None, 'minimumROISize': None, 'binWidth': 25, 'minimumROIDimensions': 1, 'symmetricalGLCM': False, 'padDistance': 5}",016951a8f9e8e5de21092d9d62b77262f92e04a5,"(0.664062, 0.664062, 2.1)",{},5aa7d57fd57e83125b605c036c40f4a0d0cfd3e4,1.2.0.post17.dev0+g46c7fb8,1,143,0.050821154965410821,0.0082145909265791753,0.12459116813839596,8.9043951456360233,0.00044669705438426832
+lung1,"(206, 347, 32, 24, 26, 3)","{'distances': [1], 'voxelArrayShift': 2000, 'additionalInfo': True, 'enableCExtensions': True, 'weightingNorm': None, 'force2D': False, 'interpolator': None, 'resampledPixelSpacing': None, 'gldm_a': 0, 'label': 1, 'normalizeScale': 1, 'normalize': False, 'force2Ddimension': 0, 'removeOutliers': None, 'minimumROISize': None, 'binWidth': 25, 'minimumROIDimensions': 1, 'symmetricalGLCM': False, 'padDistance': 5}",34dca4200809a5e76c702d6b9503d958093057a3,"(0.5703125, 0.5703125, 5.0)",{},054d887740012177bd1f9031ddac2b67170af0f3,1.2.0.post17.dev0+g46c7fb8,1,837,0.0089851482891101057,0.73822878390479196,2.7864289040763568,0.19930533973333536,0.00021955754133650324
+lung2,"(318, 333, 15, 87, 66, 11)","{'distances': [1], 'voxelArrayShift': 2000, 'additionalInfo': True, 'enableCExtensions': True, 'weightingNorm': None, 'force2D': False, 'interpolator': None, 'resampledPixelSpacing': None, 'gldm_a': 0, 'label': 1, 'normalizeScale': 1, 'normalize': False, 'force2Ddimension': 0, 'removeOutliers': None, 'minimumROISize': None, 'binWidth': 25, 'minimumROIDimensions': 1, 'symmetricalGLCM': False, 'padDistance': 5}",14f57fd04838eb8c9cca2a0dd871d29971585975,"(0.6269531, 0.6269531, 5.0)",{},e284ff05593bc6cb2747261882e452d4efbccb3a,1.2.0.post17.dev0+g46c7fb8,1,24644,0.00020294841093398597,0.046405126305385076,0.61472046408820225,2.0631991387093871,1.1469959784621812e-06

docs/customization.rst

@@ -238,6 +238,16 @@ Feature Class Level
     For more information on how weighting is applied, see the documentation on :ref:`GLCM <radiomics-glcm-label>` and
     :ref:`GLRLM <radiomics-glszm-label>`.
 
+*Distance to neighbour*
+
+- distances [[1]]: List of integers. This specifies the distances between the center voxel and the neighbor, for which
+  angles should be generated. See also :py:func:`~radiomics.imageoperations.generateAngles`
+
+.. note::
+
+    This only affects the GLCM and NGTDM feature classes. The GLSZM and GLRLM feature classes use a fixed distance of 1
+    (infinity norm) to define neighbours.
+
 Feature Class Specific Settings
 +++++++++++++++++++++++++++++++
 
@@ -250,8 +260,9 @@ Feature Class Specific Settings
 
 *GLCM*
 
-- distances [[1]]: List of integers. This specifies the distances between the center voxel and the neighbor, for which
-  angles should be generated. See also :py:func:`~radiomics.imageoperations.generateAngles`
+- symmetricalGLCM [True]: boolean, indicates whether co-occurrences should be assessed in two directions per angle,
+  which results in a symmetrical matrix, with equal distributions for :math:`i` and :math:`j`. A symmetrical matrix
+  corresponds to the GLCM as defined by Haralick et al.
 
 .. _radiomics-parameter-file-label:
 

docs/features.rst

@@ -12,6 +12,7 @@ subdivided into the following classes:
 * :ref:`radiomics-glcm-label` (26 features)
 * :ref:`radiomics-glszm-label` (16 features)
 * :ref:`radiomics-glrlm-label` (16 features)
+* :ref:`radiomics-ngtdm-label` (5 features)
 
 All feature classes, with the exception of shape can be calculated on either the original image and/or a derived image,
 obtained by applying one of several filters. The shape descriptors are independent of gray value, and are extracted from
@@ -77,5 +78,16 @@ Gray Level Run Length Matrix (GLRLM) Features
     :show-inheritance:
     :member-order: bysource
 
+.. _radiomics-ngtdm-label:
+
+Neigbouring Gray Tone Difference Matrix (NGTDM) Features
+--------------------------------------------------------
+
+.. automodule:: radiomics.ngtdm
+    :members:
+    :undoc-members:
+    :show-inheritance:
+    :member-order: bysource
+
 .. [1] Zwanenburg, A., Leger, S., Vallières, M., and Löck, S. (2016). Image biomarker
     standardisation initiative - feature definitions. In eprint arXiv:1612.07003 [cs.CV]

docs/index.rst

@@ -57,8 +57,9 @@ Currently supports the following feature classes:
 * :ref:`Gray Level Cooccurence Matrix <radiomics-glcm-label>` (GLCM)
 * :ref:`Gray Level Run Length Matrix <radiomics-glrlm-label>` (GLRLM)
 * :ref:`Gray Level Size Zone Matrix <radiomics-glszm-label>` (GLSZM)
+* :ref:`Neigbouring Gray Tone Difference Matrix <radiomics-ngtdm-label>` (NGTDM)
 
-On average, Pyradiomics extracts :math:`\approx 1300` features per image, which consist of the 16 shape descriptors and
+On average, Pyradiomics extracts :math:`\approx 1400` features per image, which consist of the 16 shape descriptors and
 features extracted from original and derived images (LoG with 5 sigma levels, 1 level of Wavelet decomposistions
 yielding 8 derived images and images derived using Square, Square Root, Logarithm and Exponential filters).
 

radiomics/glcm.py

@@ -74,6 +74,8 @@ class RadiomicsGLCM(base.RadiomicsFeaturesBase):
 
   The following class specific settings are possible:
 
+  - distances [[1]]: List of integers. This specifies the distances between the center voxel and the neighbor, for which
+    angles should be generated. See also :py:func:`~radiomics.imageoperations.generateAngles`
   - symmetricalGLCM [True]: boolean, indicates whether co-occurrences should be assessed in two directions per angle,
     which results in a symmetrical matrix, with equal distributions for :math:`i` and :math:`j`. A symmetrical matrix
     corresponds to the GLCM as defined by Haralick et al.

radiomics/ngtdm.py

@@ -0,0 +1,300 @@
+import numpy
+
+from radiomics import base, cMatrices, cMatsEnabled, imageoperations
+
+
+class RadiomicsNGTDM(base.RadiomicsFeaturesBase):
+  r"""
+  A Neighbouring Gray Tone Difference Matrix quantifies the difference between a gray value and the average gray value
+  of its neighbours within distance :math:`\delta`. The sum of absolute differences for gray level :math:`i` is stored in the matrix.
+  Let :math:`\textbf{X}_{gl}` be a set of segmented voxels and :math:`x_{gl}(j_x,j_y,j_z) \in \textbf{X}_{gl}` be the gray level of a voxel at postion
+  :math:`(j_x,j_y,j_z)`, then the average gray level of the neigbourhood is:
+
+  .. math::
+
+    \bar{A}_i &= \bar{A}(j_x, j_y, j_z) \\
+    &= \displaystyle\frac{1}{W} \displaystyle\sum_{k_x=-\delta}^{\delta}\displaystyle\sum_{k_y=-\delta}^{\delta}
+    \displaystyle\sum_{k_z=-\delta}^{\delta}{x_{gl}(j_x+k_x, j_y+k_y, j_z+k_z)}, \\
+    &\mbox{where }(k_x,k_y,k_z)\neq(0,0,0)\mbox{ and } x_{gl}(j_x+k_x, j_y+k_y, j_z+k_z) \in \textbf{X}_{gl}
+
+  Here, :math:`W` is the number of voxels in the neighbourhood that are also in :math:`\textbf{X}_{gl}`.
+
+  As a two dimensional example, let the following matrix :math:`\textbf{I}` represent a 4x4 image,
+  having 5 discrete grey levels, but no voxels with gray level 4:
+
+  .. math::
+    \textbf{I} = \begin{bmatrix}
+    1 & 2 & 5 & 2\\
+    3 & 5 & 1 & 3\\
+    1 & 3 & 5 & 5\\
+    3 & 1 & 1 & 1\end{bmatrix}
+
+  The following NGTDM can be obtained:
+
+  .. math::
+    \begin{array}{cccc}
+    i & n_i & p_i & s_i\\
+    \hline
+    1 & 6 & 0.375 & 13.35\\
+    2 & 2 & 0.125 & 2.00\\
+    3 & 4 & 0.25  & 2.63\\
+    4 & 0 & 0.00  & 0.00\\
+    5 & 4 & 0.25  & 10.075\end{array}
+
+  6 pixels have gray level 1, therefore:
+
+  :math:`s_1 = |1-10/3| + |1-30/8| + |1-15/5| + |1-13/5| + |1-15/5| + |1-11/3| = 13.35`
+
+  For gray level 2, there are 2 pixels, therefore:
+
+  :math:`s_2 = |2-15/5| + |2-15/5| = 2`
+
+  Similar for gray values 3 and 5:
+
+  :math:`s_3 = |3-14/5| + |3-18/5| + |3-20/8| + |3-5/3| = 2.63 \\
+  s_5 = |5-14/5| + |5-18/5| + |5-20/8| + |5-11/5| = 10.075`
+
+  Let:
+
+  :math:`n_i` be the number of voxels in :math:`X_{gl}` with gray level :math:`i`
+
+  :math:`N_v` be the total number of voxels in :math:`X_{gl}` and equal to :math:`\sum{n_i}`
+
+  :math:`p_i` be the gray level probability and equal to :math:`n_i/N_v`
+
+  :math:`s_i = \left\{ {\begin{array} {rcl}
+  \sum^{n_i}{|i-\bar{A}_i|} & \mbox{for} & n_i \neq 0 \\
+  0 & \mbox{for} & n_i = 0 \end{array}}\right.`
+  be the sum of absolute differences for gray level :math:`i`
+
+  :math:`N_g` be the number of discreet gray levels
+
+  :math:`N_{g,p}` be the number of gray levels where :math:`p_i \neq 0`
+
+  The following class specific settings are possible:
+
+  - distances [[1]]: List of integers. This specifies the distances between the center voxel and the neighbor, for which
+    angles should be generated. See also :py:func:`~radiomics.imageoperations.generateAngles`
+
+  References
+
+  - Amadasun M, King R; Textural features corresponding to textural properties;
+    Systems, Man and Cybernetics, IEEE Transactions on 19:1264-1274 (1989). doi: 10.1109/21.44046
+  """
+
+  def __init__(self, inputImage, inputMask, **kwargs):
+    super(RadiomicsNGTDM, self).__init__(inputImage, inputMask, **kwargs)
+
+    self.P_ngtdm = None
+
+    self._initSegmentBasedCalculation()
+
+  def _initSegmentBasedCalculation(self):
+    super(RadiomicsNGTDM, self)._initSegmentBasedCalculation()
+
+    self._applyBinning()
+
+    self.coefficients['Np'] = len(self.labelledVoxelCoordinates[0])
+    if cMatsEnabled:
+      self.P_ngtdm = self._calculateCMatrix()
+    else:
+      self.P_ngtdm = self._calculateMatrix()
+    self._calculateCoefficients()
+
+  def _calculateMatrix(self):
+    Ng = self.coefficients['Ng']
+
+    self.matrix = self.matrix.astype('float')
+
+    # Set voxels outside delineation to padding value
+    padVal = numpy.nan
+    self.matrix[(self.maskArray == 0)] = padVal
+
+    angles = imageoperations.generateAngles(self.boundingBoxSize, **self.kwargs)
+    angles = numpy.concatenate((angles, angles * -1))
+
+    dataTemp = numpy.zeros(self.matrix.shape, dtype='float')
+    countMat = numpy.zeros(self.matrix.shape, dtype='int')
+
+    with self.progressReporter(angles, desc='Calculate shifted matrices (NGTDM)') as bar:
+      for a in bar:
+        # 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 = 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:
+          angMat[:-a[0], :, :] = padVal
+
+        if a[1] > 0:
+          angMat[:, -a[1]:, :] = padVal
+        elif a[1] < 0:
+          angMat[:, :-a[1], :] = padVal
+
+        if a[2] > 0:
+          angMat[:, :, -a[2]:] = padVal
+        elif a[2] < 0:
+          angMat[:, :, :-a[2]] = padVal
+        nanmask = numpy.isnan(angMat)
+        dataTemp[~nanmask] += angMat[~nanmask]
+        countMat[~nanmask] += 1
+
+    # 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]
+
+    # Only consider voxels that are part of the Mask AND have a neighbourhood
+    validMask = nZeroMask * self.maskArray.astype('bool')
+    dataTemp[validMask] = numpy.abs(dataTemp[validMask] - self.matrix[validMask])  # 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)
+    grayLevels = self.coefficients['grayLevels']
+    with self.progressReporter(grayLevels, desc='Calculate NGTDM') as bar:
+      for i in bar:
+        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(dataTemp[i_ind])
+
+    # Fill in gray levels (needed as empty gray level slices will be deleted)
+    P_ngtdm[:, 2] = numpy.arange(1, Ng + 1)
+
+    # Delete empty grey levels
+    P_ngtdm = numpy.delete(P_ngtdm, numpy.where(P_ngtdm[:, 0] == 0), 0)
+
+    # Normalize P_ngtdm[:, 0] (= p_i)
+    P_ngtdm[:, 0] = P_ngtdm[:, 0] / self.coefficients['Np']
+
+    return P_ngtdm
+
+  def _calculateCMatrix(self):
+    angles = imageoperations.generateAngles(self.boundingBoxSize, **self.kwargs)
+    Ng = self.coefficients['Ng']
+
+    P_ngtdm = cMatrices.calculate_ngtdm(self.matrix, self.maskArray, angles, Ng)
+
+    # Delete empty grey levels
+    P_ngtdm = numpy.delete(P_ngtdm, numpy.where(P_ngtdm[:, 0] == 0), 0)
+
+    # Normalize P_ngtdm[:, 0] (= p_i)
+    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'] = self.P_ngtdm[:, 1]
+    self.coefficients['ivector'] = self.P_ngtdm[:, 2]
+
+    # Ngp = number of graylevels, for which p_i > 0
+    self.coefficients['Ngp'] = self.P_ngtdm.shape[0]
+
+  def getCoarsenessFeatureValue(self):
+    r"""
+    Calculate and return the coarseness.
+
+    :math:`Coarseness = \frac{1}{\sum^{N_g}_{i=1}{p_{i}s_{i}}}`
+
+    Coarseness is a measure of average difference between the center voxel and its neighbourhood and is an indication
+    of the spatial rate of change. A higher value indicates a lower spatial change rate and a locally more uniform texture.
+
+    N.B. :math:`\sum^{N_g}_{i=1}{p_{i}s_{i}}` potentially evaluates to 0 (in case of a completely homogeneous image).
+    If this is the case, an arbitrary value of :math:`10^6` is returned.
+    """
+    p_i = self.coefficients['p_i']
+    s_i = self.coefficients['s_i']
+    sum_coarse = numpy.sum(p_i * s_i)
+    if sum_coarse == 0:
+      return 1000000
+    else:
+      return 1 / sum_coarse
+
+  def getContrastFeatureValue(self):
+    r"""
+    Calculate and return the contrast.
+
+    :math:`Contrast = \left(\frac{1}{N_{g,p}(N_{g,p}-1)}\displaystyle\sum^{N_g}_{i=1}\displaystyle\sum^{N_g}_{j=1}{p_{i}p_{j}(i-j)^2}\right)
+    \left(\frac{1}{N_p^2}\displaystyle\sum^{N_g}_{i=1}{s_i}\right)\text{, where }p_i \neq 0, p_j \neq 0`
+
+    Contrast is a measure of the spatial intensity change, but is also dependent on the overall gray level dynamic range.
+    Contrast is high when both the dynamic range and the spatial change rate are high, i.e. an image with a large range
+    of gray levels, with large changes between voxels and their neighbourhood.
+    """
+    Ngp = self.coefficients['Ngp']
+    Np = self.coefficients['Np']
+    p_i = self.coefficients['p_i']
+    s_i = self.coefficients['s_i']
+    i = self.coefficients['ivector']
+    contrast = (numpy.sum(p_i[:, None] * p_i[None, :] * (i[:, None] - i[None, :]) ** 2) / (Ngp * (Ngp - 1))) * \
+               ((numpy.sum(s_i)) / (Np ** 2))
+
+    return contrast
+
+  def getBusynessFeatureValue(self):
+    r"""
+    Calculate and return the busyness.
+
+    :math:`Busyness = \frac{\sum^{N_g}_{i = 1}{p_{i}s_{i}}}{\sum^{N_g}_{i = 1}\sum^{N_g}_{j = 1}{|ip_i - jp_j|}}\text{, where }p_i \neq 0, p_j \neq 0`
+
+    A measure of the change from a pixel to its neighbour. A high value for busyness indicates a 'busy' image, with rapid
+    changes of intensity between pixels and its neighbourhood.
+
+    N.B. if :math:`N_{g,p} = 1`, then :math:`busyness = \frac{0}{0}`. If this is the case, 0 is returned, as it concerns
+    a fully homogeneous region.
+    """
+    p_i = self.coefficients['p_i']
+    s_i = self.coefficients['s_i']
+    i = self.coefficients['ivector']
+    absdiff = numpy.sum(numpy.abs((i * p_i)[:, None] - (i * p_i)[None, :]))
+    if absdiff == 0:
+      return 0
+    else:
+      busyness = numpy.sum(p_i * s_i) / absdiff
+    return busyness
+
+  def getComplexityFeatureValue(self):
+    r"""
+    Calculate and return the complexity.
+
+    :math:`Complexity = \frac{1}{N_p^2}\displaystyle\sum^{N_g}_{i = 1}\displaystyle\sum^{N_g}_{j = 1}{|i - j|
+    \frac{p_{i}s_{i} + p_{j}s_{j}}{p_i + p_j}}\text{, where }p_i \neq 0, p_j \neq 0`
+
+    An image is considered complex when there are many primitive components in the image, i.e. the image is non-uniform
+    and there are many rapid changes in gray level intensity.
+    """
+    Np = self.coefficients['Np']
+    p_i = self.coefficients['p_i']
+    s_i = self.coefficients['s_i']
+    i = self.coefficients['ivector']
+    complexity = numpy.sum(numpy.abs(i[:, None] - i[None, :]) * (
+      ((p_i * s_i)[:, None] + (p_i * s_i)[None, :]) / (p_i[:, None] + p_i[None, :]))) / Np ** 2
+    return complexity
+
+  def getStrengthFeatureValue(self):
+    r"""
+    Calculate and return the strength.
+
+    :math:`Strength = \frac{\sum^{N_g}_{i = 1}\sum^{N_g}_{j = 1}{(p_i + p_j)(i-j)^2}}{\sum^{N_g}_{i = 1}{s_i}}\text{, where }p_i \neq 0, p_j \neq 0`
+
+    Strenght is a measure of the primitives in an image. Its value is high when the primitives are easily defined and
+    visible, i.e. an image with slow change in intensity but more large coarse differences in gray level intensities.
+
+    N.B. :math:`\sum^{N_g}_{i=1}{s_i}` potentially evaluates to 0 (in case of a completely homogeneous image).
+    If this is the case, 0 is returned.
+    """
+    p_i = self.coefficients['p_i']
+    s_i = self.coefficients['s_i']
+    i = self.coefficients['ivector']
+    sum_s_i = numpy.sum(s_i)
+    if sum_s_i == 0:
+      return 0
+    else:
+      strength = numpy.sum((p_i[:, None] + p_i[None, :]) * (i[:, None] - i[None, :]) ** 2) / sum_s_i
+      return strength