Comparison Pyradiomics GLCM x Scikit- Image GLCM

[Guiandreis]

Hi @JoostJM, I'm new at image feature extraction so I'm exploring some open-source python packages available for it. Calculating GLCM using Pyradiomics and Scikit-Image packages provides different results (focusing on contrast and correlation) using the same image and method, could you help me understand if the different results are caused by my code? I think Scikit-Image package is more straight forward to understand but Pyradiomics provides way more final information. I'm using single slice (2D image) of the Brain1 data. Here's my code on Jupyter Notebook:

from radiomics import featureextractor,base,cMatrices,glrlm,glcm,firstorder,imageoperations,shape,glszm,shape2D,ngtdm,gldm
import radiomics
import SimpleITK as sitk
import six # six comunica python 2 e python 3
import sys
import numpy as np
import pandas as pd
from skimage.feature import greycomatrix, greycoprops
from zipfile import ZipFile
import matplotlib.pyplot as plt

imagename,maskname=radiomics.getTestCase('brain1')
imageradio=sitk.ReadImage(imagename)
mask2radio=sitk.ReadImage(maskname)
image1=sitk.GetArrayFromImage(imageradio)
mask22=sitk.GetArrayFromImage(mask2radio)
image12=image1[12,:,:]
mask12=mask22[12,:,:]
image12sitk=sitk.GetImageFromArray(image12)
mask12sitk=sitk.GetImageFromArray(mask12)

Pyradiomics GLCM
image12sitk=sitk.GetImageFromArray(image12)
mask12sitk=sitk.GetImageFromArray(mask12)
imagecalc=image12sitk
maskcalc=mask12sitk
label=1
settings={}
settings['binWidth']=25
settings['resampledPixelSpacing']=None
#settings[resampledPixelSpacing]=[3,3,3]
settings['interpolator']='stikBSpline'
settings['verbose']=True
settings['force2D']=True
settings['force2Ddimension']=0
#settings['distances']=[1]
#settings['voxelArrayShift']=300

interpolator=settings.get('interpolator')
resampledPixelSpacing=settings.get('resampledPixelSpacing')
if interpolator is not None and resampledPixelSpacing is not None:
    print('a')
    imagecalc,maskcalc=imageoperations.resampleImage(imagecalc,maskcalc,**settings)

#crop image bb is the bounding box, upon which the image and mask are corpped
bb,correctedMask=imageoperations.checkMask(imagecalc,maskcalc,label=1)
if correctedMask is not None:
    imagecalc=correctedMask
croppedImage,croppedMask=imageoperations.cropToTumorMask(imagecalc,maskcalc,bb)
print(bb)

firstOrderFeatures=firstorder.RadiomicsFirstOrder(croppedImage,croppedMask,**settings)
firstOrderFeatures.enableFeatureByName('Mean',True)

print('calculating first order features...',)
result=firstOrderFeatures.execute()
print('done')

print('calculated first order features:')
for (key,val) in six.iteritems(result):
    print('  ',key,": ",val)

shapeFeatures=shape2D.RadiomicsShape2D(croppedImage,croppedMask,**settings)
shapeFeatures.enableAllFeatures()

print('calculating shape features..',)
result=shapeFeatures.execute()
print('done')

print('Calculated shape features:')
for (key,val) in six.iteritems(result):
    print(' ',key,':',val)

# Calculate the features and print(out result)
glcmfeatures=glcm.RadiomicsGLCM(croppedImage,croppedMask,**settings)
glcmfeatures.enableAllFeatures()

print('Calculating GLCM features...',)
result = glcmfeatures.execute()
print('done')

print('Calculated GLCM features: ')
for (key, val) in six.iteritems(result):
    print('  ', key, ':', val)

Scikit-image GLCM
img12masked=image12[ 89: 147,185: 205]
greymatrixglcm=greycomatrix(img12masked,[1],[0,np.pi/4,np.pi,3*np.pi/4],levels=img12masked.max()+1,normed=True)
glcmcor=greycoprops(greymatrixglcm,prop='correlation')
glcmcont=greycoprops(greymatrixglcm,prop='contrast')
glcmdiss=greycoprops(greymatrixglcm,prop='dissimilarity')
glcmhomo=greycoprops(greymatrixglcm,prop='homogeneity')
glcmasm=greycoprops(greymatrixglcm,prop='ASM')
glcmene=greycoprops(greymatrixglcm,prop='energy')
print('correlation',glcmcor,'\ncontrast ',glcmcont,'\ndissimilarity ',glcmdiss,'\nhomogeneity',glcmhomo,'\nASM',glcmasm,'\nenergy',glcmene)
print('\ncorrelation mean',np.mean(glcmcor),'\ncontrast mean',np.mean(glcmcont),'\ndissimilarity mean ',np.mean(glcmdiss),'\nhomogeneity mean',np.mean(glcmhomo),'\nASM mean',np.mean(glcmasm),
      '\nenergy mean',np.mean(glcmene),'\nfirst order mean',np.mean(img12masked))

[JoostJM]

I'm not terribly familliar with scikit-image's calculation of GLCM, but have you looked into how scikit-image defines the neighbors? The fact that you have to input something like pi into the call makes me wonder if scikit image uses a euclidean norm for the distance function (resampling intensities not at right angles to the center voxel). PyRadiomics uses infinity norm, and does not resample voxels to get neighbor values.

in addition, are you sure you defined your angles correctly? [0,np.pi/4,np.pi,3*np.pi/4] probably overlaps angle 0 and np.pi (in degrees, this is [0, 45, 180, 135], my suggestion would be [0,np.pi/4,np.pi/2,3*np.pi/4] to also include angle 90.