skimage.feature.greycomatrix only producing diagonal values

Question

I am attempting to produce glcm on a trend-reduced digital elevation model. My current problem is that the output of skimage.feature.greycomatrix(image) only contains values in the diagonal entries of the matrix.

glcm = greycomatrix(image,distances=[1],levels=100,angles=[0] ,symmetric=True,normed=True)

The image is quantized prior with the following code:

import numpy as np
from skimage.feature import greycomatrix
def quantize(raster):
    print("\n Quantizing \n")
    raster += (np.abs(np.min(raster)) + 1)
    mean = np.nanmean(raster.raster[raster.raster > 0])
    std = np.nanstd(raster.raster[raster.raster > 0])

    raster[raster == None] = 0 # set all None values to 0
    raster[np.isnan(raster)] = 0
    raster[raster > (mean + 1.5*std)] = 0
    raster[raster < (mean - 1.5*std)] = 0 # High pass filter
    raster[raster > 0] = raster[raster > 0] - (np.min(raster[raster > 0]) - 1)
    raster[raster>101] = 0
    raster = np.rint(raster)
    flat = np.ndarray.flatten(raster[raster > 0])
    range = np.max(flat) - np.min(flat)
    print("\n\nRaster Range: {}\n\n".format(range))
    raster = raster.astype(np.uint8)
    raster[raster > 101] = 0 

How would I go about making the glcm compute values outside of the diagonal matrix (i.e. just the frequencies of the values themselves), and is there something fundamentally wrong with my approach?

Answer 1

If pixel intensities are correlated in an image, the co-occurrence of two similar levels is highly probable, and therefore the nonzero elements of the corresponding GLCM will concentrate around the main diagonal. In contrast, if pixel intensities are uncorrelated the nonzero elements of the GLCM will be spread all over the matrix. The following example makes this apparent:

import numpy as np
from skimage import data
import matplotlib.pyplot as plt
from skimage.feature import greycomatrix

x = data.brick()
y = data.gravel()

mx = greycomatrix(x, distances=[1], levels=256, angles=[0], normed=True)
my = greycomatrix(y, distances=[1], levels=256, angles=[0], normed=True)

fig, ax = plt.subplots(2, 2, figsize=(12, 8))
ax[0, 0].imshow(x, cmap='gray')
ax[0, 1].imshow(mx[:, :, 0, 0])
ax[1, 0].imshow(y, cmap='gray')
ax[1, 1].imshow(my[:, :, 0, 0])

Despite I haven't seen your raster image I'm guessing that the intensity changes very smoothly across the image returned by quantize, and hence the GLCM is mostly diagonal.