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()

I have issue with this line 48;

ReplyDeleteif(contrastraster.shape[0] - 4) or j > (contrastraster.shape[0] - 4):

plz correct me, if i'm wrong

thanks for that .. you are indeed right, this should read

Deleteif i > (contrastraster.shape[0] - 4) or j > (contrastraster.shape[0] - 4):

Hi,

ReplyDeleteVery 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

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

Deletesomething like "for i in range(0, sarraster.shape[0], 10 )"

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.

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

num_level = 5

ReplyDeleteu = 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