diff --git a/bin/Notebooks/helloRadiomics.ipynb b/bin/Notebooks/helloRadiomics.ipynb index 92042d8a..9ee45847 100644 --- a/bin/Notebooks/helloRadiomics.ipynb +++ b/bin/Notebooks/helloRadiomics.ipynb @@ -230,12 +230,90 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Getting the docstrings of the active features" + "### Optional progress reporting using a progress bar\n", + "Calculation of the GLCM and GLSZM matrices in PyRadiomics can take some time when in full-python mode. Therefore, PyRadiomics provides a variable that can hold a progress bar class, which will be used when the verbosity is INFO or DEBUG.\n", + "\n", + "This class must support usage in a `with` statement (i.e. exposes functions `__enter__` and `__exit__`), and takes an iterable as first argument (and exposes an `__iter__` function to iterate over that iterable). A description for the bar (label) is passed in a keyword argument called 'desc'. Run the cell below to illustrate this using the `tqdm` package, which has the `tqdm` class, which fits these requirements.\n", + "\n", + "**N.B. This cell will only work if the 'tqdm' package is installed, which is not included in the requirements of PyRadiomics.**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "extractor.kwargs['enableCExtensions'] = False\n", + "\n", + "# Enable the GLCM class to show the progress bar\n", + "extractor.enableFeatureClassByName('glcm')\n", + "\n", + "radiomics.setVerbosity(logging.INFO) # Verbosity must be at least INFO to enable progress bar\n", + "\n", + "import tqdm\n", + "radiomics.progressReporter = tqdm.tqdm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `click` package also exposes a progress bar, which can also take an iterable as its first argument, but the label is passed in keyword 'label'. To use `click`, we need to create a simple wrapper that redirects the 'desc' keyword to the 'label' keyword.\n", + "\n", + "The cell below is an alternative to the cell showing how to use `tqdm` as progressbar.\n", + "\n", + "**N.B. This cell will only work if the 'click' package is installed, which is not included in the requirements of PyRadiomics.**" ] }, { "cell_type": "code", "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "extractor.kwargs['enableCExtensions'] = False\n", + "\n", + "\n", + "# Enable the GLCM class to show the progress bar\n", + "extractor.enableFeatureClassByName('glcm')\n", + "\n", + "radiomics.setVerbosity(logging.INFO) # Verbosity must be at least INFO to enable progress bar\n", + "\n", + "import click\n", + "\n", + "class progressWrapper():\n", + " def __init__(self, iterable, desc=''):\n", + " # For a click progressbar, the description must be provided in the 'label' keyword argument.\n", + " # The iterable can be passed as first argument to click, no changes needed\n", + " self.bar = click.progressbar(iterable, label=desc) \n", + "\n", + " def __iter__(self):\n", + " return self.bar.__iter__() # Redirect to the __iter__ function of the click progressbar\n", + "\n", + " def __enter__(self):\n", + " return self.bar.__enter__() # Redirect to the __enter__ function of the click progressbar\n", + "\n", + " def __exit__(self, exc_type, exc_value, tb):\n", + " return self.bar.__exit__(exc_type, exc_value, tb) # Redirect to the __exit__ function of the click progressbar\n", + "\n", + "radiomics.progressReporter = progressWrapper" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Getting the docstrings of the active features" + ] + }, + { + "cell_type": "code", + "execution_count": 9, "metadata": { "collapsed": false, "scrolled": true @@ -248,56 +326,74 @@ "Active features:\n", "10Percentile\n", "\n", - " Calculate the 10\\ :sup:`th` percentile in the image array.\n", + " **5. 10th percentile**\n", + "\n", + " The 10\\ :sup:`th` percentile of :math:`\\textbf{X}`\n", " \n", "90Percentile\n", "\n", - " Calculate the 90\\ :sup:`th` percentile in the image array.\n", + " **6. 90th percentile**\n", + "\n", + " The 90\\ :sup:`th` percentile of :math:`\\textbf{X}`\n", " \n", "Energy\n", "\n", - " Calculate the Energy of the image array.\n", + " **1. Energy**\n", "\n", - " :math:`energy = \\displaystyle\\sum^{N}_{i=1}{\\textbf{X}(i)^2}`\n", + " .. math::\n", "\n", - " Energy is a measure of the magnitude of voxel values in\n", - " an image. A larger values implies a greater sum of the\n", + " \\textit{energy} = \\displaystyle\\sum^{N}_{i=1}{(\\textbf{X}(i) + c)^2}\n", + "\n", + " Here, :math:`c` is optional value, defined by ``voxelArrayShift``, which shifts the intensities to prevent negative\n", + " values in :math:`\\textbf{X}`. This ensures that voxels with the lowest gray values contribute the least to Energy,\n", + " instead of voxels with gray level intensity closest to 0.\n", + "\n", + " Energy is a measure of the magnitude of voxel values in an image. A larger values implies a greater sum of the\n", " squares of these values.\n", + "\n", + " .. note::\n", + "\n", + " This feature is volume-confounded, a larger value of :math:`c` increases the effect of volume-confounding.\n", " \n", "Entropy\n", "\n", - " Calculate the Entropy of the image array.\n", + " **3. Entropy**\n", "\n", - " :math:`entropy = -\\displaystyle\\sum^{N_l}_{i=1}{p(i)\\log_2\\big(p(i)+\\epsilon\\big)}`\n", + " .. math::\n", + "\n", + " \\textit{entropy} = -\\displaystyle\\sum^{N_l}_{i=1}{p(i)\\log_2\\big(p(i)+\\epsilon\\big)}\n", "\n", " Here, :math:`\\epsilon` is an arbitrarily small positive number (:math:`\\approx 2.2\\times10^{-16}`).\n", - " Entropy specifies the uncertainty/randomness in the\n", - " image values. It measures the average amount of\n", - " information required to encode the image values.\n", + "\n", + " Entropy specifies the uncertainty/randomness in the image values. It measures the average amount of information\n", + " required to encode the image values.\n", " \n", "InterquartileRange\n", "\n", - " Calculate the interquartile range of the image array.\n", + " **10. Interquartile Range**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{interquartile range} = \\textbf{P}_{75} - \\textbf{P}_{25}\n", "\n", - " :math:`interquartile\\ range = \\textbf{P}_{75} - \\textbf{P}_{25}`, where :math:`\\textbf{P}_{25}` and\n", - " :math:`\\textbf{P}_{75}` are the 25\\ :sup:`th` and 75\\ :sup:`th` percentile of the image array, respectively.\n", + " Here :math:`\\textbf{P}_{25}` and :math:`\\textbf{P}_{75}` are the 25\\ :sup:`th` and 75\\ :sup:`th` percentile of the\n", + " image array, respectively.\n", " \n", "Kurtosis\n", "\n", - " Calculate the Kurtosis of the image array.\n", + " **17. Kurtosis**\n", "\n", - " :math:`kurtosis = \\displaystyle\\frac{\\mu_4}{\\sigma^4}\n", - " = \\frac{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^4}}\n", - " {\\left(\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X}})^2\\right)^2}`\n", + " .. math::\n", + "\n", + " \\textit{kurtosis} = \\displaystyle\\frac{\\mu_4}{\\sigma^4} =\n", + " \\frac{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^4}}\n", + " {\\left(\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X}})^2\\right)^2}\n", "\n", " Where :math:`\\mu_4` is the 4\\ :sup:`th` central moment.\n", "\n", - " Kurtosis is a measure of the 'peakedness' of the distribution\n", - " of values in the image ROI. A higher kurtosis implies that the\n", - " mass of the distribution is concentrated towards the tail(s)\n", - " rather than towards the mean. A lower kurtosis implies the reverse:\n", - " that the mass of the distribution is concentrated towards a\n", - " spike near the Mean value.\n", + " Kurtosis is a measure of the 'peakedness' of the distribution of values in the image ROI. A higher kurtosis implies\n", + " that the mass of the distribution is concentrated towards the tail(s) rather than towards the mean. A lower kurtosis\n", + " implies the reverse: that the mass of the distribution is concentrated towards a spike near the Mean value.\n", "\n", " Related links:\n", "\n", @@ -305,70 +401,101 @@ " \n", "Maximum\n", "\n", - " Calculate the Maximum Value in the image array.\n", + " **7. Maximum**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{maximum} = \\max(\\textbf{X})\n", + "\n", + " The maximum gray level intensity within the ROI.\n", " \n", "MeanAbsoluteDeviation\n", "\n", - " Calculate the Mean Absolute Deviation for the image array.\n", + " **12. Mean Absolute Deviation (MAD)**\n", "\n", - " :math:`mean\\ absolute\\ deviation = \\frac{1}{N}\\displaystyle\\sum^{N}_{i=1}{|\\textbf{X}(i)-\\bar{X}|}`\n", + " .. math::\n", "\n", - " Mean Absolute Deviation is the mean distance of all intensity values\n", - " from the Mean Value of the image array.\n", + " \\textit{MAD} = \\frac{1}{N}\\displaystyle\\sum^{N}_{i=1}{|\\textbf{X}(i)-\\bar{X}|}\n", + "\n", + " Mean Absolute Deviation is the mean distance of all intensity values from the Mean Value of the image array.\n", " \n", "Mean\n", "\n", - " Calculate the Mean Value for the image array.\n", + " **8. Mean**\n", + "\n", + " .. math::\n", "\n", - " :math:`mean = \\frac{1}{N}\\displaystyle\\sum^{N}_{i=1}{\\textbf{X}(i)}`\n", + " \\textit{mean} = \\frac{1}{N}\\displaystyle\\sum^{N}_{i=1}{\\textbf{X}(i)}\n", + "\n", + " The average gray level intensity within the ROI.\n", " \n", "Median\n", "\n", - " Calculate the Median Value for the image array.\n", + " **9. Median**\n", + "\n", + " The median gray level intensity within the ROI.\n", " \n", "Minimum\n", "\n", - " Calculate the Minimum Value in the image array.\n", + " **4. Minimum**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{minimum} = \\min(\\textbf{X})\n", " \n", "Range\n", "\n", - " Calculate the Range of Values in the image array.\n", + " **11. Range**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{range} = \\max(\\textbf{X}) - \\min(\\textbf{X})\n", "\n", - " :math:`range = \\max(\\textbf{X}) - \\min(\\textbf{X})`\n", + " The range of gray values in the ROI.\n", " \n", "RobustMeanAbsoluteDeviation\n", "\n", - " Calculate the Robust Mean Absolute Deviation for the image array.\n", + " **13. Robust Mean Absolute Deviation (rMAD)**\n", "\n", - " :math:`robust\\ mean\\ absolute\\ deviation = \\frac{1}{N_{10-90}}\\displaystyle\\sum^{N_{10-90}}_{i=1}{|\\textbf{X}_{10-90}(i)-\\bar{X}_{10-90}|}`\n", + " .. math::\n", + "\n", + " \\textit{rMAD} = \\frac{1}{N_{10-90}}\\displaystyle\\sum^{N_{10-90}}_{i=1}\n", + " {|\\textbf{X}_{10-90}(i)-\\bar{X}_{10-90}|}\n", "\n", " Robust Mean Absolute Deviation is the mean distance of all intensity values\n", " from the Mean Value calculated on the subset of image array with gray levels in between, or equal\n", - " to the 10\\ :sub:`th` and 90\\ :sub:`th` percentile.\n", + " to the 10\\ :sup:`th` and 90\\ :sup:`th` percentile.\n", " \n", "RootMeanSquared\n", "\n", - " Calculate the Root Mean Squared of the image array.\n", + " **14. Root Mean Squared (RMS)**\n", + "\n", + " .. math::\n", "\n", - " :math:`RMS = \\sqrt{\\frac{1}{N}\\sum^{N}_{i=1}{\\textbf{X}(i)^2}}`\n", + " \\textit{RMS} = \\sqrt{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i) + c)^2}}\n", "\n", - " RMS is the square-root of the mean of all the squared\n", - " intensity values. It is another measure of the magnitude\n", - " of the image values.\n", + " Here, :math:`c` is optional value, defined by ``voxelArrayShift``, which shifts the intensities to prevent negative\n", + " values in :math:`\\textbf{X}`. This ensures that voxels with the lowest gray values contribute the least to RMS,\n", + " instead of voxels with gray level intensity closest to 0.\n", + "\n", + " RMS is the square-root of the mean of all the squared intensity values. It is another measure of the magnitude of\n", + " the image values. This feature is volume-confounded, a larger value of :math:`c` increases the effect of\n", + " volume-confounding.\n", " \n", "Skewness\n", "\n", - " Calculate the Skewness of the image array.\n", + " **16. Skewness**\n", + "\n", + " .. math::\n", "\n", - " :math:`skewness = \\displaystyle\\frac{\\mu_3}{\\sigma^3}\n", - " = \\frac{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^3}}\n", - " {\\left(\\sqrt{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^2}}\\right)^3}`\n", + " \\textit{skewness} = \\displaystyle\\frac{\\mu_3}{\\sigma^3} =\n", + " \\frac{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^3}}\n", + " {\\left(\\sqrt{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^2}}\\right)^3}\n", "\n", " Where :math:`\\mu_3` is the 3\\ :sup:`rd` central moment.\n", "\n", - " Skewness measures the asymmetry of the distribution of values about the Mean value. Depending\n", - " on where the tail is elongated and the mass of the distribution\n", - " is concentrated, this value can be positive or negative.\n", + " Skewness measures the asymmetry of the distribution of values about the Mean value. Depending on where the tail is\n", + " elongated and the mass of the distribution is concentrated, this value can be positive or negative.\n", "\n", " Related links:\n", "\n", @@ -376,41 +503,376 @@ " \n", "StandardDeviation\n", "\n", - " Calculate the Standard Deviation of the image array.\n", + " **15. Standard Deviation**\n", "\n", - " :math:`standard\\ deviation = \\sqrt{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^2}}`\n", + " .. math::\n", "\n", - " Standard Deviation measures the amount of variation\n", - " or dispersion from the Mean Value.\n", + " \\textit{standard deviation} = \\sqrt{\\frac{1}{N}\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^2}}\n", + "\n", + " Standard Deviation measures the amount of variation or dispersion from the Mean Value. By definition,\n", + " :math:\\textit{standard deviation} = \\sqrt{\\textit{variance}}\n", " \n", "TotalEnergy\n", "\n", - " Calculate the Total Energy of the image array.\n", + " **2. Total Energy**\n", + "\n", + " .. math::\n", "\n", - " :math:`total\\ energy = V_{voxel}\\displaystyle\\sum^{N}_{i=1}{\\textbf{X}(i)^2}`\n", + " \\textit{total energy} = V_{voxel}\\displaystyle\\sum^{N}_{i=1}{(\\textbf{X}(i) + c)^2}\n", + "\n", + " Here, :math:`c` is optional value, defined by ``voxelArrayShift``, which shifts the intensities to prevent negative\n", + " values in :math:`\\textbf{X}`. This ensures that voxels with the lowest gray values contribute the least to Energy,\n", + " instead of voxels with gray level intensity closest to 0.\n", "\n", " Total Energy is the value of Energy feature scaled by the volume of the voxel in cubic mm.\n", + "\n", + " .. note::\n", + "\n", + " This feature is volume-confounded, a larger value of :math:`c` increases the effect of volume-confounding.\n", " \n", "Uniformity\n", "\n", - " Calculate the Uniformity of the image array.\n", + " **19. Uniformity**\n", + "\n", + " .. math::\n", "\n", - " :math:`uniformity = \\displaystyle\\sum^{N_l}_{i=1}{p(i)^2}`\n", + " \\textit{uniformity} = \\displaystyle\\sum^{N_l}_{i=1}{p(i)^2}\n", "\n", - " Uniformity is a measure of the sum of the squares of each intensity\n", - " value. This is a measure of the heterogeneity of the image array,\n", - " where a greater uniformity implies a greater heterogeneity or a\n", - " greater range of discrete intensity values.\n", + " Uniformity is a measure of the sum of the squares of each intensity value. This is a measure of the heterogeneity of\n", + " the image array, where a greater uniformity implies a greater heterogeneity or a greater range of discrete intensity\n", + " values.\n", " \n", "Variance\n", "\n", - " Calculate the Variance in the image array.\n", + " **18. Variance**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{variance} = \\frac{1}{N}\\displaystyle\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^2}\n", + "\n", + " Variance is the the mean of the squared distances of each intensity value from the Mean value. This is a measure of\n", + " the spread of the distribution about the mean. By definition, :math:`\\textit{variance} = \\sigma^2`\n", + " \n", + "Autocorrelation\n", + "\n", + " **1. Autocorrelation**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{autocorrelation} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{p(i,j)ij}\n", + "\n", + " Autocorrelation is a measure of the magnitude of the\n", + " fineness and coarseness of texture.\n", + " \n", + "AverageIntensity\n", + "\n", + " **2. Joint Averager**\n", + "\n", + " .. math::\n", + "\n", + " \\mu_x = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{p(i,j)i}\n", + "\n", + " Returns the mean gray level intensity of the :math:`i` distribution.\n", + "\n", + " .. warning::\n", "\n", - " :math:`variance = \\sigma^2 = \\frac{1}{N}\\displaystyle\\sum^{N}_{i=1}{(\\textbf{X}(i)-\\bar{X})^2}`\n", + " As this formula represents the average of the distribution of :math:`i`, it is independent from the\n", + " distribution of :math:`j`. Therefore, only use this formula if the GLCM is symmetrical, where\n", + " :math:`p_x(i) = p_y(j) \\text{, where } i = j`.\n", + " \n", + "ClusterProminence\n", + "\n", + " **3. Cluster Prominence**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{cluster prominence} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}\n", + " {\\big( i+j-\\mu_x(i)-\\mu_y(j)\\big)^4p(i,j)}\n", + "\n", + " Cluster Prominence is a measure of the skewness and asymmetry of the GLCM. A higher values implies more asymmetry\n", + " about the mean while a lower value indicates a peak near the mean value and less variation about the mean.\n", + " \n", + "ClusterShade\n", + "\n", + " **4. Cluster Shade**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{cluster shade} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}\n", + " {\\big(i+j-\\mu_x(i)-\\mu_y(j)\\big)^3p(i,j)}\n", + "\n", + " Cluster Shade is a measure of the skewness and uniformity of the GLCM.\n", + " A higher cluster shade implies greater asymmetry about the mean.\n", + " \n", + "ClusterTendency\n", + "\n", + " **5. Cluster Tendency**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{cluster tendency} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}\n", + " {\\big(i+j-\\mu_x(i)-\\mu_y(j)\\big)^2p(i,j)}\n", + "\n", + " Cluster Tendency is a measure of groupings of voxels with similar gray-level values.\n", + " \n", + "Contrast\n", "\n", - " Variance is the the mean of the squared distances of each intensity\n", - " value from the Mean value. This is a measure of the spread\n", - " of the distribution about the mean.\n", + " **6. Contrast**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{contrast} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{(i-j)^2p(i,j)}\n", + "\n", + " Contrast is a measure of the local intensity variation, favoring values away from the diagonal :math:`(i = j)`. A\n", + " larger value correlates with a greater disparity in intensity values among neighboring voxels.\n", + " \n", + "Correlation\n", + "\n", + " **7. Correlation**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{correlation} = \\frac{\\sum^{N_g}_{i=1}\\sum^{N_g}_{j=1}{p(i,j)ij-\\mu_x(i)\\mu_y(j)}}{\\sigma_x(i)\\sigma_y(j)}\n", + "\n", + " Correlation is a value between 0 (uncorrelated) and 1 (perfectly correlated) showing the\n", + " linear dependency of gray level values to their respective voxels in the GLCM.\n", + "\n", + " .. note::\n", + "\n", + " When there is only 1 discreet gray value in the ROI (flat region), :math:`\\sigma_x` and :math:`\\sigma_y` will be\n", + " 0. In this case, the value of correlation will be a NaN.\n", + " \n", + "DifferenceAverage\n", + "\n", + " **8. Difference Average**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{difference average} = \\displaystyle\\sum^{N_g-1}_{k=0}{kp_{x-y}(k)}\n", + "\n", + " Difference Average measures the relationship between occurrences of pairs\n", + " with similar intensity values and occurrences of pairs with differing intensity\n", + " values.\n", + " \n", + "DifferenceEntropy\n", + "\n", + " **9. Difference Entropy**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{difference entropy} = \\displaystyle\\sum^{N_g-1}_{k=0}{p_{x-y}(k)\\log_2\\big(p_{x-y}(k)+\\epsilon\\big)}\n", + "\n", + " Difference Entropy is a measure of the randomness/variability\n", + " in neighborhood intensity value differences.\n", + " \n", + "DifferenceVariance\n", + "\n", + " **10. Difference Variance**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{difference variance} = \\displaystyle\\sum^{N_g-1}_{k=0}{(1-DA)^2p_{x-y}(k)}\n", + "\n", + " Difference Variance is a measure of heterogeneity that places higher weights on\n", + " differing intensity level pairs that deviate more from the mean.\n", + " \n", + "Dissimilarity\n", + "\n", + " **11. Dissimilarity**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{dissimilarity} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{|i-j|p(i,j)}\n", + "\n", + " Dissimilarity is a measure of local intensity variation defined as the mean absolute difference between the\n", + " neighbouring pairs. A larger value correlates with a greater disparity in intensity values\n", + " among neighboring voxels.\n", + " \n", + "Energy\n", + "\n", + " **12. Joint Energy**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{energy} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{\\big(p(i,j)\\big)^2}\n", + "\n", + " Energy (or Angular Second Moment)is a measure of homogeneous patterns\n", + " in the image. A greater Energy implies that there are more instances\n", + " of intensity value pairs in the image that neighbor each other at\n", + " higher frequencies.\n", + " \n", + "Entropy\n", + "\n", + " **13. Joint Entropy**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{entropy} = -\\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}\n", + " {p(i,j)\\log_2\\big(p(i,j)+\\epsilon\\big)}\n", + "\n", + " Entropy is a measure of the randomness/variability in neighborhood intensity values.\n", + " \n", + "Homogeneity1\n", + "\n", + " **14. Homogeneity 1**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{homogeneity 1} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{\\frac{p(i,j)}{1+|i-j|}}\n", + "\n", + " Homogeneity 1 is a measure of the similarity in intensity values for\n", + " neighboring voxels. It is a measure of local homogeneity that increases\n", + " with less contrast in the window.\n", + " \n", + "Homogeneity2\n", + "\n", + " **15. Homogeneity 2**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{homogeneity 2} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{\\frac{p(i,j)}{1+|i-j|^2}}\n", + "\n", + " Homogeneity 2 is a measure of the similarity in intensity values\n", + " for neighboring voxels.\n", + " \n", + "Id\n", + "\n", + " **20. Inverse Difference (ID)**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{ID} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{ \\frac{p(i,j)}{1+|i-j|} }\n", + "\n", + " ID (inverse difference) is another measure of the local homogeneity of an image.\n", + " With more uniform gray levels, the denominator will remain low, resulting in a higher overall value.\n", + " \n", + "Idm\n", + "\n", + " **18. Inverse Difference Moment (IDM)**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{IDM} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{ \\frac{p(i,j)}{1+|i-j|^2} }\n", + "\n", + " IDM (inverse difference moment) is a measure of the local\n", + " homogeneity of an image. IDM weights are the inverse of the Contrast\n", + " weights (decreasing exponentially from the diagonal i=j in the GLCM).\n", + " \n", + "Idmn\n", + "\n", + " **19. Inverse Difference Moment Normalized (IDMN)**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{IDMN} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}\n", + " { \\frac{p(i,j)}{1+\\left(\\frac{|i-j|^2}{N_g^2}\\right)} }\n", + "\n", + " IDMN (inverse difference moment normalized) is a measure of the local\n", + " homogeneity of an image. IDMN weights are the inverse of the Contrast\n", + " weights (decreasing exponentially from the diagonal :math:`i=j` in the GLCM).\n", + " Unlike Homogeneity2, IDMN normalizes the square of the difference between\n", + " neighboring intensity values by dividing over the square of the total\n", + " number of discrete intensity values.\n", + " \n", + "Idn\n", + "\n", + " **21. Inverse Difference Normalized (IDN)**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{IDN} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}\n", + " { \\frac{p(i,j)}{1+\\left(\\frac{|i-j|}{N_g}\\right)} }\n", + "\n", + " IDN (inverse difference normalized) is another measure of the local\n", + " homogeneity of an image. Unlike Homogeneity1, IDN normalizes the difference\n", + " between the neighboring intensity values by dividing over the total number\n", + " of discrete intensity values.\n", + " \n", + "Imc1\n", + "\n", + " **16. Informal Measure of Correlation (IMC) 1**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{IMC 1} = \\frac{HXY-HXY1}{\\max\\{HX,HY\\}}\n", + " \n", + "Imc2\n", + "\n", + " **17. Informal Measure of Correlation (IMC) 2**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{IMC 2} = \\sqrt{1-e^{-2(HXY2-HXY)}}\n", + " \n", + "InverseVariance\n", + "\n", + " **22. Inverse Variance**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{inverse variance} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{\\frac{p(i,j)}{|i-j|^2}},\n", + " i \\neq j\n", + " \n", + "MaximumProbability\n", + "\n", + " **23. Maximum Probability**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{maximum probability} = \\max\\big(p(i,j)\\big)\n", + "\n", + " Maximum Probability is occurrences of the most predominant pair of\n", + " neighboring intensity values.\n", + " \n", + "SumAverage\n", + "\n", + " **24. Sum Average**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{sum average} = \\displaystyle\\sum^{2N_g}_{k=2}{p_{x+y}(k)k}\n", + "\n", + " Sum Average measures the relationship between occurrences of pairs\n", + " with lower intensity values and occurrences of pairs with higher intensity\n", + " values.\n", + " \n", + "SumEntropy\n", + "\n", + " **25. Sum Entropy**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{sum entropy} = \\displaystyle\\sum^{2N_g}_{k=2}{p_{x+y}(k)\\log_2\\big(p_{x+y}(k)+\\epsilon\\big)}\n", + "\n", + " Sum Entropy is a sum of neighborhood intensity value differences.\n", + " \n", + "SumSquares\n", + "\n", + " **27. Sum of Squares**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{sum squares} = \\displaystyle\\sum^{N_g}_{i=1}\\displaystyle\\sum^{N_g}_{j=1}{(i-\\mu_x)^2p(i,j)}\n", + "\n", + " Sum of Squares or Variance is a measure in the distribution of neigboring intensity level pairs\n", + " about the mean intensity level in the GLCM.\n", + "\n", + " .. warning::\n", + "\n", + " This formula represents the variance of the distribution of :math:`i` and is independent from the distribution\n", + " of :math:`j`. Therefore, only use this formula if the GLCM is symmetrical, where\n", + " :math:`p_x(i) = p_y(j) \\text{, where } i = j`\n", + " \n", + "SumVariance2\n", + "\n", + " **26. Sum Variance**\n", + "\n", + " .. math::\n", + "\n", + " \\textit{sum variance} = \\displaystyle\\sum^{2N_g}_{k=2}{(k-SA)^2p_{x+y}(k)}\n", + "\n", + " Sum Variance is a measure of heterogeneity that places higher weights on\n", + " neighboring intensity level pairs that deviate more from the mean.\n", " \n" ] } @@ -434,7 +896,42 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Calculating features\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Disabling C extensions\n", + "Calculating features with label: 1\n", + "Loading image and mask\n", + "Adding additional extraction information\n", + "Adding image type \"Original\" with settings: {'distances': [1], 'verbose': True, 'force2Ddimension': 0, 'enableCExtensions': False, 'force2D': False, 'interpolator': 'sitkBSpline', 'resampledPixelSpacing': None, 'normalizeScale': 1, 'normalize': False, 'additionalInfo': True, 'padDistance': 5, 'removeOutliers': None, 'minimumROISize': None, 'binWidth': 25, 'label': 1, 'minimumROIDimensions': 1}\n", + "Calculating features for original image\n", + "Computing firstorder\n", + "Computing glcm\n", + "calculate GLCM: 100%|██████████████████████████████████████████████████████████████████| 33/33 [00:00<00:00, 51.08it/s]\n" + ] + } + ], + "source": [ + "print(\"Calculating features\")\n", + "featureVector = extractor.execute(imageName, maskName)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, "metadata": { "collapsed": false }, @@ -443,26 +940,25 @@ "name": "stdout", "output_type": "stream", "text": [ - "Calculating features\n", "Computed general_info_BoundingBox: (162; 84; 11; 47; 70; 7)\n", - "Computed general_info_GeneralSettings: {'normalize': False; 'additionalInfo': True; 'verbose': True; 'removeOutliers': None; 'binWidth': 25; 'enableCExtensions': True; 'label': 1; 'interpolator': 'sitkBSpline'; 'resampledPixelSpacing': None; 'padDistance': 5; 'normalizeScale': 1}\n", + "Computed general_info_GeneralSettings: {'distances': [1]; 'verbose': True; 'additionalInfo': True; 'enableCExtensions': False; 'force2D': False; 'interpolator': 'sitkBSpline'; 'resampledPixelSpacing': None; 'label': 1; 'normalizeScale': 1; 'normalize': False; 'force2Ddimension': 0; 'removeOutliers': None; 'minimumROISize': None; 'binWidth': 25; 'minimumROIDimensions': 1; 'padDistance': 5}\n", "Computed general_info_ImageHash: 5c9ce3ca174f0f8324aa4d277e0fef82dc5ac566\n", "Computed general_info_ImageSpacing: (0.7812499999999999; 0.7812499999999999; 6.499999999999998)\n", "Computed general_info_InputImages: {'Original': {}}\n", "Computed general_info_MaskHash: 9dc2c3137b31fd872997d92c9a92d5178126d9d3\n", - "Computed general_info_Version: 1.1.0.post9.dev0+g409cdcd\n", + "Computed general_info_Version: 1.1.1.post5+g0d88779\n", "Computed general_info_VolumeNum: 2\n", "Computed general_info_VoxelNum: 4137\n", "Computed original_firstorder_InterquartileRange: 253.0\n", "Computed original_firstorder_Skewness: 0.275650859086\n", "Computed original_firstorder_Uniformity: 0.0451569635559\n", "Computed original_firstorder_MeanAbsoluteDeviation: 133.447261953\n", - "Computed original_firstorder_Energy: 33122817481.0\n", + "Computed original_firstorder_Energy: 2918821481.0\n", "Computed original_firstorder_RobustMeanAbsoluteDeviation: 103.00138343\n", "Computed original_firstorder_Median: 812.0\n", - "Computed original_firstorder_TotalEnergy: 131407662126.0\n", + "Computed original_firstorder_TotalEnergy: 11579797135.3\n", "Computed original_firstorder_Maximum: 1266.0\n", - "Computed original_firstorder_RootMeanSquared: 2829.57282108\n", + "Computed original_firstorder_RootMeanSquared: 839.964644818\n", "Computed original_firstorder_90Percentile: 1044.4\n", "Computed original_firstorder_Minimum: 468.0\n", "Computed original_firstorder_Entropy: 4.6019355539\n", @@ -471,26 +967,42 @@ "Computed original_firstorder_Variance: 24527.0792084\n", "Computed original_firstorder_10Percentile: 632.0\n", "Computed original_firstorder_Kurtosis: 2.18077293939\n", - "Computed original_firstorder_Mean: 825.235436307\n" + "Computed original_firstorder_Mean: 825.235436307\n", + "Computed original_glcm_Homogeneity1: 0.276140402104\n", + "Computed original_glcm_Homogeneity2: 0.189156155892\n", + "Computed original_glcm_ClusterShade: -52.9707943386\n", + "Computed original_glcm_MaximumProbability: 0.00792784235012\n", + "Computed original_glcm_SumVariance2: 103.142793792\n", + "Computed original_glcm_Idmn: 0.957796447609\n", + "Computed original_glcm_Contrast: 52.2310659277\n", + "Computed original_glcm_DifferenceEntropy: 3.79686113536\n", + "Computed original_glcm_InverseVariance: 0.188666637795\n", + "Computed original_glcm_Dissimilarity: 5.58932678922\n", + "Computed original_glcm_SumAverage: 33.4497492152\n", + "Computed original_glcm_DifferenceVariance: 17.6107741076\n", + "Computed original_glcm_Idn: 0.866370546902\n", + "Computed original_glcm_Idm: 0.189156155892\n", + "Computed original_glcm_Correlation: 0.335214788202\n", + "Computed original_glcm_Autocorrelation: 292.684050471\n", + "Computed original_glcm_SumEntropy: 5.31547876648\n", + "Computed original_glcm_AverageIntensity: 17.1242601309\n", + "Computed original_glcm_Energy: 0.00290880217681\n", + "Computed original_glcm_SumSquares: 39.9781084143\n", + "Computed original_glcm_ClusterProminence: 26251.1709801\n", + "Computed original_glcm_Entropy: 8.79428086119\n", + "Computed original_glcm_Imc2: 0.692033706271\n", + "Computed original_glcm_Imc1: -0.091940840043\n", + "Computed original_glcm_DifferenceAverage: 5.58932678922\n", + "Computed original_glcm_Id: 0.276140402104\n", + "Computed original_glcm_ClusterTendency: 103.142793792\n" ] } ], "source": [ - "print(\"Calculating features\")\n", - "featureVector = extractor.execute(imageName, maskName)\n", - "\n", + "# Show output\n", "for featureName in featureVector.keys():\n", " print(\"Computed %s: %s\" % (featureName, featureVector[featureName]))" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/bin/helloRadiomics.py b/bin/helloRadiomics.py index 31c0275b..0406bf38 100644 --- a/bin/helloRadiomics.py +++ b/bin/helloRadiomics.py @@ -10,6 +10,62 @@ import radiomics from radiomics import featureextractor +def tqdmProgressbar(): + """ + This function will setup the progress bar exposed by the 'tqdm' package. + Progress reporting is only used in PyRadiomics for the calculation of GLCM and GLSZM in full python mode, therefore + enable GLCM and full-python mode to show the progress bar functionality + + N.B. This function will only work if the 'click' package is installed (not included in the PyRadiomics requirements) + """ + global extractor + extractor.kwargs['enableCExtensions'] = False + # Enable the GLCM class to show the progress bar + extractor.enableFeatureClassByName('glcm') + + radiomics.setVerbosity(logging.INFO) # Verbosity must be at least INFO to enable progress bar + + import tqdm + radiomics.progressReporter = tqdm.tqdm + +def clickProgressbar(): + """ + This function will setup the progress bar exposed by the 'click' package. + Progress reporting is only used in PyRadiomics for the calculation of GLCM and GLSZM in full python mode, therefore + enable GLCM and full-python mode to show the progress bar functionality. + + Because the signature used to instantiate a click progress bar is different from what PyRadiomics expects, we need to + write a simple wrapper class to enable use of a click progress bar. In this case we only need to change the 'desc' + keyword argument to a 'label' keyword argument. + + N.B. This function will only work if the 'click' package is installed (not included in the PyRadiomics requirements) + """ + global extractor + + extractor.kwargs['enableCExtensions'] = False + # Enable the GLCM class to show the progress bar + extractor.enableFeatureClassByName('glcm') + + radiomics.setVerbosity(logging.INFO) # Verbosity must be at least INFO to enable progress bar + + import click + + class progressWrapper(): + def __init__(self, iterable, desc=''): + # For a click progressbar, the description must be provided in the 'label' keyword argument. + self.bar = click.progressbar(iterable, label=desc) + + def __iter__(self): + return self.bar.__iter__() # Redirect to the __iter__ function of the click progressbar + + def __enter__(self): + return self.bar.__enter__() # Redirect to the __enter__ function of the click progressbar + + def __exit__(self, exc_type, exc_value, tb): + return self.bar.__exit__(exc_type, exc_value, tb) # Redirect to the __exit__ function of the click progressbar + + radiomics.progressReporter = progressWrapper + testCase = 'brain1' dataDir = os.path.join(os.path.abspath(""), "..", "data") @@ -61,6 +117,12 @@ # Only enable mean and skewness in firstorder extractor.enableFeaturesByName(firstorder=['Mean', 'Skewness']) +# Uncomment on of these functions to show how PyRadiomics can use the 'tqdm' or 'click' package to report progress when +# running in full python mode. Assumes the respective package is installed (not included in the requirements) + +# tqdmProgressbar() +# clickProgressbar() + print("Active features:") for cls, features in six.iteritems(extractor.enabledFeatures): if len(features) == 0: diff --git a/docs/developers.rst b/docs/developers.rst index 94ed78d2..8e7c5be1 100644 --- a/docs/developers.rst +++ b/docs/developers.rst @@ -87,6 +87,62 @@ returned by multiple yield statements, or yield statements inside a loop. Please be returned on each call to yield and that ``inputImageName`` is a unique name for each returned derived image. Derived images must have the same dimensions and occupy the same physical space to ensure compatibility with the mask. +------------------ +Progress Reporting +------------------ + +When operating in full-python mode, the calculation of the texture matrices can take some time. Therefor PyRadiomics +provides the possibility to report the progress for calculation of GLCM and GLSZM. +This is only enabled in full-python mode when the verbosity (:py:func:`~radiomics.setVerbosity()`) is set to INFO or +DEBUG. By default, none is provided and no progress of matrix calculation will be reported. + +To enable progress reporting, the ``radiomics.progressReporter`` variable should be set to a class object (NOT an +instance), which fits the following signature: + +1. Accepts an iterable as the first positional argument and a keyword argument ('desc') specifying a label to display +2. Can be used in a 'with' statement (i.e. exposes a ``__enter__`` and ``__exit__`` function) +3. Is iterable (i.e. at least specifies an ``__iter__`` function, which iterates over the iterable passed at + initialization) + +It is also possible to create your own progress reporter. To achieve this, additionally specify a function ``__next__``, +and have the ``__iter__`` function return ``self``. The ``__next__`` function takes no arguments and returns a call to +the ``__next__`` function of the iterable (i.e. ``return self.iterable.__next__()``). Any prints/progress reporting +calls can then be inserted in this function prior to the return statement. + +In ``radiomics\__init__.py`` a dummy progress reporter (``_DummyProgressReporter``) is defined, which is used when +calculating in full-python mode, but progress reporting is not enabled (verbosity > INFO) or the ``progressReporter`` +variable is not set. + +To design a custom progress reporter, the following code can be adapted and used as progressReporter:: + + class MyProgressReporter(object): + def __init__(self, iterable, desc=''): + self.desc = desc # A description is which describes the progress that is reported + self.iterable = iterable # Iterable is required + + # This function identifies the class as iterable and should return an object which exposes + # the __next__ function that should be used to iterate over the object + def __iter__(self): + return self # return self to 'intercept' the calls to __next__ to insert reporting code. + + def __next__(self): + nextElement = self.iterable.__next__() + # Insert custom progress reporting code here. This is called for every iteration in the loop + # (once for each unique gray level in the ROI for GLCM and GLSZM) + + # By inserting after the call `self.iterable.__next__()` the function will exit before the + # custom code is run when the stopIteration error is raised. + return nextElement + + # This function is called when the 'with' statement is entered + def __enter__(self): + print (self.desc) # Print out the description upon start of the loop + return self # The __enter__ function should return itself + + # This function is called when the 'with' statement is exited + def __exit__(self, exc_type, exc_value, tb): + pass # If nothing needs to be closed or handled, so just specify 'pass' + ------------------------------ Addtional points for attention ------------------------------ diff --git a/radiomics/__init__.py b/radiomics/__init__.py index ce4deeda..6621aee9 100644 --- a/radiomics/__init__.py +++ b/radiomics/__init__.py @@ -161,6 +161,58 @@ def getInputImageTypes(): return _inputImages +class _DummyProgressReporter(object): + """ + This class represents the dummy Progress reporter and is used for where progress reporting is implemented, but not + enabled (when the progressReporter is not set or verbosity level > INFO). + + PyRadiomics expects that the _getProgressReporter function returns an object that takes an iterable and 'desc' keyword + argument at initialization. Furthermore, it should be iterable, where it iterates over the iterable passed at + initialization and it should be used in a 'with' statement. + + In this class, the __iter__ function redirects to the __iter__ function of the iterable passed at initialization. + The __enter__ and __exit__ functions enable usage in a 'with' statement + """ + def __init__(self, iterable, desc=''): + self.desc = desc # A description is not required, but is provided by PyRadiomics + self.iterable = iterable # Iterable is required + + def __iter__(self): + return self.iterable.__iter__() # Just iterate over the iterable passed at initialization + + def __enter__(self): + return self # The __enter__ function should return itself + + def __exit__(self, exc_type, exc_value, tb): + pass # Nothing needs to be closed or handled, so just specify 'pass' + + +def _getProgressReporter(*args, **kwargs): + """ + This function returns an instance of the progressReporter, or, if it is not set (None), returns a dummy progress + reporter. + + To enable progress reporting, the progressReporter variable should be set to a class object (NOT an instance), which + fits the following signature: + + 1. Accepts an iterable as the first positional argument and a keyword argument ('desc') specifying a label to display + 2. Can be used in a 'with' statement (i.e. exposes a __enter__ and __exit__ function) + 3. Is iterable (i.e. at least specifies an __iter__ function, which iterates over the iterable passed at + initialization). + + It is also possible to create your own progress reporter. To achieve this, additionally specify a function `__next__`, + and have the `__iter__` function return `self`. The `__next__` function takes no arguments and returns a call to the + `__next__` function of the iterable (i.e. `return self.iterable.__next__()`). Any prints/progress reporting calls can + then be inserted in this function prior to the return statement. + """ + global progressReporter + if progressReporter is None: + return _DummyProgressReporter(*args, **kwargs) + else: + return progressReporter(*args, **kwargs) + +progressReporter = None + debugging = True logger = logging.getLogger(__name__) logger.setLevel(logging.INFO) # Set default level of logger to INFO to reflect most common setting for a log file diff --git a/radiomics/base.py b/radiomics/base.py index d9b41c61..6cf4574e 100644 --- a/radiomics/base.py +++ b/radiomics/base.py @@ -41,7 +41,11 @@ class RadiomicsFeaturesBase(object): def __init__(self, inputImage, inputMask, **kwargs): self.logger = logging.getLogger(self.__module__) self.logger.debug('Initializing feature class') - self.verbose = radiomics.handler.level <= logging.INFO # check if the handler to stderr is set to INFO or lower + # check if the handler to stderr is set and its level is INFO or lower + if logging.NOTSET < radiomics.handler.level <= logging.INFO: + self.progressReporter = radiomics._getProgressReporter + else: + self.progressReporter = radiomics._DummyProgressReporter self.kwargs = kwargs self.binWidth = kwargs.get('binWidth', 25) diff --git a/radiomics/glcm.py b/radiomics/glcm.py index 2b6e363b..91da2dd8 100644 --- a/radiomics/glcm.py +++ b/radiomics/glcm.py @@ -1,6 +1,5 @@ import numpy from six.moves import range -from tqdm import trange from radiomics import base, cMatrices, cMatsEnabled, imageoperations @@ -151,31 +150,28 @@ def _calculateMatrix(self): P_glcm = numpy.zeros((Ng, Ng, int(angles.shape[0])), dtype='float64') - if self.verbose: bar = trange(Ng, desc='calculate GLCM') - - # iterate over gray levels for center voxel - for i in range(1, Ng + 1): - # give some progress - if self.verbose: bar.update() - - # get the indices to all voxels which have the current gray level i - i_indices = numpy.where(self.matrix == i) - - # iterate over gray levels for neighbouring voxel - for j in range(1, Ng + 1): - # get the indices to all voxels which have the current gray level j - j_indices = set(zip(*numpy.where(self.matrix == j))) - - for a_idx, a in enumerate(angles): - # get the corresponding indices of the neighbours for angle a - neighbour_indices = set(zip(*(i_indices + a[:, None]))) - - # The following intersection yields the indices to voxels with gray level j - # 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)) - P_glcm[i - 1, j - 1, a_idx] = count - if self.verbose: bar.close() + # If verbosity > INFO, or no progress reporter is set in radiomics.progressReporter, _dummyProgressReporter is used, + # which just iterates over the iterator without reporting progress + with self.progressReporter(range(1, Ng + 1), desc='calculate GLCM') as bar: + # iterate over gray levels for center voxel + for i in bar: + # get the indices to all voxels which have the current gray level i + i_indices = numpy.where(self.matrix == i) + + # iterate over gray levels for neighbouring voxel + for j in range(1, Ng + 1): + # get the indices to all voxels which have the current gray level j + j_indices = set(zip(*numpy.where(self.matrix == j))) + + for a_idx, a in enumerate(angles): + # get the corresponding indices of the neighbours for angle a + neighbour_indices = set(zip(*(i_indices + a[:, None]))) + + # The following intersection yields the indices to voxels with gray level j + # 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)) + P_glcm[i - 1, j - 1, a_idx] = count P_glcm = self._applyMatrixOptions(P_glcm, angles) diff --git a/radiomics/glszm.py b/radiomics/glszm.py index ab2e14d8..827a4d81 100644 --- a/radiomics/glszm.py +++ b/radiomics/glszm.py @@ -1,6 +1,5 @@ import numpy from six.moves import range -from tqdm import trange from radiomics import base, cMatrices, cMatsEnabled, imageoperations @@ -83,56 +82,52 @@ def _calculateMatrix(self): """ self.logger.debug('Calculating GLSZM matrix in Python') + Ng = self.coefficients['Ng'] + Np = self.coefficients['Np'] size = numpy.max(self.matrixCoordinates, 1) - numpy.min(self.matrixCoordinates, 1) + 1 angles = imageoperations.generateAngles(size, **self.kwargs) # Empty GLSZ matrix - P_glszm = numpy.zeros((self.coefficients['Ng'], self.coefficients['Np'])) - - # Iterate over all gray levels in the image - numGrayLevels = self.coefficients['Ng'] + 1 - - if self.verbose: bar = trange(numGrayLevels - 1, desc='calculate GLSZM') - - for i in range(1, numGrayLevels): - # give some progress - if self.verbose: bar.update() - - ind = zip(*numpy.where(self.matrix == i)) - ind = list(set(ind).intersection(set(zip(*self.matrixCoordinates)))) + P_glszm = numpy.zeros((Ng, Np)) - while ind: # check if ind is not empty: unprocessed regions for current gray level - # Pop first coordinate of an unprocessed zone, start new stack - ind_region = [ind.pop()] + # If verbosity > INFO, or no progress reporter is set in radiomics.progressReporter, _dummyProgressReporter is used, + # which just iterates over the iterator without reporting progress + with self.progressReporter(range(1, Ng + 1), desc='calculate GLSZM') as bar: + # Iterate over all gray levels in the image + for i in bar: + ind = zip(*numpy.where(self.matrix == i)) + ind = list(set(ind).intersection(set(zip(*self.matrixCoordinates)))) - # Define regionSize - regionSize = 0 + while ind: # check if ind is not empty: unprocessed regions for current gray level + # Pop first coordinate of an unprocessed zone, start new stack + ind_region = [ind.pop()] - # Grow zone for item popped from stack of region indices, loop until stack of region indices is exhausted - # Each loop represents one voxel belonging to current zone. Therefore, count number of loops as regionSize - while ind_region: - regionSize += 1 + # Define regionSize + regionSize = 0 - # Use pop to remove next node for set of unprocessed region indices - ind_node = ind_region.pop() + # Grow zone for item popped from stack of region indices, loop until stack of region indices is exhausted + # Each loop represents one voxel belonging to current zone. Therefore, count number of loops as regionSize + while ind_region: + regionSize += 1 - # get all coordinates in the 26-connected region, 2 voxels per angle - region_full = [tuple(sum(a) for a in zip(ind_node, angle_i)) for angle_i in angles] - region_full += [tuple(sum(a) for a in zip(ind_node, angle_i)) for angle_i in angles * -1] + # Use pop to remove next node for set of unprocessed region indices + ind_node = ind_region.pop() - # get all unprocessed coordinates in the 26-connected region with same gray level - region_level = list(set(ind).intersection(set(region_full))) + # get all coordinates in the 26-connected region, 2 voxels per angle + region_full = [tuple(sum(a) for a in zip(ind_node, angle_i)) for angle_i in angles] + region_full += [tuple(sum(a) for a in zip(ind_node, angle_i)) for angle_i in angles * -1] - # Remove already processed indices to prevent reprocessing - ind = list(set(ind) - set(region_level)) + # get all unprocessed coordinates in the 26-connected region with same gray level + region_level = list(set(ind).intersection(set(region_full))) - # Add all found neighbours to the total stack of unprocessed neighbours - ind_region.extend(region_level) + # Remove already processed indices to prevent reprocessing + ind = list(set(ind) - set(region_level)) - # Update the gray level size zone matrix - P_glszm[i - 1, regionSize - 1] += 1 + # Add all found neighbours to the total stack of unprocessed neighbours + ind_region.extend(region_level) - if self.verbose: bar.close() + # Update the gray level size zone matrix + P_glszm[i - 1, regionSize - 1] += 1 # 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 diff --git a/requirements.txt b/requirements.txt index 3e8318cc..73093bf3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,6 @@ numpy>=1.9.2 SimpleITK>=0.9.1 nose>=1.3.7 nose-parameterized>=0.5.0 -tqdm>=4.7.1 PyWavelets>=0.4.0 pykwalify>=1.6.0 six>=1.10.0