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 [1] to the right [0] for this window as documented here would then be

glcm = greycomatrix(glcm_window, [1], [0],  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[0] ):
    print i,
    for j in range(sarraster.shape[1] ):
        
        #windows needs to fit completely in image
        if i <3 or j <3:        
            continue
        if i > (contrastraster.shape[0] - 4) or j > (contrastraster.shape[0] - 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, [1], [0],  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)
sarraster = sarfile.ReadAsArray()


#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[0] ):
    print i,
    for j in range(sarraster.shape[1] ):
        
        # windows needs to fit completely in image
        if i > (contrastraster.shape[0] - 4) or j > (contrastraster.shape[0] - 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, [1], [0],  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()

7 comments:

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

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

      Delete
  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

    ReplyDelete
    Replies
    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[0], 10 )"

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

    ReplyDelete
  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 [1] [0]"""
    texture = np.zeros_like(img)

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

    # don't calculate at edges
    if (i < 3) or \
    (i > (img.shape[0])) or \
    (j < 3) or \
    (j > (img.shape[0] - 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, [1], [0], 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

    ReplyDelete
  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

    ReplyDelete