[PYTHON] Texture analysis learned with pyradiomics
pyradiomics
You can calculate such image features.
- First Order Statistics
- Shape-based (2D and 3D)
- Gray Level Cooccurence Matrix (GLCM)
- Gray Level Run Length Matrix (GLRLM)
- Gray Level Size Zone Matrix (GLSZM)
- Gray Level Dependece Matrix (GLDM)
- Neighboring Gray Tone Difference Matrix (NGTDM)
For images, one-dimensional quantitative value "first order" ~~ (but also includes features that are worried about whether it should be interpreted as one-dimensional) ~~, Shape that analyzes the shape such as how round or jagged, Pixel It is a library that can analyze GLCM etc. to check the arrangement pattern. Parts such as GLCM are written in C language for speeding up, and by referring to this as a Python object, the calculation speed is increased.
Ready to use
Since it is made by a researcher, it is very simple and can be immediately converted into a CSV table. Here, I will write up to the point where the calculation result is displayed.
environment
Google colaboratory
procedure
First, install the library temporarily.
!pip install pyradiomics
You can install it like this.
Collecting pyradiomics
Downloading https://files.pythonhosted.org/packages/1d/6b/797d3fd59e16675f6d8fa3ad0f778a025e06bcb1fd595439d43aff0a522e/pyradiomics-2.2.0-cp36-cp36m-manylinux1_x86_64.whl (159kB)
|████████████████████████████████| 163kB 2.8MB/s
Collecting pykwalify>=1.6.0
Downloading https://files.pythonhosted.org/packages/36/9f/612de8ca540bd24d604f544248c4c46e9db76f6ea5eb75fb4244da6ebbf0/pykwalify-1.7.0-py2.py3-none-any.whl (40kB)
|████████████████████████████████| 40kB 5.4MB/s
Requirement already satisfied: numpy>=1.9.2 in /usr/local/lib/python3.6/dist-packages (from pyradiomics) (1.17.3)
Collecting SimpleITK>=0.9.1
Downloading https://files.pythonhosted.org/packages/bb/06/f3a67ef0e108d18840fd5e83f831d5ef1710ba46f05465fc50f9a505b518/SimpleITK-1.2.3-cp36-cp36m-manylinux1_x86_64.whl (42.5MB)
|████████████████████████████████| 42.5MB 54kB/s
Requirement already satisfied: six>=1.10.0 in /usr/local/lib/python3.6/dist-packages (from pyradiomics) (1.12.0)
Collecting PyWavelets<=1.0.0,>=0.4.0
Downloading https://files.pythonhosted.org/packages/90/d0/09b2bf3368d5bba6ee1a8868ce94eebbb105fc8bf89fa43c90348b21a7cb/PyWavelets-1.0.0-cp36-cp36m-manylinux1_x86_64.whl (4.4MB)
|████████████████████████████████| 4.4MB 35.8MB/s
Requirement already satisfied: PyYAML>=3.11 in /usr/local/lib/python3.6/dist-packages (from pykwalify>=1.6.0->pyradiomics) (3.13)
Requirement already satisfied: docopt>=0.6.2 in /usr/local/lib/python3.6/dist-packages (from pykwalify>=1.6.0->pyradiomics) (0.6.2)
Requirement already satisfied: python-dateutil>=2.4.2 in /usr/local/lib/python3.6/dist-packages (from pykwalify>=1.6.0->pyradiomics) (2.6.1)
ERROR: albumentations 0.1.12 has requirement imgaug<0.2.7,>=0.2.5, but you'll have imgaug 0.2.9 which is incompatible.
Installing collected packages: pykwalify, SimpleITK, PyWavelets, pyradiomics
Found existing installation: PyWavelets 1.1.1
Uninstalling PyWavelets-1.1.1:
Successfully uninstalled PyWavelets-1.1.1
Successfully installed PyWavelets-1.0.0 SimpleITK-1.2.3 pykwalify-1.7.0 pyradiomics-2.2.0
(At the moment, I get an error saying that you should use the lower version of imgaug, but ignore it.)
Import the required functions.
import os # needed navigate the system to get the input data
from radiomics import featureextractor # This module is used for interaction with pyradiomics
import radiomics
from radiomics import firstorder, getTestCase, glcm, glrlm, glszm, imageoperations, shape, shape2D
import numpy as np
import six
#http://simpleitk.github.io/SimpleITK-Notebooks/01_Image_Basics.html
import SimpleITK as sitk
import matplotlib.pyplot as plt
We will use the test data provided by pyradiomics. The image format is the Nrrd format image that is often used in the 3D image processing area (3DSlicer). You can use the SimpleITK library for loading. It can support various image formats. Even if Nrrd is not available, you can substitute it with tiff or png. (If you use Dicom, you can convert it to NDArray with pydicom and then convert it as a SimpleITK Image.) (https://github.com/Radiomics/pyradiomics/issues/339) Let's use the "lung2" sample.
# brain1, brain2, breast1, lung1,lung2,(prostate_phantom)
# brain,1,2 may have the wrong mask (no label area)?
imageName, maskName = getTestCase('lung2')
image = sitk.ReadImage(imageName)
mask = sitk.ReadImage(maskName)
If you uploaded an image in a common format (eg tif) to your Colab notebook, do this:
#Call from within colab's File (example assuming it is in the parent directory).
# imagePath = os.path.join(os.getcwd(), "lung2" + "_image.tif")
# maskPath = os.path.join(os.getcwd(), "lung2" + "_label.tif")
# image = sitk.ReadImage(imagePath)
# mask = sitk.ReadImage(maskPath)
Make sure it loads properly. (Slicing can be as follows, but note that the pixel sequence is numpy.adarray (h, w) for itkImage (w, y))
ndImg = sitk.GetArrayFromImage(image)
ndLbl = sitk.GetArrayFromImage(mask)
plt.imshow(ndImg[24])
plt.show()
plt.imshow(ndLbl[24])
plt.show()
This data is 3D data. Because the 2D images are lined up in the Z-axis direction. I personally think that pyradiomics is supposed to perform 3D calculations from the beginning.
Create a configuration file. It's easy here, but it goes by default.
settings = {}
settings['binWidth'] = 25
# If enabled, resample image (resampled image is automatically cropped.
settings['resampledPixelSpacing'] = None # [3,3,3] is an example for defining resampling (voxels with size 3x3x3mm)
settings['interpolator'] = sitk.sitkBSpline
settings['label'] = 1 #Since the mask area has a pixel value of 1 (otherwise it is 0).
Crop the part of the mask. This will result in a 2D image set of only the tumor area. (It is a 3D volume.)
#Check the integrity of the mask and create a bbox# 2020/1/23 Addendum
bb, correctedMask = imageoperations.checkMask(image, mask)
if correctedMask is not None:
mask = correctedMask
image, mask = imageoperations.cropToTumorMask(image, mask, bb)
The preparation is complete. The rest is the calculation.
First, calculate the one-dimensional feature from the cropped volume. It also outputs the type and meaning.
firstOrderFeatures = firstorder.RadiomicsFirstOrder(image, mask, **settings)
# firstOrderFeatures.enableFeatureByName('Mean', True)
firstOrderFeatures.enableAllFeatures()
print('Will calculate the following first order features: ')
for f in firstOrderFeatures.enabledFeatures.keys():
print(' ', f)
print(getattr(firstOrderFeatures, 'get%sFeatureValue' % f).__doc__)
print('Calculating first order features...')
results = firstOrderFeatures.execute()
print('done')
print('Calculated first order features: ')
for (key, val) in six.iteritems(results):
print(' ', key, ':', val)
Next are the characteristics of the shape.
shapeFeatures = shape.RadiomicsShape(image, mask, **settings)
shapeFeatures.enableAllFeatures()
print('Will calculate the following Shape features: ')
for f in shapeFeatures.enabledFeatures.keys():
print(' ', f)
print(getattr(shapeFeatures, 'get%sFeatureValue' % f).__doc__)
print('Calculating Shape features...')
results = shapeFeatures.execute()
print('done')
print('Calculated Shape features: ')
for (key, val) in six.iteritems(results):
print(' ', key, ':', val)
From here, we will calculate features such as GLCM.
# Show GLCM features
glcmFeatures = glcm.RadiomicsGLCM(image, mask, **settings)
glcmFeatures.enableAllFeatures()
print('Will calculate the following GLCM features: ')
for f in glcmFeatures.enabledFeatures.keys():
print(' ', f)
print(getattr(glcmFeatures, 'get%sFeatureValue' % f).__doc__)
print('Calculating GLCM features...')
results = glcmFeatures.execute()
print('done')
print('Calculated GLCM features: ')
for (key, val) in six.iteritems(results):
print(' ', key, ':', val)
# Show GLRLM features
glrlmFeatures = glrlm.RadiomicsGLRLM(image, mask, **settings)
glrlmFeatures.enableAllFeatures()
print('Will calculate the following GLRLM features: ')
for f in glrlmFeatures.enabledFeatures.keys():
print(' ', f)
print(getattr(glrlmFeatures, 'get%sFeatureValue' % f).__doc__)
print('Calculating GLRLM features...')
results = glrlmFeatures.execute()
print('done')
print('Calculated GLRLM features: ')
for (key, val) in six.iteritems(results):
print(' ', key, ':', val)
# Show GLSZM features
glszmFeatures = glszm.RadiomicsGLSZM(image, mask, **settings)
glszmFeatures.enableAllFeatures()
print('Will calculate the following GLSZM features: ')
for f in glszmFeatures.enabledFeatures.keys():
print(' ', f)
print(getattr(glszmFeatures, 'get%sFeatureValue' % f).__doc__)
print('Calculating GLSZM features...')
results = glszmFeatures.execute()
print('done')
print('Calculated GLSZM features: ')
for (key, val) in six.iteritems(results):
print(' ', key, ':', val)
#I will make it around here.
The above process can also be a 2D analysis if you input one 2D image as input, but if you know that it is a 2D analysis from the beginning, you can input it as follows. Only pass one image.
# choose 2D array
#The memory is referenced and changed, so call it again.
imageName, maskName = getTestCase('lung2')
imageStack = sitk.ReadImage(imageName)
maskSatck = sitk.ReadImage(maskName)
#Extract pixels as adarray for parametric image calculations, etc.(Example for the 24th slice)
ndimage = sitk.GetArrayFromImage(imageStack)[24]
ndmask = sitk.GetArrayFromImage(maskStack)[24]
#Basically, the data that was originally a DICOM image has a dimension of 3, so it is expanded.
ndimage = np.expand_dims(ndimage, axis=0)
ndmask = np.expand_dims(ndmask, axis=0)
print(ndimage.shape)# (1, 512, 512)Even with force2D, the shape must be 3D.
print(ndmask.shape)# (1, 512, 512)
Do this if you want to reassign the pixel ndarray as the sitk image after calculating for various reasons.
image = sitk.GetImageFromArray(ndimage)
mask = sitk.GetImageFromArray(ndmask)
Reset information such as pixel size. (Added on 2019/11/14, added on 12/21)
# copy geometric/calibration info
image.CopyInformation(imageStack[:,:,24])
mask.CopyInformation(maskStack[:,:,24])
#Alternatively, if the original image has metadata and there is no change in geometric information such as matrix size, origin, pixel spacing, etc., copy the metadata from the original data.
#However, in the case of DICOM (though this time it is tiff), it is better to update the information such as series UID and image creation time.
for i,key in enumerate(imageStack[:,:,24].GetMetaDataKeys()):
image.SetMetaData(key, imageStack[:,:,24].GetMetaData(key))
for i,key in enumerate(maskStack[:,:,24].GetMetaDataKeys()):
mask.SetMetaData(key, maskStack[:,:,24].GetMetaData(key))
In the settings, add force2D.
settings = {}
settings['binWidth'] = 25
# If enabled, resample image (resampled image is automatically cropped.
settings['resampledPixelSpacing'] = None # [3,3,3] is an example for defining resampling (voxels with size 3x3x3mm)
settings['interpolator'] = sitk.sitkBSpline
settings['label'] = 1
settings['force2D'] = True
settings['correctMask'] = True #If the geometric information such as the origin of the image to be analyzed and the mask image are different, it will be adjusted automatically.
The rest of the process is the same, so only the first order is shown here.
# Show the first order feature calculations
firstOrderFeatures = firstorder.RadiomicsFirstOrder(image, mask, **settings)
# firstOrderFeatures.enableFeatureByName('Mean', True) #When specifying a specific item
firstOrderFeatures.enableAllFeatures()
#Comment out because it is a duplicate.
# print('Will calculate the following first order features: ')
# for f in firstOrderFeatures.enabledFeatures.keys():
# print(' ', f)
# print(getattr(firstOrderFeatures, 'get%sFeatureValue' % f).__doc__)
print('Calculating first order features...')
results = firstOrderFeatures.execute()
print('done')
print('Calculated first order features: ')
for (key, val) in six.iteritems(results):
print(' ', key, ':', val)
that's all.
VisionaryImagingServices, Inc. Tatsuaki Kobayashi
Reference
- https://github.com/Radiomics/pyradiomics/tree/master/radiomics
- https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6382170/
Recommended Posts