back to lib

Function glcmdesc

Synopse

The glcmdesc function computes descriptors from the gray level co-occurence matrix.

  • glcm,v = glcmdesc(f,offset=[],mask=[])
    • Output
      • v: array with GLCM descriptors
      • glcm: GLCM matrix
    • Input
      • f: input image
      • offset: defines the distance and orientation
      • mask: defines the mask

2D offsets

  • offset = [0 D] -> zero degrees, distance D
  • offset = [D 0] -> 90 degrees, distance D
  • offset = [D D] -> 135 degrees, distance D
  • offset = [D -D] -> 45 degrees, distance D

3D offsets

-standard 2D directions:

[0 d] 0 degrees
[d -d] 45 degrees
[d 0] 90 degrees
[d d] 135 degrees
  • additional 9 directions
[d 0 -d] 0 degrees 45 degrees
[d 0 0] straight up  
[d 0 d] 0 degree 135 degrees
[d d 0] 90 degrees 45 degrees
[d -d 0] 90 degrees 135 degrees
[d d -d] 45 degrees 45 degrees
[d -d d] 45 degrees 135 degrees
[d d d] 135 degrees 45 degrees
[d -d -d] 135 degrees 135 degrees

Description

The glcmdesc function extracts from the co-ocorrence matrix a set of descriptors.

Function Code

01. def glcmdesc(f,offset=[],mask=[]):
02.     import numpy as np
03.     from iatexture.graycomatrix import glcm
04.     lev = int(f.max()+1) # levels
05.     if mask != []: # if there is a mask, it is necessary to set a invalid value to all the pixels that are not in the mask (this invalid value can be any value that it is not present inside the mask)
06.         aux = f.copy()
07.         hist,_ = np.histogram(np.ravel(f[mask]),bins=lev,range=(0,lev-1)) # look for an empty level inside the mask
08.         if np.argwhere(hist==0).shape[0]!=0: # if there is an empty levels
09.             new_v = np.argwhere(hist==0).max() # it is selected the highest empty gray level
10.         else: # # if there is not an empty levels
11.             new_v = lev # the chosen invalid value is f.max()+1
12.             lev+=1 # in this case, the number of gray levels increases
13. 
14.         aux[~mask] = new_v # all pixels that are not inside the mask receive this invalid value
15.         g = glcm(lev,offset,aux) # the glcm matrix is computed
16.         g[new_v,:] = 0 # the rows and columns that corresponds to the invalid value (not inside the mask) are nulled
17.         g[:,new_v] = 0 # it is important to notice that null rows or columns do not change the descriptors
18.     else: # if there is no mask
19.         g = glcm(lev,offset,f)
20. 
21.     ### glcm normalization ###
22.     if g.sum() != 0:
23.         g = g.astype(float)/g.sum()
24. 
25.     ### compute auxiliary variables ###
26.     (num_level, num_level2) = g.shape
27.     I, J = np.ogrid[0:num_level, 0:num_level]
28.     I = 1+ np.array(range(num_level)).reshape((num_level, 1))
29.     J = 1+ np.array(range(num_level)).reshape((1, num_level))
30.     diff_i = I - np.apply_over_axes(np.sum, (I * g), axes=(0, 1))[0, 0]
31.     diff_j = J - np.apply_over_axes(np.sum, (J * g), axes=(0, 1))[0, 0]
32.     std_i = np.sqrt(np.apply_over_axes(np.sum, (g * (diff_i) ** 2),axes=(0, 1))[0, 0])
33.     std_j = np.sqrt(np.apply_over_axes(np.sum, (g * (diff_j) ** 2),axes=(0, 1))[0, 0])
34.     cov = np.apply_over_axes(np.sum, (g * (diff_i * diff_j)),axes=(0, 1))[0, 0]
35. 
36.     gxy = np.zeros(2*g.shape[0]-1)   ### g x+y
37.     gx_y = np.zeros(g.shape[0])  ### g x-y
38.     for i in xrange(g.shape[0]):
39.         for j in xrange(g.shape[0]):
40.             gxy[i+j] += g[i,j]
41.             gx_y[np.abs(i-j)] += g[i,j]
42.     mx_y = (gx_y*np.arange(len(gx_y))).sum()
43.     v = np.zeros(11)
44.     i,j = np.indices(g.shape)+1
45.     ii = np.arange(len(gxy))+2
46.     ii_ = np.arange(len(gx_y))
47. 
48. 
49.     ### compute descriptors ###
50.     v[0] = np.apply_over_axes(np.sum, (g ** 2), axes=(0, 1))[0, 0] # Angular second moment
51.     v[1] = np.apply_over_axes(np.sum, (g * ((I - J) ** 2)), axes=(0, 1))[0, 0] # Contrast
52.     if std_i>1e-15 and std_j>1e-15: # handle the special case of standard deviations near zero
53.         v[2] = cov/(std_i*std_j)#v[2] = greycoprops(g,'correlation') # Correlation
54.     else:
55.         v[2] = 1
56.     v[3] = np.apply_over_axes(np.sum, (g* (diff_i) ** 2),axes=(0, 1))[0, 0]# Sum of squares
57.     v[4] = np.sum(g * (1. / (1. + (I - J) ** 2))) # Inverse difference moment
58.     v[5] = (gxy*ii).sum() # Sum average
59.     v[6] = ((ii-v[5])*(ii-v[5])*gxy).sum() # Sum variance
60.     v[7] = -1*(gxy*np.log10(gxy+ np.spacing(1))).sum() # Sum entropy
61.     v[8] = -1*(g*np.log10(g+np.spacing(1))).sum() # Entropy
62.     v[9] = ((ii_-mx_y)*(ii_-mx_y)*gx_y).sum() # Difference variance
63.     v[10] = -1*(gx_y*np.log10(gx_y++np.spacing(1))).sum() # Difference entropy
64.     return g,v

2D Examples

Numerical example:

1. f = array( [[2,1,1,0,0,0,0,0,0],
2.             [1,2,0,0,0,2,0,1,0],
3.             [1,0,2,1,0,2,0,1,0],
4.             [0,0,0,2,0,1,1,0,0]], dtype=uint8)
1. from glcmdesc import glcmdesc
2. import numpy as np
3. 
4. glcm,v = glcmdesc(f)
5. print 'glcm: ', glcm
6. print 'glcm descriptors',v
glcm:  [[ 0.3125   0.09375  0.125  ]
 [ 0.1875   0.0625   0.03125]
 [ 0.125    0.0625   0.     ]]
glcm descriptors [ 0.18164062  1.375      -0.16984388  0.60058594  0.6125      3.1875
  0.96484375  0.56703915  0.81387275  0.609375    0.46999155]

Numerical example with a different orientation:

1. f = array( [[2,1,1,0,0],
2.             [1,2,0,0,0],
3.             [1,0,2,1,0],
4.             [0,0,0,2,0]], dtype=uint8)
1. from glcmdesc import glcmdesc
2. 
3. offset = np.array([1,1],dtype=np.int32)
4. glcm,v = glcmdesc(f,offset)
5. print 'glcm: ', glcm
6. print 'glcm descriptors',v
glcm:  [[ 0.25        0.08333333  0.        ]
 [ 0.41666667  0.          0.        ]
 [ 0.          0.          0.25      ]]
glcm descriptors [ 0.30555556  0.5         0.71095787  0.57638889  0.75        3.5         2.25
  0.45154499  0.54938312  0.25        0.30103   ]

Examples with mask

1. f  = array( [[1,2,2,2],
2.              [2,2,1,1],
3.              [1,2,0,1],
4.              [2,1,0,1]], dtype=uint8)
5. glcm,v = glcmdesc(f,mask = f!=1)
6. print 'glcm: ', glcm
7. print 'glcm descriptors',v
glcm:  [[ 0.    0.    0.  ]
 [ 0.    0.    0.  ]
 [ 0.25  0.    0.75]]
glcm descriptors [ 0.625       1.          1.          0.          0.8         5.5         0.75
  0.24421905  0.24421905  0.75        0.24421905]

Comparison between matrices and descriptors using mask or using no mask

1. f  = array( [[2,2,2,1],
2.              [2,2,0,1],
3.              [0,2,0,1],
4.              [1,1,1,1]], dtype=uint8)
5. glcm,v = glcmdesc(f,mask = f!=1)
6. print 'glcm: ', glcm
7. print 'glcm descriptors',v
glcm:  [[ 0.          0.          0.16666667]
 [ 0.          0.          0.        ]
 [ 0.33333333  0.          0.5       ]]
glcm descriptors [ 0.38888889  2.         -0.31622777  0.55555556  0.6         5.          1.
  0.30103     0.43924729  1.          0.30103   ]
1. f  = array( [[2,2,2],
2.              [2,2,0],
3.              [0,2,0]], dtype=uint8)
4. glcm,v = glcmdesc(f)
5. print 'glcm: \n', glcm
6. print 'glcm descriptors',v
glcm: 
[[ 0.          0.          0.16666667]
 [ 0.          0.          0.        ]
 [ 0.33333333  0.          0.5       ]]
glcm descriptors [ 0.38888889  2.         -0.31622777  0.55555556  0.6         5.          1.
  0.30103     0.43924729  1.          0.30103   ]

3D examples

1. import numpy as np
2. 
3. f = np.array([[[1,2,2,0,0,1],[0,0,1,2,2,1],[1,1,0,0,0,2],[1,1,1,2,2,2],[1,1,2,2,0,0]],
4.        [[1,2,2,0,0,1],[0,0,1,2,2,1],[1,1,0,0,0,2],[1,1,1,2,2,2],[1,1,2,2,0,0]],
5.        [[1,2,2,0,0,1],[0,0,1,2,2,1],[1,1,0,0,0,2],[1,1,1,2,2,2],[1,1,2,2,0,0]]], dtype=np.uint8)
6. print 'input array \n',f
7. glcm,v = glcmdesc(f) # orientation 0 degrees
8. print 'glcm: \n', glcm
9. print 'glcm descriptors',v
input array 
[[[1 2 2 0 0 1]
  [0 0 1 2 2 1]
  [1 1 0 0 0 2]
  [1 1 1 2 2 2]
  [1 1 2 2 0 0]]

 [[1 2 2 0 0 1]
  [0 0 1 2 2 1]
  [1 1 0 0 0 2]
  [1 1 1 2 2 2]
  [1 1 2 2 0 0]]

 [[1 2 2 0 0 1]
  [0 0 1 2 2 1]
  [1 1 0 0 0 2]
  [1 1 1 2 2 2]
  [1 1 2 2 0 0]]]
glcm: 
[[ 0.2   0.08  0.04]
 [ 0.04  0.16  0.16]
 [ 0.08  0.04  0.2 ]]
glcm descriptors [ 0.1488      0.8         0.41432451  0.64        0.744       4.08        1.9136
  0.684676    0.87752801  0.4864      0.40986496]

Different orientation

1. glcm,v = glcmdesc(f,offset = np.array([0,1,-1], dtype = np.int32)) # orientation 45 degrees
2. print 'glcm: \n', glcm
3. print 'glcm descriptors',v
glcm: 
[[ 0.    0.2   0.1 ]
 [ 0.05  0.2   0.05]
 [ 0.25  0.    0.15]]
glcm descriptors [ 0.18        1.7        -0.31083494  0.69        0.57        4.1         0.89
  0.48195333  0.78379231  0.7         0.47601599]
1. glcm,v = glcmdesc(f,offset = np.array([1,0,0], dtype = np.int32)) # straigh up
2. print 'glcm: \n', glcm
3. print 'glcm descriptors',v
glcm: 
[[ 0.3         0.          0.        ]
 [ 0.          0.36666667  0.        ]
 [ 0.          0.          0.33333333]]
glcm descriptors [  3.35555556e-01   0.00000000e+00   1.00000000e+00   6.32222222e-01
   1.00000000e+00   4.06666667e+00   2.52888889e+00   4.75671184e-01
   4.75671184e-01   0.00000000e+00  -9.64327467e-17]
1. glcm,v = glcmdesc(f,offset = np.array([3,0,-3], dtype = np.int32)) # straigh up
2. print 'glcm: \n', glcm
3. print 'glcm descriptors',v
glcm: 
[[0 0 0]
 [0 0 0]
 [0 0 0]]
glcm descriptors [ 0.  0.  1.  0.  0.  0.  0. -0. -0.  0. -0.]

Using mask

1. glcm,v = glcmdesc(f,offset = np.array([0,1,-1], dtype = np.int32),mask=f!=1) # orientation 45 degrees
2. print 'glcm: \n', glcm
3. print 'glcm descriptors',v
glcm: 
[[ 0.   0.   0.2]
 [ 0.   0.   0. ]
 [ 0.5  0.   0.3]]
glcm descriptors [ 0.38        2.8        -0.5         0.64        0.44        4.6         0.84
  0.265295    0.44717262  0.84        0.265295  ]
1. glcm,v = glcmdesc(f,offset = np.array([1,0,0], dtype = np.int32),mask=f!=1) # straigh up
2. print 'glcm: \n', glcm
3. print 'glcm descriptors',v
glcm: 
[[ 0.47368421  0.          0.        ]
 [ 0.          0.          0.        ]
 [ 0.          0.          0.52631579]]
glcm descriptors [  5.01385042e-01   0.00000000e+00   1.00000000e+00   9.97229917e-01
   1.00000000e+00   4.10526316e+00   3.98891967e+00   3.00428202e-01
   3.00428202e-01   0.00000000e+00  -9.64327467e-17]

Equation

The co-occurence matrix becomes the estimate of the joint probability, , of two pixels, a distance apart along a given direction having particular (co-occurring) values and .

Let , e , denote the mean and standard deviations of the row and column sums of the co-occurrence matrix, respectively [related to the marginal distributions and ], and the number of discrete intensity levels.

Angular second moment

Contrast

Correlation

Sum of squares

Inverse difference moment

Sum average

Sum variance

Sum entropy

Entropy

Difference variance

Difference entropy

Contributions

  • Mariana Bento, october, 23th 2013: initial function.
  • Mariana Leite, september, 9th 2014: C implementation and inclusion of 3D images.

Memory test

01. import multiprocessing as mp
02. import numpy as np
03. import psutil
04. from ia636 import iacontour,iameshgrid,iaroi,iagshow,ianormalize
05. import glob, urllib, os
06. from glcmdesc import glcmdesc
07. 
08. foldername = '/awmedia/www/media/p/LesionMRI/CAINdata/Original_Data/CAIN10320013/BASELINE - CAIN1/FLAIR_LongTR_601/*'
09. folder = glob.glob(foldername)
10. folder.sort()
11. 
12. img = zeros((len(folder),560,560),dtype=int32)
13. 
14. for i in arange(len(folder)):
15.      aux = dicom.read_file(find_attachment_file(folder[i]))
16.      aux = aux.pixel_array.astype(int32)
17.      img[i] = ianormalize(aux)
18. 
19. # retalhar a imagem em retangulos
20. sizex = 40
21. sizey = 40
22. sizez = 20
23. 
24. z = arange(0,img.shape[0],sizez)
25. x = arange(0,img.shape[1],sizex)
26. y = arange(0,img.shape[2],sizey)
27. print 'initial memory use:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
28. for k in arange(z.shape[0]-1):
29.     for i in arange(x.shape[0]-1):
30.         for j in arange(y.shape[0]-1):
31.             roi = np.zeros(img.shape, dtype = bool)
32.             roi[z[k]:z[k+1],x[i]:x[i+1],y[j]:y[j+1]] =1
33.             print 'memory used before:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
34.             print 'img.max',int(img.max()+1)
35.             glcm,v = glcmdesc(img,mask = roi)
36.             print 'memory used after:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
37. 
38. 
39. 
40. print 'final memory use:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 15, in <module>
NameError: name 'dicom' is not defined

------------------------------------------------------------