## Wednesday, February 24, 2016

### Creating a Texture Image with a GLCM Co-Occurence Matrix using Scikit-Image and Python

Many clustering methods like the previously mentioned k-means cluster all the image values without considering any possible relation or patterns between neighbouring pixels. Texture is a measure considering this spatial relation --- an area with similar absolute values may visually be very different showing a different texture pattern (Source: Wikipedia)

The sci-kit image package for Python allows to calculate the Grey Level Co-Occurence Matrix (GCLM) and its texture measures as described on this page. However, this description, just as in the previous k-means example, leaves out the process of creating the texture image itself.

On this page you find a very good tutorial explaining GCLM and texture. A texture image is created by moving a window over you raster replacing the centre pixel covered by this window with the calculated texture measure of the moving window, as you see in this image taken from Myrka Hall-Beyers excellent tutorial:

A 7x7 window for any pixel would be defined by

```glcm_window = sarraster[i-3: i+4, j-3 : j+4]

```

And the GCLM for a one-pixel offset  to the right  for this window as documented here would then be

```glcm = greycomatrix(glcm_window, , ,  symmetric = True, normed = True )

```

Now you have to calculate this for every single pixel. The only way for making this specific calculation I can come up with is the slow process of looping through each pixel. Since some software implementations (like the SNAP SAR Toolbox) make this much faster, there must be some other way, so maybe a reader of this post has an idea. (I asked this question also at http://stackoverflow.com/questions/35551249/implementing-glcm-texture-feature-with-scikit-image-and-python ).

The following is only practical for small raster images, named sarraster in this case:

```for i in range(sarraster.shape ):
print i,
for j in range(sarraster.shape ):

#windows needs to fit completely in image
if i <3 or j <3:
continue
if i > (contrastraster.shape - 4) or j > (contrastraster.shape - 4):
continue

#Define size of moving window
glcm_window = sarraster[i-3: i+4, j-3 : j+4]
#Calculate GLCM and textures
glcm = greycomatrix(glcm_window, , ,  symmetric = True, normed = True )

# Calculate texture and write into raster where moving window is centered
contrastraster[i,j]      = greycoprops(glcm, 'contrast')
dissimilarityraster[i,j] = greycoprops(glcm, 'dissimilarity')
homogeneityraster[i,j]   = greycoprops(glcm, 'homogeneity')
energyraster[i,j]        = greycoprops(glcm, 'energy')
correlationraster[i,j]   = greycoprops(glcm, 'correlation')
ASMraster[i,j]           = greycoprops(glcm, 'ASM')
glcm = None
glcm_window = None

```

The following is the original image with various texture images. This is just a first test example, so no clear patterns are seen in this image:

The complete script for creating the image above. For a 1100 x 1400 raster this script takes about an hour while in the SNAP SAR Toolbox such calculations take < 1 minute, so I wonder about alternative implementations:

```# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import gdal, gdalconst
import numpy as np
from skimage.feature import greycomatrix, greycoprops

#Read SAR image into Numpy Array
filename = "//mnt//glaciology//SARimage.jpg"
sarfile = gdal.Open(filename, gdalconst.GA_ReadOnly)

#Create rasters to receive texture and define filenames
contrastraster = np.copy(sarraster)
contrastraster[:] = 0

dissimilarityraster = np.copy(sarraster)
dissimilarityraster[:] = 0

homogeneityraster = np.copy(sarraster)
homogeneityraster[:] = 0

energyraster = np.copy(sarraster)
energyraster[:] = 0

correlationraster = np.copy(sarraster)
correlationraster[:] = 0

ASMraster = np.copy(sarraster)
ASMraster[:] = 0

# Create figure to receive results
fig = plt.figure()
fig.suptitle('GLCM Textures')

# In first subplot add original SAR image
ax = plt.subplot(241)
plt.axis('off')
ax.set_title('Original Image')
plt.imshow(sarraster, cmap = 'gray')

for i in range(sarraster.shape ):
print i,
for j in range(sarraster.shape ):

# windows needs to fit completely in image
if i > (contrastraster.shape - 4) or j > (contrastraster.shape - 4):
continue

# Define size of moving window
glcm_window = sarraster[i-3: i+4, j-3 : j+4]
# Calculate GLCM and textures
glcm = greycomatrix(glcm_window, , ,  symmetric = True, normed = True )

# Calculate texture and write into raster where moving window is centered
contrastraster[i,j]      = greycoprops(glcm, 'contrast')
dissimilarityraster[i,j] = greycoprops(glcm, 'dissimilarity')
homogeneityraster[i,j]   = greycoprops(glcm, 'homogeneity')
energyraster[i,j]        = greycoprops(glcm, 'energy')
correlationraster[i,j]   = greycoprops(glcm, 'correlation')
ASMraster[i,j]           = greycoprops(glcm, 'ASM')
glcm = None
glcm_window = None

texturelist = {1: 'contrast', 2: 'dissimilarity', 3: ' homogeneity', 4: 'energy', 5: 'correlation', 6: 'ASM'}
for key in texturelist:
ax = plt.subplot(2,3,key)
plt.axis('off')
ax.set_title(texturelist[key])
plt.imshow(eval(texturelist[key] + "raster"), cmap = 'gray')

plt.show()
```

1. I have issue with this line 48;
if(contrastraster.shape - 4) or j > (contrastraster.shape - 4):
plz correct me, if i'm wrong

1. thanks for that .. you are indeed right, this should read
if i > (contrastraster.shape - 4) or j > (contrastraster.shape - 4):

2. Hi,
Very interseting, did you have any suggestion to do that with a specefic step to the pixels of the image ( i don't like to compute the glcm for all pixels but every 10 pixels for example).
Thank you for answering

1. you would add a step of ten to the range command in the loop, like here https://docs.python.org/2/library/functions.html#range

something like "for i in range(0, sarraster.shape, 10 )"

3. Probably not relevant anymore, but I believe you're doing everything right, except that pure Python is not really good for this sort of task. SNAP should be doing it in a similar way, except it uses compiled binaries, thus the speed difference.

One way to optimize it would be to use Numpy strides, where no subcopies of the arrays are created, and the data is read directly from the memory. Some good examples here: http://stackoverflow.com/questions/8070349/using-numpy-stride-tricks-to-get-non-overlapping-array-blocks

The other is Cython or F2PY. From my experience with F2PY moving windows you get speed increase x30-x300, but one would have to find a way how to call `greycomatrix` within the loop. Probably Cython would help more in this case.

4. num_level = 5
u = np.ogrid[0:num_level]
u = (1+ np.array(range(num_level)).reshape(num_level, 1))
v = 1+ np.array(range(num_level)).reshape((1, num_level))
y = np.zeros((num_level,), dtype=np.int)
XX, YY = np.meshgrid(u, y)
g = np.float_(XX + YY)

P=g.transpose()
XX, YY = np.meshgrid(v, y)
f = XX + YY
k=np.dot(f,g)
def glcm_image(img, measure="autocorrelation"):
"""TODO: allow different window sizes by parameterizing 3, 4. Also should
parameterize direction vector  """
texture = np.zeros_like(img)

# quadratic looping in python w/o vectorized routine, yuck!
for i in range(img.shape ):
for j in range(img.shape ):

# don't calculate at edges
if (i < 3) or \
(i > (img.shape)) or \
(j < 3) or \
(j > (img.shape - 4)):
continue

# calculate glcm matrix for 5 x 5 window, use dissimilarity (can swap in
# contrast, etc.)
glcm_window = img[i-2: i+3, j-2 : j+3]
glcm = np.float_(greycomatrix(glcm_window, , , symmetric = True, normed = True ))
glcm1_window = np.float_(glcm[i-2: i+3, j-2 : j+3])
print glcm[4,4].astype(float)
glcm1_window=np.reshape(glcm1_window, (5, 5))
texture[i,j] = np.float_(sum(sum(g.dot(glcm1_window))))

return texture

Sir,
I have an error in the above code. The glcm matrix that i get always contains zeros.I am not getting floating point values.Can you please tell me how i can correct
regards
srividhya

5. Sir, in the above code I am getting the error as
Traceback (most recent call last):
File "C:\Python27\ed.py", line 39, in
plt.imshow(sarraster, cmap = 'gray')
File "C:\Python27\lib\site-packages\matplotlib\pyplot.py", line 3080, in imshow
**kwargs)
File "C:\Python27\lib\site-packages\matplotlib\__init__.py", line 1710, in inner
return func(ax, *args, **kwargs)
File "C:\Python27\lib\site-packages\matplotlib\axes\_axes.py", line 5194, in imshow
im.set_data(X)
File "C:\Python27\lib\site-packages\matplotlib\image.py", line 604, in set_data
raise TypeError("Invalid dimensions for image data")
TypeError: Invalid dimensions for image data

I am using a jpg image as you have used.
Please suggest a way to resolve this error

1. maybe.. your image is not grayscale...
try input the gray image