back to lib

# Function glcmdesc

## Synopse

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

• Output
• v: array with GLCM descriptors
• glcm: GLCM matrix
• Input
• f: input image
• offset: defines the distance and orientation

## 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
 [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   ]
```

```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]
```

```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.]
```

```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.

## 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)):
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

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

```